A Finite Volume Scheme for the Solution of a Mixed Discrete-Continuous Fragmentation Model

TL;DR

This paper develops a finite volume numerical scheme for a mixed discrete-continuous fragmentation model, proving its convergence, conservation properties, and the uniqueness of solutions, supported by numerical simulations.

Contribution

It introduces a novel finite volume scheme for the mixed fragmentation equation and proves its convergence, conservation, and solution uniqueness.

Findings

Scheme conserves mass and nonnegativity

Weak solutions are unique and classical

Numerical simulations confirm convergence rate

Abstract

This paper concerns the construction and analysis of a numerical scheme for a mixed discrete-continuous fragmentation equation. A finite volume scheme is developed, based on a conservative formulation of a truncated version of the equations. The approximate solutions provided by this scheme are first shown to display conservation of mass and preservation of nonnegativity. Then, by utilising a Dunford--Pettis style argument, the sequence of approximate solutions generated is shown, under given restrictions on the model and the mesh, to converge (weakly) in an appropriate $L_1$ space to a weak solution to the problem. Additionally, by applying the methods and theory of operator semigroups, we are further able to show that weak solutions to the problem are unique and necessarily classical (differentiable) solutions. Finally, numerical simulations are performed to investigate the performance of the scheme and assess its rate of convergence.

Equations

$$\frac{\partial u_C(x,t)}{\partial t}$$

$$u_C(x,0)$$

$$\frac{\mathrm{d}{u}_{Di}(t)}{\mathrm{d}t}$$

$$u_D(0)$$

$${\int }_N^{y}xb(x|y)\mathrm{d}x+\sum _{j=1}^{N}j{b}_j(y)=y\text{for}y>N,$$

$$\sum _{j=1}^{i-1}j{b}_{j,i}=i\text{for}i=2,\dots ,N.$$

$$\frac{\partial u_C^R(x,t)}{\partial t}$$

$$u_C^R(x,0)$$

$$\frac{\partial \left(x{u}_C^R\right)}{\partial t}=\frac{\partial {\mathcal{F}}^R\left(x{u}_C^R\right)}{\partial x}-S(x{u}_C^R),u_C^R(x,0)=c_0(x),\text{for}N<x<R,t>0,$$

$${\mathcal{F}}^R(f)={\int }_x^R{\int }_N^x\frac{y}{z}a(z)b(y|z)f(z)\mathrm{d}y\mathrm{d}z,S(f)=\frac{a(x)}{x}\sum _{i=1}^{N}i{b}_i(x)f(x),\text{for}N<x<R.$$

$$\underset{x\to N,R}{\lim }\left|{\mathcal{F}}^R(f)(x)\right|=0.$$

$$\left|{\mathcal{F}}^R(f)(x)\right|\le {\int }_N^R{\chi }_{(x,R)}(z)\frac{a(z)\left|f(z)\right|}{z}\left({\int }_N^{x}yb(y|z)\mathrm{d}y\right)\mathrm{d}z,$$

$${\chi }_{(x,R)}(z)\left({\int }_N^{x}yb(y|z)\mathrm{d}y\right)\le {\int }_N^{z}yb(y|z)\mathrm{d}y\le z,$$

$${\chi }_{(x,R)}(z)\left({\int }_N^{x}yb(y|z)\mathrm{d}y\right)\le x\beta (R)(x-N).$$

$${\chi }_{(x,R)}(z)\left({\int }_N^{x}yb(y|z)\mathrm{d}y\right)\le {\chi }_{(x,R)}(z){\int }_N^{z}yb(y|z)\mathrm{d}y\le {\chi }_{(x,R)}(z)z,$$

$$\frac{\mathrm{d}{u}_{Di}^R(t)}{\mathrm{d}t}$$

$$u_{Di}^R(0)$$

$$\Vert f{\Vert }_{L_1}={\int }_0^T{\int }_N^R\left|f(x,t)\right|\mathrm{d}x\mathrm{d}t\le \frac{1}{N}{\int }_0^T{\int }_N^R\left|f(x,t)\right|x\mathrm{d}x\mathrm{d}t=\frac{1}{N}\Vert f{\Vert }_{L_1^1}.$$

$$\Vert f{\Vert }_{L_1^1}={\int }_0^T{\int }_N^R\left|f(x,t)\right|x\mathrm{d}x\mathrm{d}t\le R{\int }_0^T{\int }_N^R\left|f(x,t)\right|\mathrm{d}x\mathrm{d}t=R\Vert f{\Vert }_{L_1}.$$

$$\frac{\Phi (x)}{x}\to \infty \text{as}x\to \infty \text{and}\underset{f\in \mathcal{F}}{\sup }{\int }_{\Omega }\Phi \left(\left|f\right|\right)\mathrm{d}\mu <\infty .$$

$$\Phi (x)\ge \Phi (y)+{\Phi }^{\prime }(y)(x-y),$$

$$\underset{=0}{\underbrace{\Phi (0)}}\ge \underset{\ge 0}{\underbrace{\Phi (y)}}+{\Phi }^{\prime }(y)(0-y).$$

$$y{\Phi }^{\prime }(y)\ge \Phi (y)\ge 0.$$

$$\underset{f\in \mathcal{F}}{\sup }{\int }_{\Omega \setminus A_{\epsilon }}\left|f\right|\mathrm{d}\mu \le \epsilon .$$

$$x_{-1/2}=N,x_{I_h-1/2}=R,x_i=(x_{i-1/2}+x_{i+1/2})/2,h/k<\Delta x_i=x_{i+1/2}-x_{i-1/2}<h,$$

$$k_{1}h\le \Delta t\le k_{2}h.$$

$$A_i=\frac{1}{\Delta x_i}{\int }_{{\Lambda }_i}a(x)\mathrm{d}x\text{for}i=0,1,\dots ,I_h-1,$$

$$B_{i,j}=\frac{1}{\Delta x_i\Delta x_j}{\int }_{{\Lambda }_j}{\int }_{{\Lambda }_i}b(x|y)\mathrm{d}x\mathrm{d}y\text{for}i=0,1,...,I_h-1\text{and}j=0,1,\dots ,I_h-1,$$

$${\tilde{B}}_{i,j}=\frac{1}{\Delta x_j}{\int }_{{\Lambda }_j}{b}_i(y)\mathrm{d}y\text{for}i=1,2,\dots ,N\text{and}j=0,1,\dots ,I_h-1.$$

$$a^h(x)=\sum _{i=0}^{I_h-1}{\chi }_{{\Lambda }_i}(x)A_i,b^h(x|y)=\sum _{j=0}^{I_h-1}\sum _{i=0}^{I_h-1}{\chi }_{{\Lambda }_i}(x){\chi }_{{\Lambda }_j}(y)B_{i,j},b_i^h(y)=\sum _{j=0}^{I_h-1}{\chi }_{{\Lambda }_j}(y){\tilde{B}}_{i,j}.$$

$$x_i\frac{u_C^{n+1,i}-u_C^{n,i}}{\Delta t}=\frac{F_{i+1/2}^n-F_{i-1/2}^n}{\Delta x_i}-S_i^n,$$

$$\left({\mathcal{F}}^R(x{u}_C^R)\right)(x_{i-1/2})$$

$$=\sum _{j=i}^{I_h-1}{\int }_{{\Lambda }_j}\left(\sum _{k=0}^{i-1}{\int }_{{\Lambda }_k}ya(z)b(y|z)u_C^R(z,t)\mathrm{d}y\right)\mathrm{d}z$$

Full text

1. Introduction

Fragmentation and coagulation processes occur in many physical systems, with the associated mathematical models receiving much attention in the literature. Example application areas include colloid science [1, 2], population dynamics [3, 4] and astrophysics [5, 6]. However, analytical solutions to these models are only available for a limited number of specific cases, and we often have to resort to approximate solutions generated by an appropriate numerical scheme. A range of numerical techniques have been applied to these problems, and these broadly fall into two categories: those involving a stochastic (Monte Carlo) element, for example [7, 8, 9] and those based around various deterministic approximation schemes [10, 11, 12, 13]. The introductory chapter of [14] and the references therein provide a detailed overview of a number of these approximation methods.

In the earlier work [15], we presented a mixed discrete-continuous model of fragmentation in an attempt to resolve the issue of ‘shattering’ mass-loss observed in some purely continuous models [16]. By modelling the mass distribution amongst the smallest particles using a discrete model, whilst modelling the distribution of larger particle masses with a continuous model, the aim was to introduce a higher degree of physical fidelity thus resolving the shattering mass-loss problem, whilst also retaining the mathematical efficiency of the continuous model.

Given the similarities between this model and those existing in the literature, in addition to the added complexity of the mixed framework, we would expect in most cases to have to rely on numerical methods to obtain a solution. In this paper we present a numerical scheme for the solution of the mixed fragmentation model proposed in [15]. The basis of the scheme is a finite volume discretisation of the continuous regime equation. The use of such a method would appear a reasonable choice in this case, given its conservative nature and the motivation behind the model development. Indeed, finite volume schemes have been commonly applied to the solution of coagulation and fragmentation equations, with the first such use being [17], where the case of pure coagulation was considered. For problems involving fragmentation, the article [18] sees such a scheme employed in approximating the binary coagulation and fragmentation equation, whilst [14] and [19] examine their use for the multiple fragmentation equation, with [20] extending this to include coagulation. Further works have seen these methods applied to a number of coagulation–fragmentation model variants, for example with the inclusion of spacial diffusion [21] and additional nucleation and growth processes [22, 23]. Whilst a number of articles [24, 25, 26] cover the approximation of multi-dimensional coagulation or fragmentation, whereby particles may be classified by additional variables beyond their mass or volume.

1.1. Mixed Discrete-Continuous Model

In the mixed model of [15], a cut-off value $N\in\mathbb{N}$ is introduced; above this cut-off, particle mass is considered as a continuous variable, whilst below it, the particles are forced to take discrete integer masses. If we denote by $u_{C}(x,t)$ the particle mass density within the continuous mass regime ($x>N$), then the evolution of $u_{C}(x,t)$ is governed by the continuous multiple fragmentation equation:

[TABLE]

This equation is similar in form to the multiple fragmentation equation introduced in [16]. The function $a(x)$ provides the fragmentation rate for a particle of mass $x$, whilst $b(x|y)$ represents the distribution of particles of mass $x>N$ resulting from the break-up of a particle of mass $y>x$. The functions $a$ and $b$ are assumed to be nonnegative measurable functions, defined on $\left(N,\infty\right)$ and $\left(N,\infty\right)\times\left(N,\infty\right)$, respectively. We also require $b(x|y)=0$ for $x>y$, since no particle resulting from a fragmentation event can have a mass exceeding the original particle. The initial mass distribution within the continuous regime is given by the nonnegative function $c_{0}(x)$.

Letting $u_{Di}(t)$ denote the concentration of discrete mass $i$-mer particles ($i\leq N$) and $u_{D}(t)$ the $N$-component vector taking these values as entries, the change in the values $u_{Di}(t)$, $i=1,\dots,N$, is governed by the equation:

[TABLE]

In the case of $i=N$, the second term becomes an empty sum and is taken to be [math]. The values $a_{i}$ give the rates at which $i$-mer particles fragment, with $a_{1}=0$. The quantities $b_{i,j}$ give the expected number of $i$-mers produced from the fragmentation of a $j$-mer and the functions $b_{i}(y)$ give the expected number of $i$-mers produced from the fragmentation of a particle of mass $y>N$. The underlying physics demands that each $a_{i}$, $b_{i,j}$ and $b_{i}(y)$ be nonnegative. Finally, $d_{0}$ is the $N$-component vector of nonnegative values, specifying the initial concentrations within the discrete regime.

During each fragmentation event, mass is simply redistributed from the larger particle to the smaller resulting particles, but the total mass involved should be conserved. This gives rise to the following two conditions to supplement equations (1.1) and (1.2):

[TABLE]

The condition (1.3) is an expression of mass conservation upon the fragmentation of a particle from the continuous mass regime. The equation (1.4) comes from the conservation of mass when a particle from the discrete mass regime breaks up.

For further details on the mixed discrete-continuous model and its properties, the reader is directed to consult [15] or [27].

1.2. Truncation and Reformulation

When considering the numerical solution of equations (1.1) and (1.2), we encounter an issue in that the range of the continuous mass variable $x$ is an unbounded interval, which presents a computational problem. We therefore introduce a truncation parameter $R>N$, and restrict the continuous mass variable to the range $N<x<R$. Therefore, in place of equations (1.1), we consider the truncated version

[TABLE]

where $\chi_{(N,R)}$ denotes the characteristic function of the interval $(N,R)$. Taking our lead from the aforementioned articles, we now rewrite equation (1.5) in a conservative form, although in our case we must include an additional sink term to account for the mass leaked down to the discrete regime. Therefore we end up with the following equation for the mass quantity $xu_{C}^{R}(x,t)$:

[TABLE]

where $\mathcal{F}^{R}$ and $S$ are a truncated flux term and sink term, respectively, given by

[TABLE]

The equation (1.5) may be recovered from (1.6) by a formal application of Leibniz’s rule for differentiating under the integral. However, the equivalence of the two forms can be seen to be justified rigorously in [27, Appendix C]. Before continuing, we establish a result concerning the behaviour of the flux term $\mathcal{F}^{R}$ at the limits of our domain.

Lemma 1.1**.**

If the kernels $a$ and $b$ are assumed to belong to $L_{\infty,loc}$ on the domains $[N,\infty)$ and $[N,\infty)\times[N,\infty)$ respectively, which will be the case in the upcoming analysis, then for $f\in L_{1}(N,R)$ the flux term $\mathcal{F}^{R}(f)$ satisfies

[TABLE]

Proof.

It is a straightforward matter to bound $\mathcal{F}^{R}(f)$ as follows:

[TABLE]

which holds for $x\in(N,R)$. Recalling the mass conservation condition (1.3), we deduce that

[TABLE]

for all $z\in(N,R)$. Hence the integrand appearing in (1.7) is bounded above by $a(z)\left|f(z)\right|$, which, thanks to $a\in L_{\infty,loc}[N,\infty)$ and $f\in L_{1}(N,R)$, is integrable.

Considering the limit as $x\rightarrow N$ first, if we denote by $\beta(R)$ the essential supremum of $b$ over $[N,R]\times[N,R]$, then we have

[TABLE]

As such, the integrand in (1.7) converges pointwise to 0 over $z\in(N,R)$ as we let $x\searrow N$. An application of the Lebesgue dominated convergence theorem then gives the required convergence of $\left|\mathcal{F}^{R}(f)(x)\right|$ as $x\searrow N$. Turning now to the limit as $x\nearrow R$, another application of condition (1.3) provides us with

[TABLE]

for $z\in(N,R)$. Therefore, the integrand from (1.7) must again converge pointwise to 0 over $(N,R)$, this time as we let $x\nearrow R$. Another application of the Lebesgue dominated convergence theorem gives the convergence of $\left|\mathcal{F}^{R}(f)(x)\right|$ to 0, as $x\nearrow R$. ∎

This result will be utilised later in a number of arguments, most significantly in approximating $\mathcal{F}^{R}$ within our numerical scheme and in establishing a weak formulation of equation (1.6).

The truncation of the continuous mass interval also has an impact on our discrete regime equation; therefore, instead of equation (1.2), we consider

[TABLE]

In the case of $i=N$, the empty sum above is taken to be zero; this convention will be adopted in all similar cases which follow.

This truncation procedure is a standard approach when dealing with fragmentation and coagulation problems, having been applied for example in [28], where the theory and methods of operator semigroups were employed, and [29] where an alternative weak compactness style argument was adopted. The common approach of these works involves establishing the existence of solutions to a sequence of such truncated problems. A limit is then obtained as the truncation point is increased without bound, with this limit then being shown to satisfy the untruncated problem in some sense. Although in this article we restrict our attention to the numerical approximation of the truncated discrete–continuous problem, as given by equations (1.5) and (1.8), it can be shown that the solutions to the truncated problems converge, in an appropriate space, to give the solutions to the untruncated (1.1) and (1.2). The proof of this convergence argument follows similar lines to that set out in [30, Section 8.3.2], with the reader being directed to [27, Chapter 6] for the specific details. Additionally, the reader may find an empirical examination of this convergence in [27, Section 7.5], where the key factors influencing the convergence, and thus the selection of an appropriate $R$ are investigated.

2. Preliminaries

Having set out our problem in the previous section, we now present a brief outline of the key results which appear in the upcoming material and which may be considered nonstandard or which are particular to our case.

Theorem 2.1**.**

In the analysis pursued in subsequent results, we shall be working extensively in spaces of the type $L_{1}$. In particular we shall be working in the spaces $L_{1}=L_{1}((N,R)\times[0,T),\,{\mathrm{d}}x\,{\mathrm{d}}t)$ and $L_{1}^{1}=L_{1}((N,R)\times[0,T),x\,{\mathrm{d}}x\,{\mathrm{d}}t)$, where $N$ is a positive integer and $R>N$ is a finite real value. With the associated norms, these form equivalent spaces.

Proof.

First let us suppose that $f\in L_{1}^{1}$; then we have

[TABLE]

Therefore $f\in L_{1}$ also, with $\|f\|_{L_{1}}\leq\frac{1}{N}\|f\|_{L_{1}^{1}}$. Now let us assume that $f\in L_{1}$; then we have

[TABLE]

Hence $f\in L_{1}^{1}$ with $\|f\|_{L_{1}^{1}}\leq R\|f\|_{L_{1}}$. Taken together, the above results show us that the spaces $L_{1}((N,R)\times[0,T),\,{\mathrm{d}}x\,{\mathrm{d}}t)$ and $L_{1}((N,R)\times[0,t),x\,{\mathrm{d}}x\,{\mathrm{d}}t)$ contain the same elements and have equivalent norms. ∎

This result shall prove useful in the forthcoming analysis, allowing us to switch spaces when mathematically convenient whilst retaining convergence.

Given a sequence $\left\{f_{n}\right\}_{n=1}^{\infty}$ in a normed vector space $\left(X,\|\cdot\|\right)$, we assume the reader is familiar with the concept of weak convergence and in particular its definition in spaces of the type $L_{1}(\Omega,\mu)$. In our analysis when handling weakly convergent sequences we will usually find them appearing alongside other factors and we would like the product to converge weakly also. The following theorem gives us sufficient conditions for the product of two sequences to converge weakly and will be used extensively in the convergence proofs for our numerical schemes.

Theorem 2.2**.**

Let $\left(\Omega,\mathcal{A},\mu\right)$ be a measure space with $\mu$ finite. Suppose $f_{h}\rightharpoonup f$ in $L_{1}\left(\Omega,\mu\right)$, $g_{h}\rightarrow g$ point-wise $\mu$ a.e. in $\Omega$, and $\sup_{h}\|g_{h}\|_{L_{\infty}}<\infty$, then $f_{h}g_{h}\rightharpoonup fg$ in $L_{1}\left(\Omega,\mu\right)$.

Proof.

The reader is referred to [31, Proposition 2.61]. ∎

The main part of our convergence argument utilises the Dunford–Pettis theorem, which provides us with sufficient conditions to establish the weak convergence of our sequence of approximations. One such condition is that of equiintegrability. There are a number of equivalent characterisations of equiintegrability, which the reader may find in [31, Theorem 2.29]. For our purposes the most important characterisation of equiintegrability is given by de la Vallée Poussin’s theorem, a refined version of which is given below.

Theorem 2.3**.**

(de la Vallée Poussin’s Theorem)* Let $\mathcal{F}$ be a bounded subset of $L_{1}\left(\Omega,\mu\right)$, then $\mathcal{F}$ is equiintegrable if and only if there exists a nonnegative, convex function $\Phi\in C^{\infty}\left([0,\infty)\right)$, with $\Phi(0)=0$ and $\Phi^{\prime}(0)=1$, such that $\Phi^{\prime}$ is concave and*

[TABLE]

Proof.

The necessity of this condition can be derived easily from [32, Theorem 8], which under the assumption that $\mathcal{F}$ is equiintegrable provides us with a $\Psi$ satisfying all the stated conditions with the exception that the function $\Psi$ has derivative [math] at [math] and is not stated to be nonnegative. Given such a $\Psi$, we set $\Phi(x)=\Psi(x)+x$. Then $\Phi$ retains the required properties of $\Psi$ but additionally $\Phi^{\prime}(0)=1$. Also, by utilising the following standard inequality for $C^{1}$ convex functions

[TABLE]

with $x\geq 0$ and $y=0$ we can see that $\Phi(x)$ must be nonnegative on $[0,\infty)$. The sufficiency of our conditions comes straight from the standard version of the de la Vallée Poussin theorem [31, Theorem 2.29 (iii)]. ∎

In our analysis we shall require some properties of such a function, which we set out in the following lemma.

Lemma 2.4**.**

Let $\Phi$ be as in Theorem 2.3; then for nonnegative $x$ and $y$ we have the following:

1. (1)

$x\Phi^{\prime}(y)\leq\Phi(x)+\Phi(y)$, 2. (2)

$\Phi^{\prime}(y)\geq 0$.

Proof.

The first of these inequalities is nonstandard and the proof can be found in [32, Proposition 13 (30)]. For the second property we return to inequality (2.1), with $x=0$ and $y\geq 0$, which gives us

[TABLE]

An obvious rearrangement yields

[TABLE]

Now in the case that $y=0$ property (ii) is given by the definition of $\Phi$. Hence we may assume that $y>0$ and divide through by it to obtain the desired result that $\Phi^{\prime}(y)\geq 0$. ∎

We now come to the Dunford–Pettis theorem, one of the most significant technical tools applied in this work. The theorem provides necessary and sufficient conditions for a subset of an $L_{1}$ space to be weakly sequentially compact. That is, any sequence in the subset must have a subsequence which is weakly convergent.

Theorem 2.5**.**

(Dunford–Pettis Theorem)* Let $\left(\Omega,\mathcal{A},\mu\right)$ be a measure space and let $\mathcal{F}\subset L_{1}\left(\Omega,\mu\right)$. Then $\mathcal{F}$ is weakly sequentially compact if and only if the following conditions are satisfied:*

1. (1)

$\mathcal{F}$* is bounded in $L_{1}\left(\Omega,\mu\right)$;* 2. (2)

$\mathcal{F}$* is equiintegrable;* 3. (3)

For every $\varepsilon>0$ there exists $A_{\varepsilon}\subset\Omega$ with $A_{\varepsilon}\in\mathcal{A}$ such that $\mu\left(A_{\varepsilon}\right)<\infty$ and

[TABLE]

We note that in the case that $\mu(\Omega)<\infty$ condition (iii) is automatically satisfied by taking $A_{\varepsilon}=\Omega$ for all values of $\varepsilon$.

Proof.

See [31, Theorem 2.54]. ∎

In the later analysis of this paper we shall be relying heavily on the methods and theory of operator semigroups. In particular the concept of substochastic semigroups, the Kato–Voigt perturbation theorem and the notion of semigroup honesty. For the sake of brevity we refrain from outlining such material here, however the reader may find details of the requisite results in the preliminary sections of [15] or [27] or the text [30].

3. Development of the Numerical Scheme

We now introduce our numerical approximation scheme for the truncated system, (1.6) and (1.8). First we must discretise the continuous mass variable $x$, and so we introduce the mesh $\left\{x_{i-1/2}\right\}_{i=0}^{I_{h}}$ on the interval $(N,R)$, with

[TABLE]

where $h\in(0,1)$ and $k>1$ is some constant. Additionally we denote the interval $[x_{i-1/2},x_{i+1/2})$ by $\Lambda_{i}$, however the (left-hand-most) interval $\Lambda_{0}$ is taken to be $(x_{-1/2},x_{1/2})$.

For the time variable $t$, if $T$ is the final time up to which we wish to compute an approximate solution, then we define the time step $\Delta t=T/M$ where $M$ is some large integer. The time points are then given by $t_{n}=n\Delta t$ for $n=0,1,\dots,M$ with corresponding time intervals $\tau_{n}=[t_{n},t_{n+1})$ for $n=0,1,\dots,M-1$.

We restrict the choice of the mesh by assuming the existence of positive constants $k_{1}$ and $k_{2}$ so that the mesh sizes $h$ and $\Delta t$ satisfy

[TABLE]

The numerical scheme requires representative values for the functions $a(x)$, $b(x|y)$ and $b_{i}(y)$ over the appropriate intervals. This is done by taking their average value over each interval. Therefore we define

[TABLE]

as our approximation of $a(x)$ over the interval $\Lambda_{i}$. We approximate $b(x|y)$ over $\Lambda_{i}\times\Lambda_{j}$ by

[TABLE]

and the functions $b_{i}(y)$ are approximated over $\Lambda_{j}$ by the values

[TABLE]

We note by our initial assumption regarding the nonnegativity of $a$, $b$ and $b_{i}$, that each of the values introduced above must be nonnegative. If $\chi_{I}$ denotes the characteristic function of a set $I$, then we can construct piecewise constant approximations to the functions $a$, $b$ and $b_{i}$ as follows:

[TABLE]

Remark 3.1**.**

This is a standard means of approximation and assuming the choice of kernels is suitably restricted, the approximations will converge pointwise to the desired functions almost everywhere on the appropriate domains. In our case, the kernels $a$ and $b$ will be assumed to be $L_{\infty,loc}$ on $[N,\infty)$ and $[N,\infty)\times[N,\infty)$, respectively. In addition, the restriction (1.3) determines each $b_{i}$ as an element of $L_{\infty,loc}[N,\infty)$. Having $a$, $b$ and $b_{i}$ as $L_{\infty,loc}$ functions is sufficient to ensure that the approximations $a^{h}$, $b^{h}$ and $b_{i}^{h}$ converge pointwise to $a$, $b$ and $b_{i}$ almost everywhere in their respective domains. This is a standard result, however full details can be found in [27, Lemma 4.2.1].

We are now ready to construct the approximation scheme. Let $u_{C}^{n,i}$ denote our approximation to $u_{C}^{R}(x,t)$ over the mass interval $\Lambda_{i}$ for the time interval $\tau_{n}$. The equation (1.6) is then approximated by

[TABLE]

where $F_{i-1/2}^{n}$ is an approximation of the flux $\mathcal{F}^{R}(xu_{C}^{R})$ at the point $x=x_{i-1/2}$ over the time interval $\tau_{n}$, and is given by

[TABLE]

for $i=1,\ldots,I_{h}-1$, with $F_{-1/2}^{n}=F_{I_{h}-1/2}^{n}=0$, which can be justified by Lemma 1.1. The values $S_{i}^{n}$ approximate the sink term $S(xu_{C}^{R})$ over $\Lambda_{i}$ for the time interval $\tau_{n}$, and are computed by

[TABLE]

This gives rise to the following numerical method for the computation of the approximations $u_{C}^{n,i}$:

[TABLE]

The sequence of approximations generated by (3.3) requires us to provide an initial set of values to get started. For our starting values we simply average the initial datum over each of the mass intervals; hence

[TABLE]

Then our approximation to $u_{C}^{R}(x,t)$ over $(N,R)\times[0,T)$ is constructed as follows:

[TABLE]

Remark 3.2**.**

*The convergence proof for our numerical scheme requires the initial approximation given by (3.4) and (3.5) to converge strongly in $L_{1}(N,R)$ to the restriction of $c_{0}$ to $(N,R)$. Again, that this is the case with our definition of the discretised initial datum is a standard result and details can be found in [27, Lemma 4.2.2].

Now considering the discrete regime, let $u_{D}^{n,i}$ denote our approximation of $u_{Di}^{R}(t)$ over the time interval $\tau_{n}$. Equation (1.8) is then approximated as

[TABLE]

giving rise to the relation

[TABLE]

The initial values for the discrete approximation are simply given by the initial condition vector $d_{0}$, so that $u_{D}^{0,i}={d_{0}}_{i}$ for $i=1,\ldots,N$. Then our approximations $u_{Di}^{h}(t)$ to $u_{Di}^{R}(t)$ for $t\in[0,T)$ are given by

[TABLE]

4. Properties of Numerical Solutions: Nonnegativity and Mass Conservation

In the article [15] we proved the existence and uniqueness of a solution to our system (1.1) and (1.2). This solution was shown to possess a number of properties that we would expect given the physical nature of the model, namely the solution preserved nonnegativity and conserved total mass. In the following sections we examine whether the approximate solution provided by (3.5) and (3.7), also displays these properties. These properties, apart from being physically relevant, will also be utilised in the forthcoming proofs of the convergence of the approximations (3.5) and (3.7) to a solution to the system (1.5) and (1.8), and subsequently the uniqueness and differentiability of that solution.

4.1. Nonnegativity of the Numerical Solution

Lemma 4.1**.**

For a fixed partition $(x_{i-1/2})_{i=0}^{I_{h}}$, suppose that $\Delta t$ is sufficiently small that the following condition is satisfied:

[TABLE]

for all $i\in\left\{0,1,\dots,I_{h}-1\right\}$ such that the denominator is nonzero, and

[TABLE]

for all $i\in\left\{2,\dots,N\right\}$ such that $a_{i}\neq 0$. Then, the approximate solutions defined by (3.5) and (3.7) preserve nonnegativity.

Proof.

Starting with the approximation for the continuous regime, let us consider equation (3.3). By cancelling common terms we get that

[TABLE]

for $i=1,\dots,I_{h}-2$. Therefore we have

[TABLE]

Substituting this into (3.3) gives us

[TABLE]

for $i=1,\dots,I_{h}-2$. The cases $i=0$ and $i=I_{h}-1$ can be handled similarly to obtain the same result, where the empty sums are taken as 0.

From this it is clear that if each of the approximations $u_{C}^{n,i}$ is nonnegative, and provided $\Delta t$ is sufficiently small such that the term within the outer brackets is nonnegative, then each of the approximations $u_{C}^{n+1,i}$, for the subsequent time step, will also be nonnegative. Hence to ensure the approximations $u_{C}^{n+1,i}$ are nonnegative we can take

[TABLE]

In the case of the above denominator being zero for some $i$, such that the bound (4.3) is undefined, then $u_{C}^{n+1,i}$ can be seen from (4.2) to automatically satisfy the nonnegativity requirement, for any value of $\Delta t$.

Turning to the approximation for the discrete regime, it is immediately clear from the form of (3.6) that if all of the values $u_{C}^{n,i}$ and $u_{D}^{n,i}$ are nonnegative, then each $u_{D}^{n+1,i}$ will be nonnegative if for each $i=1,\dots,N$ we have that $1-\Delta ta_{i}$ is nonnegative. This can be ensured by taking

[TABLE]

Therefore if we choose a $\Delta t$ small enough that both (4.3) and (4.4) are satisfied, then our approximate solutions will remain nonnegative. ∎

From now on we shall assume that conditions (4.3) and (4.4) are satisfied and that $c_{0}(x)\geq 0$ and each $d_{0,i}\geq 0$ so that our approximations remain nonnegative.

Remark 4.2**.**

The bound (4.3) is dependent on the mesh and it is perhaps not immediately apparent how this bounding value might vary as we refine the mesh. In particular, it would be advantageous to confirm that it is indeed possible to find a constant $k_{1}$, such that conditions (3.1) and (4.3) can be satisfied simultaneously, whilst $h\searrow 0$. In the upcoming analysis we will place restrictions on the functions $a$ and $b$; these constraints will allow us to guarantee the existence of such a $k_{1}$.

The upcoming Theorem 5.3 imposes the restriction $a,b\in L_{\infty}$ on the restricted domains $[N,R]$ and $[N,R]\times[N,R]$ respectively, with $\alpha(R)$ and $\beta(R)$ being the essential suprema for $a$ and $b$ on said domains. This being the case, we have $A_{i}\leq\alpha(R)$ and $B_{k,i}\leq\beta(R)$ for all values of $i$ and $k$ admissible in (4.3). Furthermore, from (1.3) we may deduce that each $b_{i}(y)\leq y$, hence $\tilde{B}_{j,i}\leq R$. Finally, all mesh midpoints $x_{i}$ must clearly satisfy $x_{i}\geq N\geq 1>h$. Taken together, these bounds lead, via a simple calculation, to

[TABLE]

for $i=0,1,\dots,I_{h}-1$. Hence, we have established a possible value for $k_{1}$, which ensures (3.1) and (4.3) can be satisfied simultaneously as $h\searrow 0$.

4.2. Mass Conservation by the Numerical Solutions

In [15, Lemma 6.2], the exact solution to our system of equations (1.1) and (1.2) was shown to conserve mass between the two regimes. We now show that this property is shared by our numerical solutions.

Lemma 4.3**.**

The approximate solutions generated by (3.3) and (3.6) conserve mass.

Proof.

The mass associated with the approximate continuous regime solution, $u_{C}^{h}(x,t)$, is given by

[TABLE]

whilst the approximate solution $u_{D}^{h}(t)$ has associated mass given by

[TABLE]

Summing these two expressions gives the total mass:

[TABLE]

First let us examine the mass accounted for by the continuous regime. From the relation (3.3) we get

[TABLE]

The middle summation term is lost in going to the final line as the sum is telescoping with zero end terms. Now we consider the discrete regime mass; the generating relation (3.6) gives us

[TABLE]

The middle two terms cancel due to the mass conservation condition (1.4). Combining equations (4.2) and (4.2) we obtain

[TABLE]

From repeated application of this equality it is easily seen that the bracketed expression appearing in (4.6) is equal for all values of $n$, and hence the total mass $M^{h}(t)$ remains constant. ∎

5. Convergence of the Scheme to a Weak Solution as $h\rightarrow 0$

Having determined the nonnegativity and mass conservative properties of the approximate solutions provided by (3.5) and (3.7), in this section we set out to prove that they converge, in some sense, to a limit as the parameter $h$, and by necessity $\Delta t$, go to zero, and show that this limit itself is an ‘exact’ solution to our truncated model.

5.1. Continuous Fragmentation Regime: Convergence

Let us start with the continuous regime approximations $\left\{u^{h}_{C}\right\}$. In order to prove the (weak) convergence of this family, we employ a weak compactness argument, utilising the Dunford–Pettis theorem (Theorem 2.5), which provides necessary and sufficient conditions for weak compactness in an $L_{1}$ space. We begin by proving the equiboundedness of the set $\left\{u^{h}_{C}\right\}$.

Lemma 5.1**.**

The family of approximations $\left\{u_{C}^{h}\right\}$ is equibounded (uniformly bounded) in the space $L_{1}((N,R)\times[0,T),x\,{\mathrm{d}}x\,{\mathrm{d}}t)$.

Proof.

Recalling equation (4.2), we have for any $t\in[0,T)$ that

[TABLE]

From the analysis of Lemma 4.1, each of the values $u_{C}^{n,i}$ is nonnegative, and as such the values $S_{i}^{n}$ are nonnegative. Therefore, from the last line of equation (4.2) we deduce that

[TABLE]

Repeated application of this inequality yields

[TABLE]

The quantity $\frac{x_{i}}{x_{i-1/2}}$ can be bounded as follows:

[TABLE]

We note this bound as it will appear regularly in subsequent calculations. Substituting this within (5.1) yields

[TABLE]

for $n=0,\ldots,M-1$. Replacing this inequality in our calculation gives us the following, which holds for all $t\in[0,T)$:

[TABLE]

Integrating this inequality with respect to $t$ from [math] to $T$ we obtain the required equiboundedness of $\left\{u_{C}^{h}\right\}$ in the space $L_{1}((N,R)\times[0,T),x\,{\mathrm{d}}x\,{\mathrm{d}}t)$. ∎

We now move on to prove the second of the two required conditions for the Dunford–Pettis theorem, namely equiintegrability. However, prior to this we collect all the conditions so far imposed on our model via the functions $a$, $b$ and $b_{i}$, and the initial distributions $c_{0}$ and $d_{0}$, and also on our mesh via the parameters $h$ and $\Delta t$.

Remark 5.2**.**

Throughout the remaining analysis, the following conditions shall be assumed to be satisfied.

1. (1)

The initial mass distributions within the continuous and discrete regimes are nonnegative, that is $c_{0}(x)$ for all $x>N$ and $d_{0,i}\geq 0$ for $i=1,\ldots,N$. 2. (2)

*The functions $a$ and $b$ are nonnegative and belong to the spaces $L_{\infty,loc}\left([N,\infty)\right)$ and *

$L_{\infty,loc}\left([N,\infty)\times[N,\infty)\right)$, respectively. 3. (3)

Each of the functions $b_{i}$ is assumed to be nonnegative. This nonnegativity in conjunction with condition $2$ is sufficient to guarantee that $b_{i}\in L_{\infty,loc}\left([N,\infty)\right)$, as per Remark 3.1. 4. (4)

There exist positive constants $k_{1}$ and $k_{2}$ such that the mesh parameters $h$ and $\Delta t$ satisfy

[TABLE] 5. (5)

To ensure that the approximate solutions remain nonnegative, the time step $\Delta t$ is assumed to satify the following constraints:

[TABLE]

for all cases of the denominator being nonzero, and

[TABLE] 6. (6)

There exists a constant $\theta>0$ such that

[TABLE]

where $K(R)=\alpha(R)\beta(R)R$, with $\alpha(R)$ and $\beta(R)$ being the essential suprema for $a$ and $b$ on the restricted domains $[N,R]$ and $[N,R]\times[N,R]$, respectively.

Theorem 5.3**.**

Under the assumptions outlined in Remark 5.2, the family $\left\{u_{C}^{h}\right\}$ is equiintegrable in $L_{1}((N,R)\times[0,T),x\,{\mathrm{d}}x\,{\mathrm{d}}t)$.

Proof.

Consider the constant sequence comprising solely of the initial data $c_{0}\in L_{1}\left((N,R),x\,{\mathrm{d}}x\right)$. Clearly this sequence is convergent, therefore $\left\{c_{0}\right\}$ forms a weakly sequentially compact set in $L_{1}\left((N,R),x\,{\mathrm{d}}x\right)$. Hence by the de la Vallee Poussin theorem (Theorem 2.3) there exists a nonnegative, convex function $\Phi\in C^{\infty}([0,\infty))$, with $\Phi(0)=0$ and $\Phi^{\prime}(0)=1$ such that $\Phi^{\prime}$ is concave and satisfies

[TABLE]

A standard inequality (2.1), for $C^{1}$ convex functions gives us

[TABLE]

Multiplying this by $x_{i}\Delta x_{i}$ and summing over all $i$ gives

[TABLE]

Utilising equation (3.3) we can rewrite this as

[TABLE]

Recalling the definition of $S_{i}^{n}$ from (3.2), we see that it must be nonnegative. Additionally, Lemma 2.4(ii) and Lemma 4.1 give $\Phi^{\prime}\left(u_{C}^{n+1,i}\right)\geq 0$, hence we can drop the term involving $S_{i}^{n}$ from (5.1) and the inequality will still remain valid, giving us

[TABLE]

With some easy modification, equation (4.1) becomes the inequality

[TABLE]

which, if placed in the previous inequality, results in

[TABLE]

Utilising Lemma 2.4(i) with $x=u_{C}^{n,j}$ and $y=u_{C}^{n+1,i}$ and noting that the constants $\alpha(R)$ and $\beta(R)$ bound the average values $A_{j}$ and $B_{i,j}$, we get

[TABLE]

As $j$ is restricted to be greater than $i$ we have $x_{j}>x_{i}$ for admissible $j$ and $i$. This allows us to switch $x_{i}$ for $x_{j}$ in the second term and take this within the inner summation. Following this we expand the summation over $j$ to give

[TABLE]

If we change the index variable from $j$ to $i$ in the second summation and re-arrange then we obtain

[TABLE]

Some further manipulations produce

[TABLE]

By the final assumption of Remark 5.2, we have $1-K(R)\Delta t>0$ allowing us to divide through to get

[TABLE]

Repeated application of this inequality yields

[TABLE]

For values of $t$ in the interval $\tau_{n}=[t_{n},t_{n+1})$ this gives us

[TABLE]

An application of Jensen’s inequality [33, Theorem 2.2] allows us to switch the order of $\Phi$ and integration to get

[TABLE]

By assumption $6$ of Remark 5.2, that $K(R)\Delta t\leq\theta<1$, we deduce that

[TABLE]

which holds for all $t\in[0,T)$. Integrating the inequality with respect to $t$ from [math] to $T$ confirms the equiintegrability of the family $\left\{u_{C}^{h}\right\}$ in the space $L_{1}((N,R)\times[0,T),x\,{\mathrm{d}}x\,{\mathrm{d}}t)$. ∎

By Theorem 2.5 (Dunford–Pettis theorem), the sequence $\left\{u_{C}^{h}\right\}$ forms a weakly sequentially compact set in the space $L_{1}\left((N,R)\times[0,T),x\,{\mathrm{d}}x\,{\mathrm{d}}t\right)$. This implies the existence of a subsequence $\left\{u_{C}^{hj}\right\}$ and a function $u_{C}^{R}\in L_{1}\left((N,R)\times[0,T),x\,{\mathrm{d}}x\,{\mathrm{d}}t\right)$ such that $u_{C}^{hj}\rightharpoonup u_{C}^{R}$ in $L_{1}\left((N,R)\times[0,T),x\,{\mathrm{d}}x\,{\mathrm{d}}t\right)$ as $j\rightarrow\infty$ and $h_{j}\rightarrow 0$.

Remark 5.4**.**

From now on this convergent subsequence will be considered implicitly, unless otherwise stated; as such we now use the notation $\left\{u_{C}^{h}\right\}$ to denote such a convergent subsequence, the choice of which, we note, may not be unique.

5.2. Continuous Fragmentation Regime: Weak Solution

Having shown that our sequence of approximations converges (weakly) to a limit, we now aim to show that this limit provides a solution to our truncated equation (1.6). Precisely, we intend to show that the function $u_{C}^{R}$ satisfies the following criterion.

Definition 5.5**.**

The function $u_{C}^{R}$ is a weak solution of equation (1.6), if it satisfies

[TABLE]

for all functions $\varphi$, which are twice continuously differentiable functions on $(N,R)\times(0,T)$, such that $\varphi$ and each of its derivatives up to order $2$ may be continuously extended to $[N,R]\times[0,T)$, and such that for each fixed $x\in[N,R]$, the support of $\varphi$ with respect to $t$ is a compact subset of $[0,T)$. We denote the set of such extended functions by $C_{c}^{2}\left([N,R]\times[0,T)\right)$. Finally, we note that the weak formulation (5.5) was obtained from (1.6) in the usual manner, recalling the zero boundary conditions established in Lemma 1.1.

Remark 5.6**.**

We now make note of a property of the function $\varphi$ and its derivatives, which we will make use of in our analysis. As $\varphi$ has compact support and is identically zero outwith this support, its derivatives, both first and second, must also be zero outwith the support. Now within this compact support, $\varphi$ and its derivatives are continuous and so must be bounded functions.

Definition 5.7**.**

In the analysis which follows we make use of the following three approximations to $x$ over the domain $(N,R)$. First we have the left endpoint approximation, defined by

[TABLE]

Secondly we consider the midpoint approximation, defined by

[TABLE]

and finally we introduce the right endpoint approximation given by

[TABLE]

Remark 5.8**.**

It is a simple exercise to show that the three approximations, introduced above, converge point-wise (uniformly) to $x$ over the domain $(N,R)$ as the mesh parameter $h$ goes to [math]. The reader may find details given in [27, Lemma 5.2.3].

We are now in a position to proceed with our proof that $u_{C}^{R}$ is a weak solution to (1.6).

Definition 5.9**.**

Let $\varphi\in C_{c}^{2}\left([N,R]\times[0,T)\right)$ , then for sufficiently small $\Delta t$, the support of $\varphi$ with respect to $t$ lies within $[0,t_{M-1}]$. We define $\varphi_{i}^{n}$ as an approximation of $\varphi$ on $\Lambda_{i}\times\tau_{n}$ by

[TABLE]

with $\varphi_{i}^{M-1}=\varphi_{i}^{M}=0$ for admissible $i$ and define $\varphi_{I_{h}}^{n}=0$ for all $n$.

Rearranging equation (3.3), multiplying by $\varphi_{i}^{n}$ and summing over $n=0,\ldots,M-1$ and $i=0,\ldots,I_{h}-1$, gives us

[TABLE]

Rearrangement of the summations and utilising the compact support of $\varphi$ and the zero boundary flux gives us the following equality:

[TABLE]

The above equality can be seen as the discrete equivalent of the weak formulation (5.5). Our approach now involves taking the limit as $h\rightarrow 0$ of (5.2) and showing that we do indeed obtain (5.5) with $u_{C}^{R}$ as a weak solution. Observing the terms of (5.5) we see that the integrals are with respect to the measure ${\mathrm{d}}x\,{\mathrm{d}}t$ whilst we have shown that convergence occurs in the space with weighted measure $x\,{\mathrm{d}}x\,{\mathrm{d}}t$. At this point we highlight the use of Theorem 2.1 to switch spaces but retain convergence.

Theorem 5.10**.**

Under the assumptions outlined in Remark 5.2, the function $u_{C}^{R}$ obtained as the limit of the sequence $\left\{u_{C}^{h}\right\}$, is a weak solution of our equation, satisfying (5.5).

Proof.

Looking initially at the first two terms of (5.2), we can express them as follows:

[TABLE]

Considering the first of the double integrals, let $(x,t)\in(N,R)\times(0,T)$, then in the case that $0<\Delta t\leq t$, Taylor expansions of the $\varphi$ terms about the point $(x,t)$ give us

[TABLE]

[TABLE]

Simple cancellations and recalling the condition (3.1) relating $h$ and $\Delta t$ give us

[TABLE]

The expression on the left-hand side of (5.5) can thus be seen to converge pointwise to $\varphi_{t}$ on $(N,R)\times(0,T)$ as the mesh size goes to 0. Furthermore, since the derivatives of $\varphi$ are bounded as per Remark 5.6, we can bound the left-hand side of (5.5) on $(N,R)\times(0,T)$, with a bound that is uniform w.r.t $h$, as $h\searrow 0$. As noted in Remark 5.8, the functions $X^{h}(x)$ converge pointwise to $x$ on $(N,R)$ as $h\searrow 0$, and are clearly bounded by $R$ for all values of $x$ and $h$. Noting that $(N,R)\times[0,T)\backslash(N,R)\times(0,T)$ is of measure [math] with respect to the measure $\,{\mathrm{d}}x\,{\mathrm{d}}t$, we see that the terms accompanying $u_{C}^{h}$ in our double integral satisfy the conditions for $\left\{g_{h}\right\}$ from Theorem 2.2. We have shown previously that $u_{C}^{h}\rightharpoonup u_{C}^{R}$ in $L_{1}\left((N,R)\times[0,T),\,x\,{\mathrm{d}}x\,{\mathrm{d}}t\right)$ and by Theorem 2.1 this is also the case in $L_{1}\left((N,R)\times[0,T),\,{\mathrm{d}}x\,{\mathrm{d}}t\right)$ and so an application of Theorem 2.2 gives us

[TABLE]

as the mesh size parameter $h$ goes to 0. Next, we consider the second term appearing above. Since $\varphi$ is $C_{c}^{2}\left([N,R]\times[0,T)\right)$, its derivatives are bounded, allowing us to deduce that

[TABLE]

for all $x\in(N,R)$, as we now demonstrate. Consider the following:

[TABLE]

Expressing $\varphi(\xi^{h}(x),t)$ using a Taylor expansion about $(x,t)\in(N,R)\times(0,T)$, and recalling Remark 5.6 about the derivatives of $\varphi$ we get

[TABLE]

Hence bounding the derivative $\partial\varphi/\partial x$ and noting that $\left|\xi^{h}(x)-x\right|\leq h$ gives us

[TABLE]

for some constant $C_{1}$ independent of $h$ and $\Delta t$. Similarly, expanding $\varphi(x,t)$ about $(x,0)$, where $x\in(N,R)$, produces

[TABLE]

for $t\geq 0$. The use of the notation $\tfrac{\partial\varphi}{\partial t}(x,0_{+})$ signifies we are considering the right derivative of $\varphi$ with respect to $t$ at $t=0$. The expansion (5.9) then leads to

[TABLE]

for some other constant $C_{2}$, independent of $h$ and $\Delta t$. Returning to (5.2) we have

[TABLE]

Hence (5.7) does indeed hold for $x\in(N,R)$, furthermore the convergence is uniform with respect to $x$. Together with the pointwise convergence of $X^{h}(x)$ to $x$, and from Remark 3.2, the $L_{1}(N,R)$ strong convergence of $u^{h}_{C}(x,0)$ to the restriction of $c_{0}$ to $(N,R)$, another application of Theorem 2.2 yields

[TABLE]

Moving on to the third term of equation (5.2), for $t\in\tau_{n}$ and $x\in\Lambda_{i}$ we can write the numerical flux as an integral as follows:

[TABLE]

Then the third term of equation (5.2) is given by

[TABLE]

Expressed in full this gives us the following, after a switch in the order of integration:

[TABLE]

Due to the boundedness of the partial derivative $\varphi_{x}$ and the $L_{\infty,loc}$ property of the functions $a$ and $b$, for almost all fixed $(x,w)\in(N,R)\times(N,R)$, the product $\chi_{[\Xi^{h}(x),R]}(v)a^{h}(v)b^{h}(w|v)\frac{\partial\varphi}{\partial x}(x,t)$ is a bounded (uniformly w.r.t. $h$) function of $v$ and $t$. Also, as a consequence of Remarks 3.1 and 5.8, it converges pointwise almost everywhere (w.r.t the measure $\,{\mathrm{d}}v\,{\mathrm{d}}t$) on $(v,t)\in(N,R)\times[0,T)$ to $\chi_{[x,R]}(v)a(v)b(w|v)\frac{\partial\varphi}{\partial x}(x,t)$ as $h\rightarrow 0$. Since $u_{C}^{h}\rightharpoonup u_{C}^{R}$ in $L_{1}((N,R)\times[0,T),\,{\mathrm{d}}x\,{\mathrm{d}}t)$, an application of Theorem 2.2 gives us

[TABLE]

Using the local boundedness of $a$ and $b$, along with the boundedness of the partial derivative $\frac{\partial\varphi}{\partial x}$, and the boundedness of the sequence $\left\{u_{C}^{h}\right\}$ in $L_{1}((N,R)\times[0,T),\,{\mathrm{d}}x\,{\mathrm{d}}t)$, the left–hand side above can be bounded by a constant. It is easily seen that $\chi_{[N,\Xi^{h}(x)]}(w)X^{h}(w)$ converges pointwise to $\chi_{[N,x]}(w)w$ as $h\rightarrow 0$ and can be bounded by $R$ over our domain of interest. Therefore applying the Lebesgue dominated convergence theorem [33, Theorem 1.8], we get that (5.2) converges to

[TABLE]

Therefore, in the limit as $h\rightarrow 0$, the third term of (5.2) coincides with the third term of (5.5).

Now the fourth term from equation (5.2) is given fully by

[TABLE]

The pointwise convergence of $a^{h}$, $b_{j}^{h}$ and $\xi^{h}$ along with the continuity of $\varphi$ means that

[TABLE]

for all $t\in[0,T)$ and almost all $v\in(R,N)$ as $h\rightarrow 0$. Since $a$ and $b_{i}$ are in $L_{\infty,loc}([N,\infty))$ and $\varphi$ is $C^{2}$ on $[N,R]\times[0,T)$ with compact support (hence is a bounded function), the expressions on either side belong to $L_{\infty}((N,R)\times[0,T))$, with the left-hand side being uniformly bounded w.r.t. $h$. Hence, with $u_{c}^{h}\rightharpoonup u_{C}^{R}$ in $L_{1}((N,R)\times[0,T),\,{\mathrm{d}}v\,{\mathrm{d}}t)$, applying Theorem 2.2, yields

[TABLE]

Taken together (5.2), (5.10), (5.2) and (5.2) show that $u_{C}^{R}$ satisfies (5.5) for all $\varphi\in C_{c}^{2}([N,R]\times[0,T))$, and hence $u_{C}^{R}$ is a weak solution, as set out in Definition 5.5. ∎

5.3. Discrete Fragmentation Regime: Convergence

Now let us consider the discrete regime approximations. This is treated by a similar approach to the one we adopted for the continuous regime equation, but as a first step we establish a bound on the values $u_{D}^{n,i}$.

Lemma 5.11**.**

There exists a constant $C$, independent of $h$ and $R$, such that for all values of $n$ and $i$ we have

[TABLE]

Proof.

The nonnegativity of $u_{D}^{n,i}$ follows from Lemma 4.1. We shall therefore concentrate on the upper bound. From Lemma 4.3 we have, for all admissible $n$, that the following holds:

[TABLE]

Therefore, for all $n$ and $i$ we have that

[TABLE]

where the constant $C$ is independent of the mesh parameter $h$ and the truncation parameter $R$. ∎

Theorem 5.12**.**

For each $i=1,\ldots,N$, the family $\left\{u^{h}_{Di}\right\}$ forms a sequentially weakly compact set in $L_{1}\left([0,T)\right)$, hence must have a weakly convergent subsequence.

Proof.

The bound obtained in Lemma 5.11 allows us to easily establish equiboundedness and equiintegrability in $L_{1}\left([0,T)\right)$ for each of the families $\left\{u^{h}_{Di}\right\}$ as follows:

[TABLE]

Now let $\Phi$ be any function of the nature described in Theorem 2.3. Since $\Phi$ is increasing, the established bound for $u^{n,i}_{D}$ allows us to deduce that

[TABLE]

Hence each of the families $\left\{u^{h}_{Di}\right\}$ is equiintegrable. By the Dunford–Pettis theorem (Theorem 2.5), each of the families form a weakly sequentially compact set in $L_{1}\left([0,T)\right)$. As such, they all contain some weakly convergent subsequence. ∎

Remark 5.13**.**

We note that what we seek is a collection of values $\left\{h^{j}\right\}$, such that all of the sequences $\left\{u^{h^{j}}_{Di}\right\}$, for $i=1,\ldots,N$, converge weakly, as $j\rightarrow\infty$ and $h_{j}\rightarrow 0$. We achieve this by means of a diagonal argument, which we now outline. Knowing that the family $\left\{u^{h}_{D1}\right\}$ has a weakly convergent subsequence, let us denote the corresponding sequence of $h$-values by $\left\{h^{j}\right\}_{j=1}^{\infty}$ and consider the family $\left\{u^{h^{j}}_{D2}\right\}_{j=1}^{\infty}$. As this set satisfies the equiboundedness and equiintegrability conditions of the Dunford–Pettis theorem, it too must have a weakly convergent subsequence. Extracting this subsequence and denoting the corresponding $h$-values by $\left\{h^{j_{n}}\right\}_{n=1}^{\infty}$, we then have both $\left\{u^{h^{j_{n}}}_{D1}\right\}$ and $\left\{u^{h^{j_{n}}}_{D2}\right\}$ converging (weakly) as $n\rightarrow\infty$ and $j_{n}\rightarrow\infty$ . We can continue this process, working through each of the families $\left\{u^{h}_{Di}\right\}$, until we have a set of common $h$-values, $\left\{h_{j}^{\prime}\right\}_{j=1}^{\infty}$, for which all the subsequences $\left\{u^{h_{j}^{\prime}}_{Di}\right\}_{j=1}^{\infty}$ are (weakly) convergent as $j\rightarrow\infty$ and $h_{j}^{\prime}\rightarrow 0$.

From now on these convergent subsequences are considered implicitly and we use $\left\{u^{h}_{Di}\right\}$ to denote said subsequences, unless otherwise stated. Let us denote the weak limit of $\left\{u^{h}_{Di}\right\}$ by $u_{Di}^{R}$ (note the upper case superscript notation for the limit).

5.4. Discrete Fragmentation Regime: Weak Solution

Having established the convergence of our sequence of approximations we now aim to determine whether the limit produced provides a solution to the equation (1.8) and if so in what sense. As such, following on from Definition 5.5, we introduce

Definition 5.14**.**

We say that the function $u_{Di}^{R}$ is a weak solution of equation (1.8) if it satisfies

[TABLE]

for any $\phi\in C_{c}^{2}([0,T))$, where $C_{c}^{2}([0,T))$ is defined in an analogous fashion to $C_{c}^{2}\left([N,R]\times[0,T)\right)$ from Definition 5.5.

Theorem 5.15**.**

The functions $u_{Di}^{R}$ obtained as weak limits of the sequences $\left\{u^{h}_{Di}\right\}$ are indeed weak solutions of (1.8), satisfying equation (5.14) for any $\phi\in C_{c}^{2}([0,T))$.

Proof.

For such a function $\phi$, let us denote its approximation over $\tau_{n}$ by $\phi^{n}$, which is defined as

[TABLE]

and $\phi^{M}=0$. Multiplying (3.6) by $\phi^{n}$ and summing over $n$ from [math] to $M-1$, gives us the following equality:

[TABLE]

Since $\phi$ is compactly supported, for sufficiently small $\Delta t$ we have $\phi^{M-1}=0$; then, further manipulation of the first term yields

[TABLE]

Looking more closely at the first term above we can rewrite it as

[TABLE]

Assuming that $t\in(0,T)$ and $\Delta t\leq t$, then a Taylor series expansion of $\phi(t-\Delta t)$ about $t$ gives

[TABLE]

Therefore, we have

[TABLE]

As such, the left-hand side of (5.17) can be seen to converge pointwise to $\phi_{t}$ on $(0,T)$, as $h$, and by condition (3.1), $\Delta t$ goes to [math]. By an analogous argument to that used for $\varphi$ and its derivatives, $\phi$ and its derivative $\phi_{t}$ must be bounded, therefore we can bound the left-hand side above, with the bound being uniform w.r.t $h$. Then, as ${u_{D}^{h}}_{i}\rightharpoonup{u_{D}^{R}}_{i}$ in $L_{1}\left([0,T)\right)$, applying Theorem 2.2, as before, gives us

[TABLE]

By definition, $u_{Di}^{n}={d_{0}}_{i}$, and since $\phi$ is $C^{2}$ with compact support, its derivative must be bounded, from which we deduce that

[TABLE]

as $h$ goes to 0. Therefore

[TABLE]

By defining $b_{i,i}$ to be $-1$, we can combine the third and fourth terms of (5.4), writing them as

[TABLE]

and since $u_{Dj}^{h}\rightharpoonup u_{Dj}^{R}$ in $L_{1}\left([0,T)\right)$, for each $j$, we have

[TABLE]

giving us the third and fourth terms of our weak formulation (5.14). Rewriting the final term of our discrete relation, we get

[TABLE]

From Remark 3.1 we have $a^{h}(y)$ and $b_{i}^{h}(y)$ converging pointwise to $a(y)$ and $b_{i}(y)$ respectively, and along with $\phi$ are bounded (uniformly with respect to $h$), a final application of Theorem 2.2 allows us to deduce that

[TABLE]

as the mesh size parameter $h\rightarrow 0$. Taking the results (5.18), (5.19), (5.4) and (5.21), we see that by letting $h\rightarrow 0$ in (5.4) we obtain the weak formulation (5.14), hence $u_{Di}^{R}$ is indeed a weak solution of (1.8). ∎

In this section we established the weak convergence of a subsequence of our sequence of approximate solutions as the mesh parameter was decreased to zero. The limits were shown to provide a set of weak solutions to the truncated equations (1.6) and (1.8). However, there are a number of questions which remain unanswered which we seek to address in the following section.

6. Uniqueness and Differentiability of Solutions

In the previous section, we formed approximate solutions to a truncated version of our system. A subsequence of these approximations was shown to converge to a weak solution to our problem, as the underlying mesh was refined. This convergence of subsequences, rather than the full sequence, raises the possibility of nonunique solutions, with each convergent subsequence possibly offering a different solution. In this section we seek to address this, showing that any limits must coincide, providing a unique solution. Further, we would like to establish whether this solution may in fact display extra regularity, as we might expect from the results in [15].

6.1. Continuous Regime

Returning to equation (1.5), we introduce the space $X_{C}^{R}=L_{1}\left((N,R),x\,{\mathrm{d}}x\right)$ with the aim of recasting the equation as an abstract Cauchy problem, as was carried out for (1.1) in [15]. Motivated by the terms appearing on the right-hand side, of (1.5) we define the following linear operators $A_{C}^{R}$ and $B_{C}^{R}$ on the space $X_{C}^{R}$:

[TABLE]

with the respective domains

[TABLE]

Assuming that the functions $a$ and $b$ retain the properties imposed in Remark 5.2, in particular $a\in L_{\infty,loc}\left([N,\infty)\right)$ and $b\in L_{\infty,loc}\left([N,\infty)\times[N,\infty)\right)$, then, the following property holds for the operators $A_{C}^{R}$ and $B_{C}^{R}$.

Lemma 6.1**.**

The operators $A_{C}^{R}$ and $B_{C}^{R}$ are bounded linear operators on the space $X_{C}^{R}$, with $\left\|B_{C}^{R}f\right\|_{X_{C}^{R}}\leq\left\|A_{C}^{R}f\right\|_{X_{C}^{R}}$ for all $f\in X_{C}^{R}$.

Proof.

We will first consider the operator $A_{C}^{R}$. Let $f\in X_{C}^{R}$; then we have

[TABLE]

Therefore $A_{C}^{R}$ is a bounded operator on the space $X_{C}^{R}$. The boundedness of $B_{C}^{R}$ in $X_{C}^{R}$ can be seen as follows. Let $f\in X_{C}^{R}$; then

[TABLE]

The change in the order of integration can be justified by the nonnegativity of the integrand along with Tonelli’s theorem. The inequality in going from the third to the fourth line comes as a result of the mass conservation condition (1.3). ∎

Equation (1.5) is then recast as the following abstract Cauchy problem in the space $X_{C}^{R}$:

[TABLE]

Here $u_{C}^{R}$ denotes an $X_{C}^{R}$-valued function rather than the scalar-valued function of two variables from the previous section. However due to the relationship between the spaces $L_{1}\left(I,L_{1}(\Omega,\,{\mathrm{d}}\mu)\right)$ and $L_{1}\left(\Omega\times I,\,{\mathrm{d}}\mu\,{\mathrm{d}}t\right)$, we may switch between the two, with each $L_{1}$-valued solution to (6.1) providing us with a scalar-valued solution to (1.5) and vice versa.

Lemma 6.2**.**

The operator $\left(A_{C}^{R}+B_{C}^{R}\right)$ generates a uniformly continuous semigroup of positive contractions on $X_{C}^{R}$.

Proof.

As a bounded linear operator on the space $X_{C}^{R}$, the sum $A_{C}^{R}+B_{C}^{R}$ generates a uniformly continuous semigroup $\left(T_{R}(t)\right)_{t\geq 0}$ on $X_{C}^{R}$, [34, Chapter 1, Theorem 3.7 and Chapter 2, Corollary 1.5]. Furthermore, the Kato–Voigt perturbation theorem, [30, Corollary 5.17] as applied in [15, Theorem 3.2], is readily utilised in the case $\left(A_{C}^{R}+B_{C}^{R},X_{C}^{R}\right)$ to give us an ‘extension’ of $\left(A_{C}^{R}+B_{C}^{R},X_{C}^{R}\right)$ as a generator of a substochastic semigroup. Now as $A_{C}^{R}+B_{C}^{R}$ is defined and bounded on all of $X_{C}^{R}$, this extension must be $A_{C}^{R}+B_{C}^{R}$ itself and by [35, Theorem 2.6], the substochastic semigroup generated must be $\left(T_{R}(t)\right)_{t\geq 0}$. Hence $\left(T_{R}(t)\right)_{t\geq 0}$ is a uniformly continuous semigroup of positive contractions. ∎

From standard results concerning strongly continuous semigroups, [36, Theorem 2.40 and Theorem 2.41], the existence of the semigroup $\left(T_{R}(t)\right)_{t\geq 0}$ on $X_{C}^{R}$, provides a unique strong solution to equation (6.1), given by $u_{C}^{R}(t)=T_{R}(t)c_{0}^{R}$. Additionally, by [34, Chapter 2, Proposition 6.4] this is also a unique mild solution, satisfying an equation of the form

[TABLE]

We are able to take the operator $A_{C}^{R}+B_{C}^{R}$ inside the integral as a consequence of its boundedness by applying [37, Proposition 1.1.7.]. Hence our equation (6.2) corresponds with the mild solution form of [38, Definition 1.10]. The article [38] also provides a notion of an $X_{C}^{R}$-valued weak solution of equations of type (6.1), which we outline here for our specific example.

Definition 6.3**.**

The function $u_{C}^{R}:[0,T)\longrightarrow X_{C}^{R}$ is a weak solution of (6.1), if for all $\phi$ in $L_{\infty}\left((N,R)\right)$ (the dual space of $X_{C}^{R}$), we have that $t\rightarrow\langle u_{C}^{R}(t),\phi\rangle$ is locally integrable in $(0,T)$ and

[TABLE]

for all $\psi\in C^{\infty}\left([0,T)\right)$ with compact support, where $\langle g,\phi\rangle$ denotes the duality pairing of $g$ and $\phi$.

Remark 6.4**.**

If $D\subseteq{X_{C}^{R}}^{\prime}$ is dense in the weak-$\ast$ topology, then it is sufficient to show that (6.3) holds for all $\phi\in D$ to establish $u_{C}^{R}:[0,T)\rightarrow X_{C}^{R}$ as a weak solution of (6.1); see [38, Definition 1.9].

The results so far have provided us with the existence of unique strong and mild solutions to (6.1). We now show that, for our case, any mild solution satisfying (6.2) must necessarily be a weak solution as in Definition 6.3 and vice-versa, providing the existence and uniqueness of a weak solution.

Theorem 6.5**.**

The function $u_{C}^{R}:[0,\infty)\longrightarrow X_{C}^{R}$ provided by the semigroup $\left(T_{R}(t)\right)_{t\geq 0}$ is the unique weak solution to equation (6.1), satisfying Definition 6.3 over any time interval $[0,T)$ where $T<\infty$.

Proof.

From the analysis above we have the semigroup $\left(T_{R}(t)\right)_{t\geq 0}$ providing a unique mild (strong) solution $u_{C}^{R}:[0,\infty)\longrightarrow X_{C}^{R}$ to equation (6.1). The result [38, Theorem 1.20] tells us that if the right-hand side of (6.1) is integrable then, given either a weak or mild solution to our evolution equation (6.1), the solution can be modified on a set of measure zero to obtain a solution of the other form. Therefore in our case, if the conditions of [38, Theorem 1.20] hold, the existence of a unique mild solution will provide a unique weak solution (as given in Definition 6.3) on each finite time interval $[0,T)$, with this solution being given by the semigroup $\left(T_{R}(t)\right)_{t\geq 0}$.

Let $0<T<\infty$; and let us consider the following integral:

[TABLE]

In going to the final line, the bounding of the inner integral of the second term from the previous line comes from the mass conservation condition (1.3), whilst the subsequent norm inequality relies on the fact that the semigroup $\left(T_{R}(t)\right)_{t\geq 0}$ consists of contractions.

By [39, Chapter 2, Theorem 2 and Theorem 4], the above bound implies that the right-hand side of equation (6.1) is integrable over the interval $[0,T)$ and so by [38, Theorem 1.20], any mild or weak solutions must agree (up to sets of measure zero). Therefore equation (6.1) has a unique weak solution given by $u_{C}^{R}(t)=T_{R}(t)c_{0}^{R}$, which is in fact a strong solution. ∎

Let $u_{C}^{R}(x,t)$ be a scalar representation of the semigroup solution $(u_{C}^{R}(t))(x)=\left(T_{R}(t)c_{0}^{R}\right)(x)$. Then applying Definition 6.3 to our example and noting the equivalence of (1.5) and (1.6) as detailed fully in [27, Appendix C], we get that equation (6.3) is equivalent to

[TABLE]

for all $\varphi$ of the form $\varphi(x,t)=\phi(x)\psi(t)$ where $\phi\in C^{\infty}_{c}\left((N,R)\right)$ and $\psi\in C^{\infty}_{c}\left([0,T)\right)$, due to the weak-$\ast$ density of $C^{\infty}_{c}\left((N,R)\right)$ in $L_{\infty}\left((N,R)\right)$, via Remark 6.4.

Having determined a one-to-one correspondence between scalar-valued weak solutions satisfying (6.1) and $X_{C}^{R}$-valued strong solutions of the abstract Cauchy problem (6.1), we are now in a position to establish the uniqueness and differentiability of the weak solutions of the previous section.

Theorem 6.6**.**

The weak solution obtained as the limit of the sequence of approximate solutions for the continuous regime in Theorems 5.3 and  5.10 is unique, continuously differentiable with respect to $t$ on any interval $[0,T)$ and satisfies equation (1.5) directly, except perhaps on a set of measure zero.

Proof.

It is easily seen that any scalar-valued function $u_{C}^{R}(x,t)$ satisfying Definition 5.5 will immediately satisfy the equation (6.1) above. From any such scalar-valued function we may define a function $u_{C}^{R}:[0,T)\longrightarrow X_{C}^{R}$ via $\left(u_{C}^{R}(t)\right)(x):=u_{C}^{R}(x,t)$, for almost all $(x,t)\in(N,R)\times[0,T)$. Since the scalar function satisfies equation (6.1), the $X_{C}^{R}$-valued function must provide a weak solution, as defined in Definition 6.3, to the abstract Cauchy problem (6.1). By Theorem 6.5, this weak solution must necessarily also be the unique (up to sets of measure zero) strong solution of (6.1). Hence for each scalar weak solution $u_{C}^{R}(x,t)$ satisfying Definition 5.5 there is a unique corresponding $X_{C}^{R}$-valued strong solution $u_{C}^{R}$ to the associated abstract Cauchy problem, whereby the original function $u_{C}^{R}(x,t)$ is one scalar representation of the strong solution. Since the strong solution is unique, as are its scalar representations (up to measure zero), it follows that the weak solution must be unique, up to a set of measure zero. Further, as a representation of a strongly differentiable $X_{C}^{R}$-valued function, by [30, Theorem 2.40], $u_{C}^{R}(x,t)$ is continuously differentiable with respect to $t$, except perhaps on a set of zero measure. Further, by a similar argument as applied at the end of [30, Theorem 8.3] or directly as in [27, Theorem 3.2.7.], $u_{C}^{R}(x,t)$ can be seen to directly satisfy equation (1.5) almost everywhere. ∎

This result greatly strengthens those of the previous section, where before we had only the existence of a (weakly) convergent subsequence and the possibility of the numerical scheme converging to multiple weak solutions. We now know that the limit solution must necessarily be unique, continuously differentiable with respect to $t$ and a solution, in the classical sense, of the truncated fragmentation equation.

6.2. Discrete Regime

Recalling the truncated discrete regime equation (1.8), we have for $i=1,\ldots,N$:

[TABLE]

With the aim of recasting these equations as an inhomogeneous abstract Cauchy problem as in [15], we introduce the space $X_{D}=\mathbb{R}^{N}$, equipped with the weighted norm:

[TABLE]

The equations (6.5) then become

[TABLE]

where $A_{D}$ and $B_{D}$ are defined on $X_{D}$, by

[TABLE]

$u_{C}^{R}$ is the truncated semigroup solution of (6.1) and where $C_{R}:D(C_{R})\subseteq X_{C}^{R}\rightarrow X_{D}$ is given by

[TABLE]

for $i=1,\ldots,N$. The fact that $D(C_{R})=X_{C}^{R}$ is a consequence of the $L_{\infty,loc}$ boundedness of $a$ and each $b_{i}$. Recall that since the space $X_{D}$ is finite-dimensional, the operators $A_{D}$ and $B_{D}$ are bounded and therefore by [34, Chapter 1, Theorem 3.7 and Chapter 2, Corollary 1.5] their sum must generate a uniformly continuous semigroup $\left(T(t)\right)_{t\geq 0}$ on $X_{D}$. We then consider the term $C_{R}[u_{C}^{R}(t)]$ as a perturbation and (6.6) as an inhomogeneous abstract Cauchy problem. The following lemma establishes the differentiability of this perturbation term, a property that we will require in showing the existence of solutions of equation (6.6), along with determining their nature.

Lemma 6.7**.**

The term $C_{R}[u_{C}^{R}(t)]$ from (6.6) is strongly differentiable (in the space $X_{D}$) with respect to $t$, at almost all points of $[0,T]$. Furthermore its derivative is given by $C_{R}[\frac{{\mathrm{d}}}{{\mathrm{d}}t}u_{C}^{R}(t)]$.

Proof.

Recalling the $L_{\infty,loc}$ boundedness of $a$ and the bound $b_{i}(y)\leq y$ for each $i\in\left\{1,2,\ldots,N\right\}$, which is easily derived from (1.3), we have that

[TABLE]

where $\alpha(R)$ is the essential supremum of $a$ over $[N,R]$. From the differentiability of $u_{C}^{R}$, by letting $h\rightarrow 0$ on both sides of the above calculation, we may deduce that $C_{R}[u_{C}^{R}(t)]$ is differentiable (in the space $X_{D}$) at almost all points of $[0,T]$, with derivative $C_{R}[\frac{{\mathrm{d}}}{{\mathrm{d}}t}u_{C}^{R}(t)]$. ∎

Having shown the differentiability of the perturbation term $C_{R}[u_{C}^{R}(t)]$, we now look at its derivative more closely, showing that it is integrable, belonging to the space $L_{1}((0,T),X_{D})$ and in doing so establish the existence of a unique strong solution to equation (6.6).

Theorem 6.8**.**

The derivative of $C_{R}[u_{C}^{R}(t)]$ belongs to the space $L_{1}((0,T),X_{D})$. As such equation (6.6) has a unique strong solution, which is given by

[TABLE]

Proof.

If we take the $X_{D}$-norm of the derivative $\frac{{\mathrm{d}}}{{\mathrm{d}}t}C_{R}[u_{C}^{R}(t)]$ established in the previous lemma and integrate from [math] to $T$, then we obtain

[TABLE]

Recalling the abstract Cauchy problem (6.1), since the operators $A_{C}^{R}$ and $B_{C}^{R}$ are bounded and since $u_{C}^{R}$ is given by a contraction semigroup, we have

[TABLE]

where $M$ is a constant such that $\left\|A_{C}^{R}+B_{C}^{R}\right\|\leq M$. Inserting this into (6.2) gives us

[TABLE]

Therefore the derivative of $C_{R}[u_{C}^{R}(t)]$ belongs to the space $L_{1}((0,T),X_{D})$; hence, by [35, Chapter 4, Corollary 2.2 and Corollary 2.10], the equation (6.6) has a unique strong solution $u_{D}^{R}:[0,T)\rightarrow X_{D}$ given by (6.7). ∎

Having established the existence of a unique strong solution to equation (6.6) given by (6.7), this solution must also provide us with a unique mild solution to our equation. Now we consider the possibility of weak solutions, as defined in Definition 6.3. We aim to show that any weak solution of equation (6.6) must also be a mild solution (permitting changes on sets of measure zero), and hence the weak solution must be unique and differentiable.

Theorem 6.9**.**

Given an integrable weak solution $u_{D}^{R}:[0,T)\rightarrow X_{D}$ of equation (6.6) as defined in Definition 6.3, then it must also be a strong solution. Therefore any integrable weak solution must be unique up to sets of measure zero and differentiable in $X_{D}$.

Proof.

Considering the right-hand side of equation (6.6), taking the norm in $X_{D}$ and integrating from [math] to $T$ gives us

[TABLE]

The assumption that $u_{D}^{R}(t)$ is integrable and [39, Chapter 2, Theorem 2 and Theorem 4] allow us to deduce that the first of these integrals must be finite, whilst recalling that $u_{C}^{R}(t)$ was given by a contraction semigroup immediately enables us to bound the second integral above. Hence, [39, Chapter 2, Theorem 2 and Theorem 4], the right-hand side of equation (6.6) is integrable and therefore [38, Theorem 1.20] tells us that the weak solution $u_{D}^{R}(t)$ (allowing for changes on sets on measure zero) must also be a mild solution. Since, by Theorem 6.8, equation (6.6) has a unique mild solution which is in fact a strong solution, the weak solution $u_{D}^{R}(t)$ we started with must agree with the strong solution (up to sets of measure zero) and therefore is unique and differentiable. ∎

Having established that any integrable weak solution of equation (6.6), in the sense of Definition 6.3, is also a strong solution, we now set out to prove that the solutions of equation (6.5) obtained previously as the limit of our numerical scheme, provide us with such a weak solution and in the process establish their uniqueness and differentiability.

Lemma 6.10**.**

The weak solutions $u_{Di}^{R}(t)$ to the equations (6.5), obtained from our numerical scheme, when taken as the components of $u_{D}^{R}:[0,T)\rightarrow X_{D}$, produce a $u_{D}^{R}$ which is integrable.

Proof.

The weak solutions $u_{Di}^{R}$ to the equations (6.5), constructed in the previous section, were obtained as the weak limits in $L_{1}(0,T)$ of the sequences $\left\{u^{h}_{Di}\right\}$ as we let $h\searrow 0$. By the weak lower semicontinuity of the norm [33, Theorem 2.11], and using the bound (5.14), we obtain

[TABLE]

where $C$ denotes the constant from (5.14). If we take the components of the function $u_{D}^{R}:[0,T)\rightarrow X_{D}$ to be given by $u_{Di}^{R}$ for $i=1,\ldots,N$ then we get

[TABLE]

Therefore the function $u_{D}^{R}:[0,T)\rightarrow X_{D}$ formed by taking $u_{Di}^{R}$ as its $i^{th}$ component is integrable. ∎

Theorem 6.11**.**

The weak solutions $u_{Di}^{R}:[0,T)\rightarrow\mathbb{R}$ for $i=1,\ldots,N$ of the equations (5.14) obtained in the previous section agree with the components of the strong solution established in Theorem 6.8 and hence are unique (up to sets of measure zero) and differentiable.

Proof.

Let us consider our abstract equation (6.6) with the aim of rewriting it in a weak formulation as in (6.3). First let us note that the dual space of $X_{D}$ is $\mathbb{R}^{N}$ and the duality pairing $\langle\cdot,\cdot\rangle$ appearing in (6.3) is given by the standard inner product on $\mathbb{R}^{N}$. Let $u_{D}^{R}:[0,T)\rightarrow X_{D}$ be a weak solution of equation (6.6) with the components $u_{Di}^{R}:[0,T)\rightarrow\mathbb{R}$ for $i=1,\ldots,N$. Then, in this case, the function $F:(0,T)\rightarrow X_{D}$ appearing in (6.3) is given componentwise by

[TABLE]

for $i=1,\ldots,N$. Therefore the weak formulation of equation (6.6) becomes:

[TABLE]

for any $\phi=\left(\phi_{1},\ldots,\phi_{N}\right)\in\mathbb{R}^{N}$ and $\psi\in C_{c}^{\infty}\left([0,T)\right)$, where $u_{C}^{R}(\cdot,t)$ is the (unique) real-valued representation of the $X_{C}^{R}$-valued $u_{C}^{R}:[0,T)\rightarrow X_{C}^{R}$.

Comparing this with equation (5.14) of the previous section, it is easily seen that the $u_{Di}^{R}$ obtained there provide us with a solution to the above equation and so taking these $u_{Di}^{R}:[0,T)\rightarrow\mathbb{R},$ $i=1,\ldots,N$, as the components of an $X_{D}$-valued function $u_{D}^{R}:[0,T)\rightarrow X_{D}$ provides us with an integrable weak solution to equation (6.6). As any such $u_{D}^{R}:[0,T)\rightarrow X_{D}$ must be unique and differentiable in $X_{D}$, the components $u_{Di}^{R}:[0,T)\rightarrow\mathbb{R}$ must be unique (up to sets of measure zero) and differentiable in the traditional sense. ∎

In this article we have shown the convergence of our approximate solutions to a weak solution of the truncated problem given by equations (1.5) and (1.8), and the equivalence of this weak solution to the unique strong/classical solution of the truncated problem. However, it is possible using standard arguments along the lines of [30, Section 8.3.2] to show these truncated solutions converge, in the sense of the appropriate space, to the unique strong/classical solution of the untruncated problem as given by (1.1) and (1.2), whose existence was established in [15]. For further details of this convergence in the specific case of the mixed discrete-continuous model, the reader is directed to consult [27, Chapter 6]. Furthermore, for an experimental study of this convergence, the factors influencing it and therefore the selection of a suitable value for the truncation parameter $R$, the reader is directed to [27, Section 7.5].

7. Numerical Experiments

To assess the efficacy of our numerical scheme, we tested it on the power law model as set out in [15, Section 7], where the continuous equation was defined by

[TABLE]

The discrete equation was specified by the following values for $a_{i}$ and $b_{i,j}$

[TABLE]

and the continuous to discrete distribution functions $b_{i}(y)$ were given by

[TABLE]

It can be easily verified that conditions (1.3) and (1.4) are satisfied by these choices. The parameters $N$ and $R$ were set at $5$ and $15$ respectively, and an initial state assumed, with

[TABLE]

and $d_{0}$ being the $N-$vector consisting entirely of $1^{\prime}$s. The parameters $\alpha$ and $\nu$ were varied, taking all possible combinations of $\alpha\in\left\{0.5,0.1,-0.5,-1,-2\right\}$ and $\nu\in\left\{0,-0.5,-1,-1.5\right\}$, with the final time $T$ selected in each case to allow the system to reach a near equilibrium state.

The approximate solutions generated by the numerical scheme were compared to the exact solutions derived in [15, Section 7.1], with the discrepancy being measured by taking the relative error with respect to the norm on $L_{1}\left([0,T),X_{D}\right)\times L_{1}\left([0,T),X_{C}^{R}\right)$. That is, supposing $u^{h}=(\underline{u}_{D}^{h}(t),u_{C}^{h}(x,t))$ is our approximation of an exact solution $u=(\underline{u}_{D}(t),u_{C}^{R}(x,t))$, then we measure the error via

[TABLE]

where the norm $\|\cdot\|$, is given by

[TABLE]

For each model configuration, we computed approximate solutions over a sequence of uniform meshes, refining at each step by halving the mesh parameter $h$. The charts in Figure 1 plot the observed relative error against the mesh parameter $h$, for all possible parameter configurations. From even the briefest examination of the charts it is clear that as the mesh is refined, the relative error of the approximations is reduced. Whilst if we were to examine the gradients of the lines appearing in Figure 1, then they appear generally to be getting closer to 1, as the mesh is refined. With the gradients between the most refined mesh pairings having a mean value of 1.0301, across all configurations. This would suggest that our numerical scheme has order $\gamma\approx 1$, with the error in the approximations being $\mathcal{O}(h)$. The full numerical details of the errors and the associated convergence rates underlying Figure 1 may be found in [27, Appendix A].

8. Conclusions

In this article we introduced a numerical scheme for the approximate solution of a truncated version of a mixed discrete-continuous fragmentation model. The scheme was based upon a finite volume discretisation of the modelling equation for the continuous component.

The resulting numerical approximations were first shown to be nonnegative and to conserve mass, provided the underlying mesh satisfied certain constraints. Following which we established the weak convergence of a subsequence of our approximations, as the mesh size parameter $h$ was decreased to zero. The resulting limits were then shown to provide a weak solution to the truncated model.

By relating the scalar-valued weak formulation of our truncated model to an equivalent weak formulation within a Banach space setting, we were able to establish a one-to-one relationship between any scalar and Banach-space-valued weak solutions. Under suitable constraints, these Banach-space-valued weak solutions were shown to provide the unique strong solution to the associated abstract Cauchy problem, in the process establishing the uniqueness of the original scalar-valued weak solutions. Additionally, as a further consequence of this linkage, the scalar weak solutions were shown to be differentiable classical solutions.

Finally, by conducting a range of experiments with a test model, under varying model parameter choices and mesh refinements, we experimentally established that the error in our numerical solutions was $\mathcal{O}(h)$.

Acknowledgments

This work was supported by the UK Engineering and Physical Sciences Research Council [EP/J500495/1 03].