# A finite-volume scheme for dynamic reliability models

Christiane Cocozza-Thivent, Robert Eymard, Sophie Mercier

► To cite this version:

Christiane Cocozza-Thivent, Robert Eymard, Sophie Mercier. A finite-volume scheme for dynamic reliability models. IMA Journal of Numerical Analysis, 2005, 26 (3), pp.446-471. 10.1093/imanum/drl001 . hal-00693102

HAL Id: hal-00693102

https://hal.science/hal-00693102

Submitted on 31 Jul 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A finite-volume scheme for dynamic reliability models

C. COCOZZA-THIVENT†, R. EYMARD‡ AND S. MERCIER§

Laboratoire d'Analyse et de Mathématiques Appliquées (CNRS UMR 8050),
Université de Marne-la-Vallée, 5, boulevard Descartes, Champs-sur-Marne,
77454 Marne-la-Vallée Cedex 2, France

In a model arising in the dynamic reliability study of a system, the probability of the state of the system is completely described by the Chapman–Kolmogorov equations, which are scalar linear hyperbolic partial differential equations coupled by their right-hand side, the solution of which are probability measures. We propose in this paper a finite-volume scheme to approximate these measures. We show, thanks to the proof of the tightness of the approximate solution, that the conservation of the probability mass leads to a compactness property. The convergence of the scheme is then obtained in the space of continuous functions with respect to the time variable, valued in the set of probability measures on $\mathbb{R}^d$. We finally show on a numerical example the accuracy and efficiency of the approximation method.

Keywords: dynamic reliability; Markov processes; finite-volume method; weak convergence.

## 1. Introduction

We consider a system which is described, at each time $t \in \mathbb{R}_+$, by its physical state $I_t \in E$, where $E$ is a finite set, and by environmental variables $X_t \in \mathbb{R}^d$ (e.g. the pressure of some part of the system or the temperature). We assume that if the system remains in the same state for some time, then the environmental variables satisfy the ordinary differential equation

$$ \frac{dX_t}{dt} = \mathbf{v}(i, X_t) \quad \forall t \in \mathbb{R}_+ \text{ s.t., } \quad I_t = i \quad \forall i \in E, \qquad (1.1) $$

where $\mathbf{v}$ is a mapping from $E \times \mathbb{R}^d$ to $\mathbb{R}^d$. Stochastic transitions from state $(i, x) \in E \times \mathbb{R}^d$ to $(j, y) \in E \times \mathbb{R}^d$ are defined, thanks to a rate $a(i, j, x)$ and a probability measure $\mu(i, j, x)(dy)$. We denote by $\rho(t)(i, dx)$ the probability distribution of $(I_t, X_t)$, for all $t \in \mathbb{R}_+$. For such processes $(I_t, X_t)_{t \in \mathbb{R}_+}$ (called ‘piecewise deterministic Markov processes’ by Davis (1984, 1993)), the family of marginal distributions $(\rho(t)(i, dx))_{t \in \mathbb{R}_+}$ is solution of the Chapman–Kolmogorov equations

$$ \begin{aligned} & \int_0^t \sum_{i \in E} \int_{\mathbb{R}^d} \sum_{j \in E} a(i, j, x) \left( \int_{\mathbb{R}^d} \varphi(j, y) \mu(i, j, x)(dy) - \varphi(i, x) \right) \rho(s)(i, dx) ds \\ & + \int_0^t \sum_{i \in E} \int_{\mathbb{R}^d} \mathbf{v}(i, x) \cdot \nabla \varphi(i, x) \rho(s)(i, dx) ds - \sum_{i \in E} \int_{\mathbb{R}^d} \varphi(i, x) \rho(t)(i, dx) \\ & + \sum_{i \in E} \int_{\mathbb{R}^d} \varphi(i, x) \rho_0(i, dx) = 0 \quad \forall i \in E, \quad \forall t \in \mathbb{R}_+, \quad \forall \varphi \in C_c^1(\mathbb{R}^d)^E, \end{aligned} \qquad (1.2) $$

<footnote>†Email: christiane.cocozza@univ-mlv.fr</footnote>

<footnote>‡Email: robert.eymard@univ-mlv.fr</footnote>

<footnote>§Email: sophie.mercier@univ-mlv.fr</footnote>
where we denote by $C_c^1(\mathbb{R}^d)$ the set of continuously differentiable functions from $\mathbb{R}^d$ to $\mathbb{R}$ with a compact support, and, for a space of functions on some set $\Omega$, say $\mathcal{F}$, we denote by $\mathcal{F}^E$ the set of functions $f$ on $E \times \Omega$, such that for all $i \in E$, $f(i, \cdot) \in \mathcal{F}$. We assume the following hypothesis on the data, denoted by (H) in this paper.

1. The set $E$ is finite, $d \in \mathbb{N}^*$.

2. The transition rate $a$ is such that, for all $i, j \in E$, $a(i, j, \cdot) \in C_b(\mathbb{R}^d, \mathbb{R}_+)$, the set of continuous and bounded functions from $\mathbb{R}^d$ to $\mathbb{R}_+$. We denote by $A > 0$ an upper bound on these functions $a(i, j, \cdot)$. We denote by $b(i, \cdot) = \sum_{j \in E} a(i, j, \cdot)$, for all $i \in E$, and by $B = A \text{Card}E$.

3. The velocity $v$ is assumed to be such that, for all $i \in E$,

* the function $v(i, \cdot) \in C(\mathbb{R}^d, \mathbb{R}^d)$, the set of continuous functions from $\mathbb{R}^d$ to $\mathbb{R}^d$, and $v(i, \cdot)$ is locally Lipschitz continuous,

* the function $\text{div } v(i, \cdot)$ is almost everywhere bounded in absolute value by some real value $D > 0$ (independent of $i$),

* the function $v(i, \cdot)$ is sub-linear, i.e. there are some $V_1 > 0$ and $V_2 > 0$ such that

$$ \forall i \in E, \quad \forall x \in \mathbb{R}^d, \quad |v(i, x)| \le V_1|x| + V_2. $$

4. The function $\mu$ is such that

* for all $i, j \in E$ and $x \in \mathbb{R}^d$, $\mu(i, j, x) \in \mathcal{P}(\mathbb{R}^d)$, the set of probability measures on $\mathbb{R}^d$,

* for all $i, j \in E$ and all $\psi \in C_b(\mathbb{R}^d, \mathbb{R})$, the function $x \to \int \psi(y)\mu(i, j, x)(dy)$ is continuous from $\mathbb{R}^d$ to $\mathbb{R}$.

5. We assume that $\rho_0 \in \mathcal{P}(E \times \mathbb{R}^d)$, where $\mathcal{P}(E \times \mathbb{R}^d)$ is the set of probability measures on $E \times \mathbb{R}^d$.

In Cocozza-Thivent et al., we have proved that under (H), there exists one and only one solution $\rho(\cdot)$ to equation (1.2), namely, a single function $\rho(\cdot): \mathbb{R}_+ \mapsto \mathcal{P}(E \times \mathbb{R}^d)$ such that the function $t \mapsto \sum_{i \in E} \int_{\mathbb{R}^d} \psi(i, x)\rho(t)(i, dx)$ is Lebesgue-integrable and bounded on $[0, T]$, for all $T \in \mathbb{R}_+$ and for all $\psi \in C_b(\mathbb{R}^d)^E$ and such that (1.2) is fulfilled. Moreover, the unique solution $\rho(t)$ is actually such that the function $t \mapsto \sum_{i \in E} \int_{\mathbb{R}^d} \psi(i, x) \rho(t)(i, dx)$ is continuous from $\mathbb{R}_+$ to $\mathbb{R}$ for all $\psi \in C_b(\mathbb{R}^d)^E$.

**REMARK 1.1** Hypothesis (H-4) includes the important practical cases $\mu(i, j, x)(dy) = \delta_0(dy)$ (the Dirac mass in 0), which corresponds e.g. to the case of a semi-Markov process (also called Markov renewal process) where the elapsed time in the current state stands for the environmental condition (and is set equal to 0 at transition time) and $\mu(i, j, x)(dy) = \delta_x(dy)$, which corresponds to the continuity of the environmental variables or any combination of these two cases with respect to the coordinates of $x$.

Let us give two interesting cases where a functional interpretation of (1.2) can be given. In the case where there exist some regular functions $u_0 \in L^1(\mathbb{R}^d)^E$ and $u \in C(\mathbb{R}_+, L^1(\mathbb{R}^d))^E$ such that $\rho_0(i, dx) = u_0(i, x)dx$ and $\rho(t)(i, dx) = u(i, x, t)dx$, and assuming that $\mu(i, j, x)(dy) = M(i, j, x, y)dy$ for all $i, j \in E$ and $x, y \in \mathbb{R}^d$, these functions satisfy in a weak sense the following system of linear scalar
hyperbolic equations, which are only coupled by their right-hand side

$$
\begin{align}
\partial_t u(i, x, t) + \operatorname{div}(u(i, x, t)\mathbf{v}(i, x)) &= \sum_{j \in E} \int_{\mathbb{R}^d} a(j, i, y)u(j, y, t)M(j, i, y, x)dy - u(i, x, t)b(i, x), \nonumber \\
&\qquad \text{for a.e. } x \in \mathbb{R}^d, \quad \forall t \in \mathbb{R}_+, \forall i \in E, \tag{1.3} \\
u(i, x, 0) &= u_0(i, x), \quad \text{for a.e. } x \in \mathbb{R}^d, \quad \forall i \in E. \nonumber
\end{align}
$$

Let us now take the more particular case, where $d=1$; there exist regular functions $u_0 \in L^1(\mathbb{R}_+)^E$ and $u \in C(\mathbb{R}_+, L^1(\mathbb{R}_+))^E$ such that $\rho_0(i, dx) = u_0(i, x)dx$ and $\rho(t)(i, dx) = u(i, x, t)dx$, assuming that $\mu(i, j, x)(dy) = \delta_0(dy)$ for all $i, j \in E$ and $x \in \mathbb{R}_+$ (see Remark 1.1) and $\mathbf{v}(i, \cdot) \ge 0$ (in order to ensure the existence of a solution, the support of which remains a subset of $\mathbb{R}_+$). Then the term involving $\mu$ in (1.2) can be seen as resulting from the weak formulation of a boundary condition on $\{0\} \times \mathbb{R}_+$ (see Cocozza-Thivent & Eymard, 2004). Hence, the functions $u$ and $u_0$ satisfy in a weak sense the following system of linear scalar hyperbolic equations, only coupled by the boundary condition on $\{0\} \times \mathbb{R}_+$:

$$
\begin{align}
\partial_t u(i, x, t) + \operatorname{div}(u(i, x, t)\mathbf{v}(i, x)) &= -u(i, x, t)b(i, x), \quad \text{for a.e. } x \in \mathbb{R}_+, \forall t \in \mathbb{R}_+, \forall i \in E, \nonumber \\
\mathbf{v}(i, 0) u(i, 0, t) &= \sum_{j \in E} \int_{\mathbb{R}^d} a(j, i, y)u(j, y, t)dy, \quad \forall t \in \mathbb{R}_+, \forall i \in E, \tag{1.4} \\
u(i, x, 0) &= u_0(i, x), \quad \text{for a.e. } x \in \mathbb{R}^d, \forall i \in E. \nonumber
\end{align}
$$

Hence, in the particular case $E = \{1\}$, System (1.4) is completely identical, dropping the condition $b(1, \cdot) = a(1, 1, \cdot)$, to a population dynamics model called the McKendrick-Von Foerster model, improved by Sinko & Streifer (1967) and Bell & Anderson (1967). In this model, $u(1, x, t)$ represents the density of population of age $x$ at the time $t$, $\mathbf{v}(1, x)$ is the specific growth rate, $b(1, x)$ is the specific mortality at the age $x$, $a(1, 1, x)$ is the specific fertility at the age $x$ and $u_0(1, x)$ is the initial density of the population. This model allows for the computation of all the informations concerning the time evolution of the population and has been the object of numerous studies and adaptations to more complex situations (see Mischler et al., 2002, for non-linear improvements, and references therein, and see Feller, 1966, for the probabilistic notions).

For such models, whether dynamic reliability or population dynamics model, a challenging problem is to calculate the marginal distribution $\rho(t)(\cdot, dx)$ of the process $(I_t, X_t)$. Actually, as can be seen in the literature (see, e.g. Labeau, 1996), such a distribution is analytically calculable in only very simple cases and only numerical approximations are usually possible. One of the most up-to-date methods for such numerical computations appeals to Monte Carlo simulation (see, e.g. Labeau, 1996, and other papers by P. E. Labeau).

We construct here an approximation of $\rho(t)(i, dx)$ using a finite-volume method and we prove that such an approximation weakly converges towards the unique solution of (1.2). Such a method has already been studied and implemented in some restrictive cases of dynamic reliability (see Cocozza-Thivent & Eymard and Cocozza-Thivent & Eymard (2004) for the ‘semi-Markov’ case). Note that similar numerical approximations for (1.4) have been studied in the 1D case (see Ackleh et al. (2002), Ackleh & Ito (1997) and Abia et al. (2004), and references therein for the population dynamics problem). The study of the scheme on unstructured meshes in 2D (or more) is much more difficult, and the proofs of convergence for the approximation of a linear or non-linear scalar hyperbolic problem, using a finite-volume method, are more recent (they rely on the notion of weak bounded variation inequality, see, e.g. Eymard et al., 2000). Nevertheless, it appears that the suitable functional framework
for the study of (1.2) is the set of continuous functions of the time variable valued in $P(E \times \mathbb{R}^d)$, and not in $(L^1(\mathbb{R}^d) \cap L^\infty(\mathbb{R}^d))^E$. This difference makes it necessary to use some tools which are new in the framework of the study of convergence of numerical schemes for partial differential equations. Indeed, the relative compactness of a family of bounded functions of $L^\infty(\mathbb{R}^d)$ must be replaced by that of a family of probability measures in $P(\mathbb{R}^d)$, which implies the study of the tightness of this family (see Lemma 3.4). This study appears to be quite complex, due to the redistribution of probability induced by the measure $\mu(i, j, x)$ (the transport part, due to $v$, does not lead to an actual difficulty). Note also that the general convergence result is obtained in this paper under the assumptions $\delta t \to 0$ and $|M|/\delta t \to 0$, where $|M|$ is the space step and $\delta t$ is the time step. Such a condition instead of a simple bound on $|M|/\delta t$ seems to be necessary in the general setup, due to the apparent impossibility to prove in the general case a certain weak bounded variation inequality (see Remark 3.2). This creates a thorough difference with the $L^\infty$ framework, where the CFL condition for the convergence of the explicit scheme is expressed as a bound on $\delta t/|M|$.

Hence, the outline of this paper is the following. In Section 2, we give the numerical scheme which is deduced from a finite-volume approximation of (1.3). We then derive the convergence analysis from the conservation of the measure on the whole space and from a tightness estimate (see Section 3), which allow for some convergence properties, continuity properties and compactness properties, leading to the convergence of the scheme to the unique weak solution of (1.2). We then conclude this paper with a numerical example showing the efficiency and the accuracy of the method (Section 4).

## 2. The numerical scheme

We first denote by $g: E \times \mathbb{R}^d \times \mathbb{R} \to \mathbb{R}^d$ the solution of

$$ \frac{\partial g}{\partial t}(i, x, t) = v(i, g(i, x, t)) \quad \forall (i, x, t) \in E \times \mathbb{R}^d \times \mathbb{R}, \qquad (2.1) $$

with

$$ g(i, x, 0) = x \quad \forall (i, x) \in E \times \mathbb{R}^d. \qquad (2.2) $$

Thanks to the fact that the Cauchy-Lipschitz theorem holds within Hypothesis (H-3), the single solution $g$ of (2.1)-(2.2) fulfills the following properties:

1. the function $g$ is such that $g \in C^1(\mathbb{R}; C(\mathbb{R}^d))^E$, where $C^1(\mathbb{R}; C(\mathbb{R}^d))$ is the set of continuously differentiable functions with respect to time, valued in $C(\mathbb{R}^d)$,

2. the function $x \mapsto g(i, x, t)$ is locally Lipschitz continuous with respect to $x$ on $\mathbb{R}^d$, for all $t > 0$, all $i \in E$,

3. the following property holds:

$$ g(i, g(i, x, s), t) = g(i, x, s + t) \quad \forall (i, x, s, t) \in E \times \mathbb{R}^d \times \mathbb{R} \times \mathbb{R}. \qquad (2.3) $$

4. denoting by $J(i, y, s)$ the Jacobian in the change of space variable $t = s$ and $x = g(i, y, s)$, one has

$$ J(i, y, 0) = 1 \quad \forall (i, y) \in E \times \mathbb{R}^d $$
and

$$ \frac{\partial}{\partial s} J(i, y, s) = \operatorname{div} v(i, g(i, y, s)) J(i, y, s) \quad \forall (i, y, s) \in E \times \mathbb{R}^d \times \mathbb{R}, $$

which leads to

$$ J(i, y, s) = \exp\left(\int_{0}^{s} \operatorname{div} v(i, g(i, y, \tau)) d\tau\right) \quad \forall (i, y, s) \in E \times \mathbb{R}^d \times \mathbb{R}. $$

We then introduce the notion of admissible mesh $\mathcal{M}$ of $\mathbb{R}^d$, i.e. a family of measurable subsets of $\mathbb{R}^d$ such that

$$ 1. \bigcup_{K \in \mathcal{M}} K = \mathbb{R}^d, $$

$$ 2. \forall (K, L) \in \mathcal{M}^2, K \neq L \Rightarrow K \cap L = \emptyset (\mathcal{M} \text{ is a partition of } \mathbb{R}^d), $$

$$ 3. m_K := \int_K dx > 0 \text{ for all } K \in \mathcal{M}. $$

$$ 4. \sup_{K \in \mathcal{M}} \operatorname{diam}(K) < +\infty \text{ where } \operatorname{diam}(K) = \sup_{((x,y) \in K^2)} |x - y|. \text{ We then set } |\mathcal{M}| := \sup_{K \in \mathcal{M}} \operatorname{diam}(K). $$

Let a real number $\delta t > 0$ and an admissible mesh $\mathcal{M}$ be given. In order to define the numerical scheme, we first apply the principles of the finite-volume method to the initial condition, defining $p_0^{(i,K)}$ by

$$ m_K p_0^{(i,K)} = \int_K \rho_0(i, dx) \quad \forall K \in \mathcal{M} \quad \forall i \in E, \qquad (2.4) $$

and in order to determine the discrete rates of state transition, we set

$$ m_L a_{KL}^{(ij)} = \int_L a(j, i, y) \int_K \mu(j, i, y)(dx)dy \quad \forall (K, L) \in \mathcal{M}^2, \forall (i, j) \in E^2, \qquad (2.5) $$

$$ m_K b_K^{(i)} = \int_K b(i, x)dx \quad \forall K \in \mathcal{M}, \forall i \in E. $$

We define the following coefficients for the convection part of the scheme, using the notations given at the beginning of this section:

$$ m_{KL}^{(i)} = \int_{\{x \in K / g(i, x, -\delta t) \in L\}} J(i, x, -\delta t) dx = \int_{\{y \in L / g(i, y, \delta t) \in K\}} dy, \qquad (2.6) $$

which thus satisfy

$$ \sum_{K \in \mathcal{M}} m_{KL}^{(i)} = m_L \quad \forall L \in \mathcal{M}. \qquad (2.7) $$

For $n \in \mathbb{N}$ fixed, assuming $p_n^{(i,K)}$ to be constructed (all $i, K$), let us construct $p_{n+1}^{(i,K)}$. We first apply a finite-volume scheme to the left-hand side of the partial differential equations (1.3), which corresponds to the transport part of the equations. Hence, defining the function $P_n$ by

$$ P_n(i, x) = p_n^{(i,K)} \quad \forall x \in K, \forall K \in \mathcal{M}, \forall i \in E, \qquad (2.8) $$
we consider the weak solution $\tilde{u}$ of the partial differential equation

$$ \partial_t \tilde{u}(x, t) + \text{div}(\tilde{u}(x, t)v(i, x)) = 0, \quad \text{for a.e. } (x, t) \in \mathbb{R}^d \times (0, \delta t), \\ \tilde{u}(x, 0) = P_n(i, x), \quad \text{for a.e. } x \in \mathbb{R}^d. $$

This function $\tilde{u}$ is then such that $\tilde{u}(x, t) = J(i, x, -t)P_n(i, g(i, x, -t))$, for a.e. $x \in \mathbb{R}^d$ and all $t \in (0, \delta t)$. We then define $m_K \tilde{p}_n^{(i,K)} := \int_K \tilde{u}(x, \delta t)dx$, which thus satisfies, thanks to the above definitions,

$$ m_K \tilde{p}_n^{(i,K)} = \sum_{L \in \mathcal{M}} m_{KL}^{(i)} p_n^{(i,L)} \quad \forall K \in \mathcal{M}, \forall i \in E. \quad (2.9) $$

We now take into account the right-hand side of the partial differential equations (1.3), by setting

$$ m_K p_{n+1}^{(i,K)} = \frac{1}{\delta t b_K^{(i)} + 1} m_K \tilde{p}_n^{(i,K)} + \delta t \sum_{j \in E} \sum_{L \in \mathcal{M}} \frac{a_{KL}^{(ij)}}{\delta t b_L^{(j)} + 1} m_L \tilde{p}_n^{(j,L)} \quad \forall K \in \mathcal{M}, \forall i \in E. \quad (2.10) $$

Indeed, the first term of the right-hand side of (2.10) accounts for the probability that the state does not change during the time step, the denominator of this term corresponding to the negative part of the right-hand side of (1.3), i.e. the transition of the state $i$ to another state $j$. The second term of the right-hand side of (2.10) takes into account the transitions of all the other states $j$ to the state $i$. We then construct an approximation $P_t^{M, \delta t}(\cdot, x)dx$ of the measure $\rho(t)(\cdot, dx)$ such that $P_t^{M, \delta t}(\cdot, x)$ is constant on each $\{i\} \times [n\delta t, (n+1)\delta t[ \times K$:

$$ P_t^{M, \delta t}(i, x) := p_n^{(i,K)} = P_n(i, x) \quad \forall i \in E, \forall x \in K, \\ \forall t \in [n\delta t, (n+1)\delta t[, \forall K \in \mathcal{M}, \forall n \in \mathbb{N}. \quad (2.11) $$

**REMARK 2.1 (IMPLEMENTATION OF THE SCHEME)** For a practical implementation of the numerical scheme, the number of non-zero terms in $(p_n^{(i,K)})_{K \in \mathcal{M}}$ has to be finite at each time step $n$ ($n \in \mathbb{N}$). This can be ensured e.g. by the following assumptions: $\rho_0$ has compact support and $\inf_{K \in \mathcal{M}} m_K > 0$, which implies that the cardinality of $\{K \in \mathcal{M}: K \subset B(0, R)\}$ is finite for all $R > 0$, where $B(0, R) = \{x \in \mathbb{R}^d: |x| < R\}$. These additional assumptions are not required for the definition and the study of the numerical scheme but have to be kept in mind for its implementation.

**REMARK 2.2 (VARIABLE TIME STEPS)** For the sake of simplicity, we only consider in this paper a constant time step $\delta t$, although it would be easy to generalise our results to the case of a variable time step.

# 3. Convergence study

## 3.1 Estimates

We begin with a technical lemma, which will be useful in the following.

**LEMMA 3.1** Under Hypothesis (H-3), we have

$$ |g(i, x, t) - x| \le \left(|x| + \frac{V_2}{V_1}\right) (\exp(V_1t) - 1) \quad \forall x \in \mathbb{R}^d, \forall t \in \mathbb{R}_+, \forall i \in E. \quad (3.1) $$
Proof. Let $t \ge 0$ and $x \in \mathbb{R}^d$ be given. We can write

$$
\begin{aligned}
|g(i, x, t) - x| &= \left| \int_0^t \mathbf{v}(i, g(i, x, s)) \mathrm{d}s \right| \\
&\le \int_0^t (V_1(|g(i, x, s) - x| + |x|) + V_2) \mathrm{d}s \\
&= (V_1|x| + V_2)t + V_1 \int_0^t |g(i, x, s) - x| \mathrm{d}s.
\end{aligned}
$$

We set $h(t) = |g(i, x, t) - x|$. We then have $h(t) \le (V_1|x| + V_2)t + V_1 \int_0^t h(\tau)\mathrm{d}\tau$. Defining the continuous function $\tilde{h}$ by $\tilde{h}(t) = (V_1|x| + V_2)t + V_1 \int_0^t \tilde{h}(\tau)\mathrm{d}\tau$, we then get, using Gronwall's lemma, that $\tilde{h}(t) - h(t) \ge 0$ for all $t \ge 0$. Since $\tilde{h}(t) = \frac{V_1|x|+V_2}{V_1}(\exp(V_1t) - 1)$, we get (3.1). $\square$

LEMMA 3.2 (MASS CONSERVATION) Under Hypotheses (H), let $\mathcal{M}$ be an admissible mesh on $\mathbb{R}^d$ in the sense given in Section 2 and let $\delta t > 0$. Then there exists one and only one family of non-negative real values $(p_n^{(i,K)})_{n \in \mathbb{N}, i \in E, K \in \mathcal{M}}$ which is the solution of (2.4)–(2.10). Moreover, this family is such that

$$ \sum_{i \in E} \sum_{K \in \mathcal{M}} m_K p_n^{(i,K)} = 1 \quad \forall n \in \mathbb{N}. \qquad (3.2) $$

Proof. The existence and uniqueness of the solution $(p_n^{(i,K)})_{n \in \mathbb{N}, i \in E, K \in \mathcal{M}}$ to (2.4)–(2.10) is clear. We prove (3.2) by induction on $n$. The property is clear for $n = 0$. Suppose that $\sum_{i \in E} \sum_{K \in \mathcal{M}} m_K p_n^{(i,K)} = 1$. Then, $P_n, \tilde{p}_n^{(\cdot, \cdot)}$ and consequently $p_{n+1}^{(\cdot, \cdot)}$ are uniquely defined by (2.8)–(2.10).

Moreover, we have

$$ \sum_{i \in E} \sum_{K \in \mathcal{M}} m_K p_{n+1}^{(i,K)} = \sum_{i \in E} \sum_{K \in \mathcal{M}} \frac{1}{\delta t b_K^{(i)} + 1} m_K \tilde{p}_n^{(i,K)} + \delta t \sum_{j \in E} \sum_{L \in \mathcal{M}} \frac{m_L \tilde{p}_n^{(j,L)}}{\delta t b_L^{(j)} + 1} \sum_{i \in E} \sum_{K \in \mathcal{M}} a_{KL}^{(ij)}, $$

with

$$ \sum_{i \in E} \sum_{K \in \mathcal{M}} a_{KL}^{(ij)} = \frac{1}{m_L} \sum_{i \in E} \int_L a(j, i, y) \sum_{K \in \mathcal{M}} \int_K \mu(j, i, y)(dx)dy = b_L^{(j)}. $$

This leads to

$$
\begin{aligned}
\sum_{i \in E} \sum_{K \in \mathcal{M}} m_K p_{n+1}^{(i,K)} &= \sum_{i \in E} \sum_{K \in \mathcal{M}} \frac{1}{\delta t b_K^{(i)} + 1} m_K \tilde{p}_n^{(i,K)} + \delta t \sum_{j \in E} \sum_{L \in \mathcal{M}} \frac{b_L^{(j)} m_L \tilde{p}_n^{(j,L)}}{\delta t b_L^{(j)} + 1} \\
&= \sum_{i \in E} \sum_{K \in \mathcal{M}} m_K \tilde{p}_n^{(i,K)} = \sum_{i \in E} \sum_{K \in \mathcal{M}} \sum_{L \in \mathcal{M}} m_{KL}^{(i)} p_n^{(i,L)} = \sum_{i \in E} \sum_{L \in \mathcal{M}} m_L p_n^{(i,L)}, \quad (3.3)
\end{aligned}
$$

which completes the proof. $\square$
We now build a sequence of radii which permits us to control the propagation of mass by both mechanisms: transport and change of state. For this purpose, we first notice that for all $R > 0$ and all $i, j \in E$, the family $\{\mu(i, j, x)(dy): |x| \le R\}$ is weakly compact, which is a direct consequence of (H-4). The family $\{\mu(i, j, x)(dy)\}_{x \in B(0, R)}$ is then tight (see, e.g., Billingsley (1968)), which allows us to introduce the function $f_{i,j}: \mathbb{R}_+^* \times \mathbb{R}_+^* \to \mathbb{R}_+^*$, given by

$$ \forall R \in \mathbb{R}_+^*, \forall \varepsilon \in \mathbb{R}_+^*, \quad f_{i,j}(R, \varepsilon) = \inf \left\{ r \ge R, \sup_{x \in B(0,R)} \int_{\mathbb{R}^d \setminus B(0,r)} \mu(i, j, x)(dy) \le \varepsilon \right\}. \quad (3.4) $$

We get the following result.

**LEMMA 3.3 (PROPAGATION OF THE PROBABILITY MASS)** Under Hypotheses (H), let $\varepsilon > 0$, $T > 0$, $\theta > 0$ and $R_0 > 0$ be given. Let the sequences $(R^{(k)})_{k \in \mathbb{N}}$ and $(\overline{R}^{(k)})_{k \in \mathbb{N}}$ be defined by

$$ R^{(0)} = R_0, \quad (3.5) $$

$$ \overline{R}^{(k)} = \exp(V_1 T) R^{(k)} + \frac{V_2 + \theta \exp(V_1 T)}{V_1} (\exp(V_1 T) - 1) \quad \forall k \in \mathbb{N}, \quad (3.6) $$

$$ R^{(k+1)} = \sup_{i,j \in E} f_{i,j}(\overline{R}^{(k)} + \theta T, \varepsilon) \quad \forall k \in \mathbb{N}. \quad (3.7) $$

For all $\delta t \in (0, T)$, defining $N \in \mathbb{N}$ such that $N \delta t \le T < (N+1)\delta t$, let the sequence $(R_n^{(k)})_{n=0,\dots,N}$ be defined for all $k \in \mathbb{N}$ by

$$ R_0^{(k)} = R^{(k)}, $$

$$ R_{n+1}^{(k)} = (R_n^{(k)} + \theta \delta t) \exp(V_1 \delta t) + \frac{V_2}{V_1} (\exp(V_1 \delta t) - 1) \quad \forall n = 0, \dots, N-1. \quad (3.8) $$

Then the following properties hold:

$$ R_n^{(k)} = \exp(V_1 n \delta t) R^{(k)} + \left( \frac{V_2}{V_1} + \frac{\theta \delta t \exp(V_1 \delta t)}{\exp(V_1 \delta t) - 1} \right) (\exp(V_1 n \delta t) - 1) \quad \forall n = 0, \dots, N, \quad \forall k \in \mathbb{N}, \quad (3.9) $$

$$ R_n^{(k)} < R_{n+1}^{(k)}, \quad R^{(k)} \le R_n^{(k)} \le \overline{R}^{(k)} < R^{(k+1)} \le R_n^{(k+1)} \quad \forall n = 0, \dots, N, \quad \forall k \in \mathbb{N}, \quad (3.10) $$

$$ \sup_{x \in B(0, R_{n+1}^{(k-1)} + \theta \delta t), i, j \in E} \int_{\mathbb{R}^d \setminus B(0, R_{n+1}^{(k)})} \mu(i, j, x)(dy) \le \varepsilon \quad \forall n = 0, \dots, N-1, \quad \forall k \in \mathbb{N}^*, \quad (3.11) $$

$$ \begin{gathered} \forall k \in \mathbb{N}, \forall n = 0, \dots, N-1 \quad \forall i \in E, \forall x \in \mathbb{R}^d, \quad (3.12) \\ |x| < R_n^{(k)} + \theta \delta t \Rightarrow |g(i, x, \delta t)| < R_{n+1}^{(k)}. \end{gathered} $$

*Proof.* Using (3.8), we easily get the proof of (3.9) by induction on $n$. The proof of (3.10) is a consequence of $\exp(V_1 \delta t) - 1 \ge V_1 \delta t$ and therefore $\frac{\theta \delta t}{\exp(V_1 \delta t) - 1} \le \frac{\theta}{V_1}$ and of $f_{i,j}(\overline{R}^{(k)} + \theta T, \varepsilon) \ge \overline{R}^{(k)} + \theta T > \overline{R}^{(k)}$. Then, using that $f_{i,j}(R, \varepsilon)$ is nondecreasing with respect to $R$, the relations (3.10)
and (3.7) yield

$$
\begin{align*}
f_{i,j}(R_{n+1}^{(k-1)} + \theta \delta t, \varepsilon) &\le f_{i,j}(\bar{R}^{(k-1)} + \theta T, \varepsilon) \\
&\le R^{(k)} < R_{n+1}^{(k)} \quad \forall n = 0, \dots, N-1, \forall k \in \mathbb{N}^*, \forall i, j \in E
\end{align*}
$$

which gives (3.11) (note that it is important that the last inequality is strict since the infimum in (3.4) is
not necessarily reached). Finally, we get (3.12) using (3.8) and Lemma 3.1.
□

LEMMA 3.4 (TIGHTNESS) Under Hypotheses (H), let $\delta t > 0$ be given and let $\mathcal{M}$ be an admissible mesh on $\mathbb{R}^d$ in the sense given in Section 2. Let $\theta > |\mathcal{M}|/\delta t$ and $T > \delta t$. Let $(p_n^{(i,K)})_{i \in E, n \in \mathbb{N}, K \in \mathcal{M}}$ be the solution of (2.4)–(2.10).

Then, for all $\varepsilon_0 > 0$, there exists $R > 0$ only depending on $\varepsilon_0$, $T$, $V_1$, $V_2$, $A$, Card$E$, $\rho_0$, $\mu$ and $\theta$, such that

$$
\sum_{K \in \mathcal{M}} \sum_{K \cap B(0, R) \neq \emptyset} m_K p_n^{(i,K)} \geq 1 - \varepsilon_0 \quad \forall n \in \mathbb{N} \text{ s.t. } n\delta t \leq T, \quad (3.13)
$$

whence the tightness of the family of probability measures $(P_t^{\mathcal{M},\delta t}(\cdot, x)dx)_{\{(t,\mathcal{M},\delta t):|\mathcal{M}|/\delta t<\theta;0<t,\delta t<T\}}$ on $E \times \mathbb{R}^d$.

Proof. Let $\varepsilon > 0$ be given (this real number will be chosen as a function of $\varepsilon_0$ at the end of the proof).
We denote by $R_0 > 0$ a real number such that

$$
\sum_{i \in E} \int_{\mathbb{R}^d \setminus B(0, R_0)} \rho_0(i, dx) \le \varepsilon. \quad (3.14)
$$

Hence, $R_0$ only depends on $\varepsilon$ and $\rho_0$. We then denote by $(R^{(k)})_{k \in \mathbb{N}}$ and $(\bar{R}^{(k)})_{k \in \mathbb{N}}$ the sequences defined for these values $\varepsilon > 0$ and $R_0 > 0$ by (3.5), (3.6) and (3.7), which therefore only depend on $\varepsilon$, $\rho_0$, $T$, $V_1$, $V_2$, $\mu$ and $\theta$. We finally denote by $N \in \mathbb{N}$ the integer, such that $N\delta t \le T < (N+1)\delta t$, and by $(R_n^{(k)})_{n=0,\dots,N}$ the sequence defined for all $k \in \mathbb{N}$ by (3.8). For all $n=0,\dots,N$ and $k \in \mathbb{N}$, we denote by $\mathcal{M}_n^{(k)} = \{K \in \mathcal{M}, K \cap B(0, R_n^{(k)}) = \emptyset\}$, and we define

$$
q_n^{(k)} = \sum_{K \in \mathcal{M}_n^{(k)}} \sum_{i \in E} m_K p_n^{(i,K)} \quad \forall n = 0, \dots, N, \quad \forall k \in \mathbb{N}.
$$

which leads, using (3.14), to

$$
q_0^{(k)} \le \varepsilon \quad \forall k \in \mathbb{N}. \tag{3.15}
$$

The principle of this proof is to find some $k_0 \in \mathbb{N}$ and some $C > 0$, independent of $\mathcal{M}$, $\delta t$ and $\varepsilon$,
such that $q_n^{(k_0)} \le C\varepsilon$, for all $n \in \mathbb{N}$ such that $n\delta t \le T$, which will be sufficient to conclude, choosing
$\varepsilon = \varepsilon_0/C$ and $R \ge \bar{R}^{(k_0)}$ for all $n \in \mathbb{N}$ such that $n\delta t \le T$. Let $k \in \mathbb{N}^*$ and $n = 0, \dots, N-1$ be given.
We have

$$
q_{n+1}^{(k)} = \sum_{K \in \mathcal{M}_{n+1}^{(k)}} \sum_{i \in E} \left( \frac{1}{\delta t b_K^{(i)} + 1} m_K \tilde{p}_n^{(i,K)} + \delta t \sum_{j \in E} \sum_{L \in \mathcal{M}} \frac{a_{KL}^{(ij)}}{\delta t b_L^{(j)} + 1} m_L \tilde{p}_n^{(j,L)} \right).
$$
We thus get

$$
\begin{equation}
\begin{aligned}
q_{n+1}^{(k)} = & \sum_{K \in \mathcal{M}_{n+1}^{(k)}} \sum_{i \in E} \frac{1}{\delta t b_K^{(i)} + 1} m_K \tilde{p}_n^{(i,K)} + \delta t \sum_{L \in \mathcal{M}_{n+1}^{(k-1)}} \sum_{j \in E} \sum_{K \in \mathcal{M}_{n+1}^{(k)}} \sum_{i \in E} \frac{a_{KL}^{(ij)}}{\delta t b_L^{(j)} + 1} m_L \tilde{p}_n^{(j,L)} \\
& + \delta t \sum_{L \in \mathcal{M} \setminus \mathcal{M}_{n+1}^{(k-1)}} \sum_{j \in E} \sum_{K \in \mathcal{M}_{n+1}^{(k)}} \sum_{i \in E} \frac{a_{KL}^{(ij)}}{\delta t b_L^{(j)} + 1} m_L \tilde{p}_n^{(j,L)}.
\end{aligned}
\tag{3.16}
\end{equation}
$$

We first look for a bound for the third term of the right-hand side of (3.16). Let $j \in E$ and $L \in \mathcal{M} \setminus \mathcal{M}_{n+1}^{(k-1)}$ be given. We have

$$
\sum_{i \in E} \sum_{K \in \mathcal{M}_{n+1}^{(k)}} a_{KL}^{(ij)} = \sum_{i \in E} \frac{1}{m_L} \int_L a(j, i, y) \sum_{K \in \mathcal{M}_{n+1}^{(k)}} \int_K \mu(j, i, y)(dx)dy. \quad (3.17)
$$

Since we have $L \cap B(0, R_{n+1}^{(k-1)}) \neq \emptyset$ and $\text{diam}(L) \le |\mathcal{M}| < \theta \delta t$, thus we have $y \in B(0, R_{n+1}^{(k-1)} + \theta \delta t)$.
Since $K \subset \mathbb{R}^d \setminus B(0, R_{n+1}^{(k)})$ we can then write, using (3.11),

$$
\forall y \in L \quad \sum_{K \in \mathcal{M}_{n+1}^{(k)}} \int_K \mu(j, i, y)(dx) \le \int_{\mathbb{R}^d \setminus B(0, R_{n+1}^{(k)})} \mu(j, i, y)(dx) \le \varepsilon. \quad (3.18)
$$

Thus, we get from (3.17) and (3.18) (recall that $B = A \text{Card}E$)

$$
\forall n = 0, \dots, N-1, \forall k \in \mathbb{N}^*, \forall j \in E, \forall L \in \mathcal{M} \setminus \mathcal{M}_{n+1}^{(k-1)}, \sum_{i \in E} \sum_{K \in \mathcal{M}_{n+1}^{(k)}} a_{KL}^{(ij)} \le B\varepsilon. \quad (3.19)
$$

Since we get a bound of the second term at the right-hand side of (3.16), using

$$
\sum_{K \in \mathcal{M}_{n+1}^{(k)}} \sum_{i \in E} a_{KL}^{(ij)} \le b_L^{(j)} \quad \forall j \in E, \quad \forall L \in \mathcal{M},
$$

we thus get from (3.16) and (3.19)

$$
q_{n+1}^{(k)} \le \sum_{K \in \mathcal{M}_{n+1}^{(k)}} \sum_{i \in E} \frac{1}{\delta t b_K^{(i)} + 1} m_K \tilde{p}_n^{(i,K)} + \delta t \sum_{L \in \mathcal{M}_{n+1}^{(k-1)}} \sum_{j \in E} \frac{b_L^{(j)}}{\delta t b_L^{(j)} + 1} m_L \tilde{p}_n^{(j,L)} + \delta t B\varepsilon,
$$

using $\sum_{L \in \mathcal{M}} \sum_{j \in E} m_L \tilde{p}_n^{(j,L)} = 1$ (see (3.3)), which yields

$$
q_{n+1}^{(k)} \le \sum_{K \in \mathcal{M}_{n+1}^{(k)}} \sum_{i \in E} m_K \tilde{p}_n^{(i,K)} + \delta t \sum_{L \in \mathcal{M}_{n+1}^{(k-1)} \setminus \mathcal{M}_{n+1}^{(k)}} \sum_{j \in E} \frac{b_L^{(j)}}{\delta t b_L^{(j)} + 1} m_L \tilde{p}_n^{(j,L)} + \delta t B\varepsilon.
$$

Therefore, since $b_L^{(j)} \le B$, and thus $\frac{b_L^{(j)}}{\delta t b_L^{(j)} + 1} \le \frac{B}{\delta t B + 1}$, we get

$$
q_{n+1}^{(k)} \le \sum_{K \in \mathcal{M}_{n+1}^{(k)}} \sum_{i \in E} m_K \tilde{p}_n^{(i,K)} + \frac{\delta t B}{\delta t B + 1} \sum_{L \in \mathcal{M}_{n+1}^{(k-1)} \setminus \mathcal{M}_{n+1}^{(k)}} \sum_{j \in E} m_L \tilde{p}_n^{(j,L)} + \delta t B \varepsilon.
$$
Since $R_{n+1}^{(k-1)} \le R_{n+1}^{(k)}$ implies $M_{n+1}^{(k)} \subset M_{n+1}^{(k-1)}$, we get

$$q_{n+1}^{(k)} \le \frac{1}{\delta t B + 1} \sum_{K \in \mathcal{M}_{n+1}^{(k)}} \sum_{i \in E} m_K \tilde{p}_n^{(i,K)} + \frac{\delta t B}{\delta t B + 1} \sum_{L \in \mathcal{M}_{n+1}^{(k-1)}} \sum_{j \in E} m_L \tilde{p}_n^{(j,L)} + \delta t B \varepsilon. \quad (3.20)$$

Let us now remark that, thanks to the relation (3.12), we obtain

$$\forall k \in \mathbb{N}, \forall n = 0, \dots, N-1, \forall i \in E, \forall K \in \mathcal{M}_{n+1}^{(k)}, \forall L \in \mathcal{M}, m_{KL}^{(i)} \neq 0 \Rightarrow L \in \mathcal{M}_n^{(k)}.$$

This implies using (2.9) and (2.7)

$$\sum_{i \in E} \sum_{K \in \mathcal{M}_{n+1}^{(k)}} m_K \tilde{p}_n^{(i,K)} \le q_n^{(k)} \quad \forall n = 0, \dots, N, \quad \forall k \in \mathbb{N}.$$

Hence, we get from (3.20)

$$q_{n+1}^{(k)} \le \frac{1}{\delta t B + 1} q_n^{(k)} + \frac{\delta t B}{\delta t B + 1} q_n^{(k-1)} + \delta t B \varepsilon.$$

Since $q_0^{(k)} \le \varepsilon$ for all $k \in \mathbb{N}$ and $q_n^{(0)} \le 1$ for all $n \in \mathbb{N}$, we then easily prove by induction on $n$ and on $k$ that, for all $k, n \in \mathbb{N}$, $q_n^{(k)} \le \varepsilon + \tilde{q}_n^{(k)} + n\delta t B \varepsilon$, where $(\tilde{q}_n^{(k)})_{n,k \in \mathbb{N}}$ is defined by

$$(\tilde{q}_n^{(0)} = 1 \quad \forall n \in \mathbb{N}) \quad \text{and} \quad (\tilde{q}_0^{(k)} = 0 \quad \forall k \in \mathbb{N}^*) \quad \text{and}$$

$$(\tilde{q}_{n+1}^{(k)} = \tilde{q}_n^{(k)} + \delta t B \tilde{q}_n^{(k-1)} \quad \forall n = 0, \dots, N, \forall k \in \mathbb{N}^*).$$

Then an easy inductive argument shows that, for all $n = 0, \dots, N$, we have $\tilde{q}_n^{(k)} = \frac{n!}{k!(n-k)!} (\delta t B)^k$ for all $k = 0, \dots, n$ (with the usual convention $0! = 1$) and $\tilde{q}_n^{(k)} = 0$ for all $k > n$. This gives

$$\tilde{q}_n^{(k)} \le \frac{n^k}{k!} (\delta t B)^k \le \frac{(BT)^k}{k!} \quad \forall n \in \mathbb{N} \text{ s.t. } n\delta t \le T, \forall k \in \mathbb{N}.$$

We thus obtain

$$q_n^{(k)} \le \varepsilon + \frac{(BT)^k}{k!} + B\varepsilon T \quad \forall n \in \mathbb{N} \text{ s.t. } n\delta t \le T, \forall k \in \mathbb{N}.$$

We then choose $k_0 \in \mathbb{N}$ such that $\frac{(BT)^{k_0}}{k_0!} \le \varepsilon$. Thus, $k_0$ only depends on $\varepsilon$, $B$ (with $B = A \text{Card} E$) and $T$. We then get

$$q_n^{(k_0)} \le \varepsilon(2 + BT) \quad \forall n \in \mathbb{N} \text{ s.t. } n\delta t \le T.$$

The conclusion of the above lemma follows by setting $R = \bar{R}^{(k_0)}$ (defined by (3.6)) and $\varepsilon = \varepsilon_0/(2 + BT)$. $\square$

## 3.2 Convergence lemmas

Let $f \in C_b(\mathbb{R}^d)^E$ be a given function, and let $t_1, t_2 \in \mathbb{R}_+$ be such that $0 \le t_1 \le t_2$. For a given admissible mesh in the sense given in Section 2, and for a given $\delta t > 0$, using all the notations of
Section 2, we define the expressions $T_1^{\mathcal{M},\delta t}(f, t_1, t_2)$ and $T_2^{\mathcal{M},\delta t}(f, t_1, t_2)$ by

$$T_1^{\mathcal{M},\delta t}(f, t_1, t_2) = \delta t \sum_{i \in E} \sum_{n=n_1}^{n_2-1} \sum_{K \in \mathcal{M}} m_K \bar{f}^{(i,K)} \frac{\tilde{p}_n^{(i,K)}}{\delta t b_K^{(i)} + 1}, \quad (3.21)$$

$$T_2^{\mathcal{M},\delta t}(f, t_1, t_2) = \sum_{i \in E} \int_{[t_1, t_2]} \int_{\mathbb{R}^d} f(i, x) P_t^{\mathcal{M},\delta t}(i, x) dx dt, \quad (3.22)$$

where

$$\bar{f}^{(i,K)} = \frac{1}{m_K} \int_K f(i, x) dx \quad \forall i \in E, \forall K \in \mathcal{M}, \quad (3.23)$$

represents the mean value of $f$ on $K$ and where we denote by $n_1, n_2 \in \mathbb{N}$ the integers such that $n_1 \delta t \le t_1 < (n_1 + 1) \delta t$ and $n_2 \delta t \le t_2 < (n_2 + 1) \delta t$.

REMARK 3.1 Throughout this paper, we use the convention that, if $n_1 \ge n_2$, all the sums $\sum_{n=n_1}^{n_2-1} \cdot$) are set equal to zero.

Similarly, for $\varphi \in C_c^1(\mathbb{R}^d)^E$, we define the expressions $T_3^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$, $T_4^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$, $T_5^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$, $T_6^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$, $T_7^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$ and $T_8^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$ by

$$T_3^{\mathcal{M},\delta t}(\varphi, t_1, t_2) = \sum_{i \in E} \sum_{n=n_1}^{n_2-1} \sum_{K \in \mathcal{M}} m_K \tilde{p}_n^{(i,K)} \frac{\delta t b_K^{(i)}}{\delta t b_K^{(i)} + 1} \bar{\varphi}^{(i,K)}, \quad (3.24)$$

$$T_4^{\mathcal{M},\delta t}(\varphi, t_1, t_2) = \sum_{i \in E} \int_{t_1}^{t_2} \int_{\mathbb{R}^d} b(i, x) \varphi(i, x) P_t^{\mathcal{M},\delta t}(i, x) dx dt, \quad (3.25)$$

$$T_5^{\mathcal{M},\delta t}(\varphi, t_1, t_2) = \sum_{i \in E} \sum_{n=n_1}^{n_2-1} \sum_{K \in \mathcal{M}} \delta t \sum_{j \in E} \sum_{L \in \mathcal{M}} a_{KL}^{(ij)} \frac{m_L \tilde{p}_n^{(j,L)}}{\delta t b_L^{(j)} + 1} \bar{\varphi}^{(i,K)}, \quad (3.26)$$

$$T_6^{\mathcal{M},\delta t}(\varphi, t_1, t_2) = \sum_{i \in E} \int_{t_1}^{t_2} \int_{\mathbb{R}^d} \sum_{j \in E} a(i, j, x) \left( \int \varphi(j, y) \mu(i, j, x)(dy) \right) P_t^{\mathcal{M},\delta t}(i, x) dx dt, \quad (3.27)$$

$$T_7^{\mathcal{M},\delta t}(\varphi, t_1, t_2) = \sum_{i \in E} \sum_{n=n_1}^{n_2-1} \sum_{K \in \mathcal{M}} [m_K p_{n+1}^{(i,K)} - m_K \tilde{p}_n^{(i,K)}] \bar{\varphi}^{(i,K)}, \quad (3.28)$$

$$T_8^{\mathcal{M},\delta t}(\varphi, t_1, t_2) = \sum_{i \in E} \left[ - \int_{t_1}^{t_2} \int_{\mathbb{R}^d} \mathbf{v}(i, x) \cdot \nabla \varphi(i, x) P_t^{\mathcal{M},\delta t}(i, x) dx dt - \int_{\mathbb{R}^d} \varphi(i, x) P_{t_1}^{\mathcal{M},\delta t}(i, x) dx + \int_{\mathbb{R}^d} \varphi(i, x) P_{t_2}^{\mathcal{M},\delta t}(i, x) dx \right]. \quad (3.29)$$

Note that considering Scheme (2.10) multiplied by $\bar{\varphi}^{(i,K)}$ and summed over $i \in E, K \in \mathcal{M}$ and $n = n_1, \dots, n_2 - 1$, we exactly get

$$-T_3^{\mathcal{M},\delta t}(\varphi, t_1, t_2) + T_5^{\mathcal{M},\delta t}(\varphi, t_1, t_2) - T_7^{\mathcal{M},\delta t}(\varphi, t_1, t_2) = 0. \quad (3.30)$$
On the other hand, using the above definitions, we can write

$$
\begin{align}
& -T_4^{\mathcal{M},\delta t}(\varphi, t_1, t_2) + T_6^{\mathcal{M},\delta t}(\varphi, t_1, t_2) - T_8^{\mathcal{M},\delta t}(\varphi, t_1, t_2) \nonumber \\
&= \sum_{i \in E} \left[ \int_{t_1}^{t_2} \int_{\mathbb{R}^d} \left[ \sum_{j \in E} a(i, j, x) \int \varphi(j, y) \mu(i, j, x) (dy) \right. \right. \nonumber \\
&\qquad \left. \left. -b(i, x) \varphi(i, x) + \mathbf{v}(i, x) \cdot \nabla \varphi(i, x) \right] P_t^{\mathcal{M},\delta t}(i, x) dx dt \right. \nonumber \\
&\qquad \left. + \int_{\mathbb{R}^d} \varphi(i, x) P_{t_1}^{\mathcal{M},\delta t}(i, x) dx - \int_{\mathbb{R}^d} \varphi(i, x) P_{t_2}^{\mathcal{M},\delta t}(i, x) dx \right]. \tag{3.31}
\end{align}
$$

The following results will lead to the convergence of the scheme, by proving that the differences be-
tween $T_3^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$ and $T_4^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$, between $T_5^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$ and $T_6^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$ and between
$T_7^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$ and $T_8^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$ can be made as small as desired, letting $\delta t \to 0$ and $|\mathcal{M}|/\delta t \to 0$.

LEMMA 3.5 Under assumptions (H), let $\theta > 0$, $f \in C_b(\mathbb{R}^d)^E$, $T > 0$ and $t_1, t_2 \in \mathbb{R}_+$ be such that $0 \le t_1 \le t_2 \le T$. Then, for each $\varepsilon_0 > 0$, there exists $\delta t_0$, which only depends on $\varepsilon_0$, $T$, $f$, $V_1$, $V_2$, $A$, Card $E$, $\rho_0$, $\mu$ and $\theta$, such that for all $\delta t \in [0, \delta t_0]$, for all admissible meshes $\mathcal{M}$ in the sense given in Section 2, such that $|\mathcal{M}| \le \theta \delta t$, with $P^{\mathcal{M},\delta t}$ given by (2.4–2.11) and $T_1^{\mathcal{M},\delta t}(f, t_1, t_2)$ and $T_2^{\mathcal{M},\delta t}(f, t_1, t_2)$ defined, respectively, by (3.21) and (3.22), we have

$$
|T_1^{\mathcal{M},\delta t}(f, t_1, t_2) - T_2^{\mathcal{M},\delta t}(f, t_1, t_2)| \le \varepsilon_0, \quad (3.32)
$$

which means precisely that $\lim_{\delta t \to 0, |\mathcal{M}| \le \theta \delta t} (T_1^{\mathcal{M},\delta t}(f, t_1, t_2) - T_2^{\mathcal{M},\delta t}(f, t_1, t_2)) = 0$.

*Proof.* Let $\varepsilon > 0$. For a given admissible mesh in the sense given in Section 2, and for a given $\delta t > 0$, using all the notations of Section 2, we define the expression $T_9^{\mathcal{M},\delta t}(f, t_1, t_2)$ by

$$
T_9^{\mathcal{M},\delta t}(f, t_1, t_2) = \delta t \sum_{i \in E} \sum_{n=n_1}^{n_2-1} \sum_{K \in \mathcal{M}} m_K \bar{f}^{(i,K)} \tilde{p}_n^{(i,K)}.
$$

We then have

$$
\begin{align*}
|T_1^{\mathcal{M},\delta t}(f, t_1, t_2) - T_9^{\mathcal{M},\delta t}(f, t_1, t_2)| &\le \delta t B \delta t \sum_{i \in E} \sum_{n=n_1}^{n_2-1} \sum_{K \in \mathcal{M}} m_K |\bar{f}^{(i,K)}| \tilde{p}_n^{(i,K)} \\
&\le \delta t B T \sup_{i \in E} \|f(i, \cdot)\|_{L^\infty(\mathbb{R}^d)},
\end{align*}
$$

and therefore there exists $\delta t_1 > 0$ only depending on $\varepsilon$, $f$, $B$ and $T$ such that if $\delta t \le \delta t_1$, then
$|T_1^{\mathcal{M},\delta t}(f, t_1, t_2) - T_9^{\mathcal{M},\delta t}(f, t_1, t_2)| \le \varepsilon$.
We now consider $T_{10}^{\mathcal{M},\delta t}(f, t_1, t_2)$ defined by

$$T_{10}^{\mathcal{M},\delta t}(f, t_1, t_2) = \sum_{i \in E} \int_{n_1 \delta t}^{n_2 \delta t} \int_{\mathbb{R}^d} f(i, x) P_t^{\mathcal{M},\delta t}(i, x) dx dt = \delta t \sum_{i \in E} \sum_{n=n_1}^{n_2-1} \sum_{K \in \mathcal{M}} m_K \bar{f}^{(i,K)} p_n^{(i,K)}.$$

We then get that

$$|T_{10}^{\mathcal{M},\delta t}(f, t_1, t_2) - T_2^{\mathcal{M},\delta t}(f, t_1, t_2)| \le \sum_{i \in E} \int_{\mathbb{R}^d} \left( \int_{n_1 \delta t}^{(n_1+1)\delta t} |f(i, x)| P_t^{\mathcal{M},\delta t}(i, x) dt + \int_{n_2 \delta t}^{(n_2+1)\delta t} |f(i, x)| P_t^{\mathcal{M},\delta t}(i, x) dt \right) dx,$$

and therefore $|T_{10}^{\mathcal{M},\delta t}(f, t_1, t_2) - T_2^{\mathcal{M},\delta t}(f, t_1, t_2)| \le 2\delta t \sup_{i \in E} \|f(i, \cdot)\|_{L^{\infty}(\mathbb{R}^d)}$. Therefore, there exists $\delta t_2$, only depending on $\varepsilon$ and $f$, such that $|T_{10}^{\mathcal{M},\delta t}(f, t_1, t_2) - T_2^{\mathcal{M},\delta t}(f, t_1, t_2)| \le \varepsilon$.

Setting $T_{11}^{\mathcal{M},\delta t}(f, t_1, t_2) = T_9^{\mathcal{M},\delta t}(f, t_1, t_2) - T_{10}^{\mathcal{M},\delta t}(f, t_1, t_2)$, we get

$$T_{11}^{\mathcal{M},\delta t}(f, t_1, t_2) = \delta t \sum_{i \in E} \sum_{n=n_1}^{n_2-1} \sum_{K \in \mathcal{M}} p_n^{(i,K)} \sum_{L \in \mathcal{M}} m_{LK}^{(i)} (\bar{f}^{(i,L)} - \bar{f}^{(i,K)}).$$

We then get that, for a given $R > 0$, $T_{11}^{\mathcal{M},\delta t}(f, t_1, t_2) = T_{12}^{\mathcal{M},\delta t,R}(f, t_1, t_2) + T_{13}^{\mathcal{M},\delta t,R}(f, t_1, t_2)$, with

$$T_{12}^{\mathcal{M},\delta t,R}(f, t_1, t_2) = \delta t \sum_{i \in E} \sum_{n=n_1}^{n_2-1} \sum_{K \in \mathcal{M} \cap B(0,R)=\emptyset} p_n^{(i,K)} \sum_{L \in \mathcal{M}} m_{LK}^{(i)} (\bar{f}^{(i,L)} - \bar{f}^{(i,K)})$$

and

$$T_{13}^{\mathcal{M},\delta t,R}(f, t_1, t_2) = \delta t \sum_{i \in E} \sum_{n=n_1}^{n_2-1} \sum_{K \in \mathcal{M} \cap B(0,R) \neq \emptyset} p_n^{(i,K)} \sum_{L \in \mathcal{M}} m_{LK}^{(i)} (\bar{f}^{(i,L)} - \bar{f}^{(i,K)}).$$

Since

$$|T_{12}^{\mathcal{M},\delta t,R}(f, t_1, t_2)(f, \delta t, \mathcal{M}, R)| \le 2 \sup_{i \in E} \|f(i, \cdot)\|_{L^{\infty}(\mathbb{R}^d)} \delta t \sum_{i \in E} \sum_{n=n_1}^{n_2-1} \sum_{K \in \mathcal{M} \cap B(0,R)=\emptyset} p_n^{(i,K)} m_K,$$

thanks to Lemma 3.4, letting $\delta t \le T$ and therefore $|\mathcal{M}| \le \theta T$, it is possible to choose an $R$ which only depends on $\varepsilon, T, V_1, V_2, A, \text{Card}E, \rho_0, \mu, f$ and $\theta$ such that $|T_{12}^{\mathcal{M},\delta t,R}(f, t_1, t_2)| \le \varepsilon T$.

Now let $K \in \mathcal{M}$ be such that $K \cap B(0, R) \neq \emptyset$ and $L$ such that $m_{LK}^{(i)} \neq 0$. Thanks to Lemma 3.1, we can choose $R' \ge R + \theta T \ge R + \text{diam}(K)$, only depending on $T, V_1, V_2, \theta$ and $R$ such that $L \subset B(0, R')$ and $K \subset B(0, R')$. Since $f$ is continuous, it is uniformly continuous on $B(0, R')$. Let $\alpha \in (0, 1)$ be such that $|f(i, x) - f(i, y)| \le \varepsilon$ for all $(i, x, y) \in E \times B(0, R') \times B(0, R')$ with $|x - y| \le \alpha$.
Thanks again to Lemma 3.1, let $\delta t_3 > 0$, which only depends on $\alpha, \theta, V_1, V_2$ and $T$, be such that, for all $\delta t \le \delta t_3$, we have $|g(i, z, \delta t) - z| \le \alpha/3$ for all $z \in B(0, R')$.

As $m_{LK}^{(i)} \ne 0$, there is a $z \in K$ such that $g(i, z, \delta t) \in L$. For $\delta t \le \min(\delta t_3, \frac{\alpha}{3\theta})$, we deduce that, for all $x \in K$ and $y \in L$,

$$|x - y| \le |x - z| + |z - g(i, z, \delta t)| + |g(i, z, \delta t) - y| \le \alpha$$

and consequently $|\bar{f}^{(i,L)} - \bar{f}^{(i,K)}| \le \varepsilon$. This proves that $|T_{13}^{M,\delta t,R}(f, t_1, t_2)| \le \varepsilon T$.

Choosing at the beginning of the proof $\varepsilon$ such that $\varepsilon + \varepsilon + \varepsilon T + \varepsilon T = \varepsilon_0$, and then taking $\delta t_0 = \min(\delta t_1, \delta t_2, T, \delta t_3, \frac{\alpha}{3\theta})$ completes the proof. $\square$

We now deduce the following result.

LEMMA 3.6 Under assumptions (H), let $\theta > 0$, $\varphi \in C_c^1(\mathbb{R}^d)^E$, $T > 0$ and $t_1, t_2 \in \mathbb{R}_+$ be such that $0 \le t_1 \le t_2 \le T$. Then, for each $\varepsilon_0 > 0$, there exists $\delta t_0$, which only depends on $\varepsilon_0$, $T$, $\varphi$, $V_1$, $V_2$, $A$, Card$E$, $\rho_0$, $\mu$ and $\theta$, such that for all $\delta t \in [0, \delta t_0]$, for all admissible meshes $\mathcal{M}$ in the sense given in Section 2, such that $|\mathcal{M}| \le \theta \delta t$, with $P^{\mathcal{M},\delta t}$ be given by (2.4)–(2.11) and $T_3^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$, $T_4^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$, $T_5^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$ and $T_6^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$ defined by (3.24)–(3.27), we have

$$|T_3^{\mathcal{M},\delta t}(\varphi, t_1, t_2) - T_4^{\mathcal{M},\delta t}(\varphi, t_1, t_2)| \le \varepsilon_0, \quad (3.33)$$

and

$$|T_5^{\mathcal{M},\delta t}(\varphi, t_1, t_2) - T_6^{\mathcal{M},\delta t}(\varphi, t_1, t_2)| \le \varepsilon_0, \quad (3.34)$$

which means precisely that $\lim_{\delta t \to 0, |\mathcal{M}| \le \theta \delta t} (T_3^{\mathcal{M},\delta t}(f, t_1, t_2) - T_4^{\mathcal{M},\delta t}(f, t_1, t_2)) = 0$ and $\lim_{\delta t \to 0, |\mathcal{M}| \le \theta \delta t} (T_5^{\mathcal{M},\delta t}(f, t_1, t_2) - T_6^{\mathcal{M},\delta t}(f, t_1, t_2)) = 0$.

*Proof.* We consider the function $f$ defined by

$$f(i, x) = b(i, x)\varphi(i, x) \quad \forall (i, x) \in E \times \mathbb{R}^d.$$

Setting $C_\varphi = \sup_{i \in E} \|\nabla \varphi(i, \cdot)\|_{L^\infty(\mathbb{R}^d)}$, we then have

$$|b_K^{(i)} \bar{\varphi}^{(i,K)} - \bar{b} \varphi^{(i,K)}| \le C_\varphi |\mathcal{M}| B,$$

which yields that

$$|T_3^{\mathcal{M},\delta t}(\varphi, t_1, t_2) - T_1^{\mathcal{M},\delta t}(f, t_1, t_2)| \le C_\varphi |\mathcal{M}| B T.$$

Noting that $T_2^{\mathcal{M},\delta t}(f, t_1, t_2) = T_4^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$, inequality (3.33) is then a direct consequence of Lemma 3.5 applied to the function $f$.

We now consider the function $f$ defined by

$$f(i, x) = \sum_{j \in E} a(i, j, x) \left( \int \varphi(j, y) \mu(i, j, x)(dy) \right) \quad \forall (i, x) \in E \times \mathbb{R}^d.$$
Note that, although the function $\varphi$ has a compact support, we cannot in general expect the same property for the function $f$ under the Hypothesis (H-4), and we can only deduce that $f \in C_b(\mathbb{R}^d)^E$. We then have

$$
\begin{aligned}
& \left| m_K \bar{f}^{(i,K)} - \sum_{j \in E} \sum_{L \in \mathcal{M}} m_K a_{LK}^{(ji)} \bar{\varphi}^{(j,L)} \right| \\
& = \left| \int_{\{x \in K\}} \sum_{j \in E} a(i, j, x) \left( \sum_{L \in \mathcal{M}} \int_{\{y \in L\}} \varphi(j, y) \mu(i, j, x)(dy) \right) dx \right. \\
& \qquad \left. - \sum_{j \in E} \sum_{L \in \mathcal{M}} \int_{\{x \in K\}} a(i, j, x) \int_{\{y \in L\}} \mu(i, j, x)(dy) dx \frac{1}{m_L} \int_{\{z \in L\}} \varphi(j, z) dz \right| \\
& \le \sum_{j \in E} \sum_{L \in \mathcal{M}} \int_{\{x \in K\}} \int_{\{y \in L\}} a(i, j, x) \left| \varphi(j, y) - \frac{1}{m_L} \int_{\{z \in L\}} \varphi(j, z) dz \right| \mu(i, j, x)(dy) dx \\
& \le C_{\varphi} |\mathcal{M}| B m_K.
\end{aligned}
$$

We thus get that

$$|T_5^{\mathcal{M}, \delta t}(\varphi, t_1, t_2) - T_1^{\mathcal{M}, \delta t}(f, t_1, t_2)| \le C_{\varphi} |\mathcal{M}| BT.$$

Noting again that $T_2^{\mathcal{M}, \delta t}(f, t_1, t_2) = T_6^{\mathcal{M}, \delta t}(\varphi, t_1, t_2)$, the conclusion follows from the application of Lemma 3.5 to $f$. $\square$

We now study the difference between $T_7^{\mathcal{M}, \delta t}(\varphi, t_1, t_2)$ and $T_8^{\mathcal{M}, \delta t}(\varphi, t_1, t_2)$ defined, respectively, by (3.28) and (3.29).

LEMMA 3.7 Under assumptions (H), let $\varphi \in C_c^1(\mathbb{R}^d)^E$, $T > 0$ and $t_1, t_2 \in \mathbb{R}_+$ be such that $0 \le t_1 \le t_2 \le T$. Then, for each $\varepsilon_0 > 0$, there exist $\delta t_0$ and $\theta_0$ which only depends on $\varepsilon_0$, $T$, $\varphi$, $v$, $V_1$, $V_2$, $A$, $\text{Card}E$, $\rho_0$, and $\mu$ such that for all $\delta t \in [0, \delta t_0]$, for all admissible meshes $\mathcal{M}$ in the sense given in Section 2, such that $|\mathcal{M}| \le \theta_0 \delta t$, with $P^{\mathcal{M}, \delta t}$ given by (2.4)–(2.11) and $T_7^{\mathcal{M}, \delta t}(\varphi, t_1, t_2)$ and $T_8^{\mathcal{M}, \delta t}(f, t_1, t_2)$ defined by (3.28)–(3.29), we have that

$$|T_7^{\mathcal{M}, \delta t}(\varphi, t_1, t_2) - T_8^{\mathcal{M}, \delta t}(\varphi, t_1, t_2)| \le \varepsilon_0, \quad (3.35)$$

which means precisely that $\lim_{\delta t \to 0, |\mathcal{M}|/\delta t \to 0} (T_7^{\mathcal{M}, \delta t}(f, t_1, t_2) - T_8^{\mathcal{M}, \delta t}(f, t_1, t_2)) = 0$.

Proof. Let $\varepsilon > 0$ be given. Let $\mathcal{M}$ be an admissible mesh in the sense given in Section 2, and let $\delta t > 0$ be given. Using the notations from Section 2, we define the function $\bar{\varphi}_{\mathcal{M}}$ by $\bar{\varphi}_{\mathcal{M}}(i, x) = \bar{\varphi}^{(i,K)} = \frac{1}{m_K} \int_K \varphi(i, y) dy$ for all $x \in K$ and all $K \in \mathcal{M}$. We also note that, setting $C_{\varphi} = \sup_{i \in E} \|\nabla \varphi(i, \cdot)\|_{L^\infty(\mathbb{R}^d)}$, we have

$$|\bar{\varphi}_{\mathcal{M}}(i, x) - \varphi(i, x)| \le C_{\varphi} |\mathcal{M}| \quad \forall x \in \mathbb{R}^d, \forall i \in E. \quad (3.36)$$
We get from (2.6), (2.8), (2.9), (2.11) and (3.28)

$$
\begin{align*}
T_7^{\mathcal{M},\delta t}(\varphi, t_1, t_2) & = \sum_{i \in E} \sum_{K \in \mathcal{M}} \sum_{n=n_1}^{n_2-1} \int_K [P_{n+1}(i, x)\varphi(i, x) - P_n(i, g(i, x, -\delta t))J(i, x, -\delta t)\bar{\varphi}_{\mathcal{M}}(i, x)]dx \, ds \\
& = \frac{1}{\delta t} \sum_{i \in E} \sum_{n=n_1}^{n_2-1} \int_{n\delta t}^{(n+1)\delta t} \int_{\mathbb{R}^d} [P_{s+\delta t}^{\mathcal{M},\delta t}(i, x)\varphi(i, x) \\
& \phantom{=} - P_s^{\mathcal{M},\delta t}(i, g(i, x, -\delta t))J(i, x, -\delta t)\bar{\varphi}_{\mathcal{M}}(i, x)]dx \, ds \\
& = \frac{1}{\delta t} \sum_{i \in E} \int_{\mathbb{R}^d} \left( \int_{(n_1+1)\delta t}^{(n_2+1)\delta t} P_t^{\mathcal{M},\delta t}(i, y)\varphi(i, y)dt - \int_{n_1\delta t}^{n_2\delta t} P_s^{\mathcal{M},\delta t}(i, y)\bar{\varphi}_{\mathcal{M}}(i, g(i, y, \delta t))ds \right) dy.
\end{align*}
$$

This leads to $T_7^{\mathcal{M},\delta t}(\varphi, t_1, t_2) = -T_{14}^{\mathcal{M},\delta t}(\varphi, t_1, t_2) - T_{15}^{\mathcal{M},\delta t}(\varphi, t_1) + T_{16}^{\mathcal{M},\delta t}(\varphi, t_2)$, where we define

$$
T_{14}^{\mathcal{M},\delta t}(\varphi, t_1, t_2) = \frac{1}{\delta t} \sum_{i \in E} \int_{(n_1+1)\delta t}^{n_2\delta t} \int_{\mathbb{R}^d} P_s^{\mathcal{M},\delta t}(i, x)[\bar{\varphi}_{\mathcal{M}}(i, g(i, x, \delta t)) - \varphi(i, x)]dx \, ds,
$$

$$
T_{15}^{\mathcal{M},\delta t}(\varphi, t_1) = \frac{1}{\delta t} \sum_{i \in E} \int_{n_1 \delta t}^{(n_1+1) \delta t} \int_{\mathbb{R}^d} P_s^{\mathcal{M},\delta t}(i, y) \bar{\varphi}_{\mathcal{M}}(i, g(i, y, \delta t)) dy \, ds
$$

and

$$
T_{16}^{\mathcal{M},\delta t}(\varphi, t_2) = \frac{1}{\delta t} \sum_{i \in E} \int_{n_2 \delta t}^{(n_2+1) \delta t} \int_{\mathbb{R}^d} P_s^{\mathcal{M},\delta t}(i, y) \varphi(i, y) dy \, ds.
$$

We first remark that

$$
T_{15}^{\mathcal{M},\delta t}(\varphi, t_1) = \sum_{i \in E} \int_{\mathbb{R}^d} P_{t_1}^{\mathcal{M},\delta t}(i, y) \bar{\varphi}_{\mathcal{M}}(i, g(i, y, \delta t)) dy,
$$

and that

$$
T_{16}^{\mathcal{M},\delta t}(\varphi, t_2) = \sum_{i \in E} \int_{\mathbb{R}^d} P_{t_2}^{\mathcal{M},\delta t}(i, y) \varphi(i, y) dy.
$$

Hence, $T_{16}^{\mathcal{M},\delta t}(\varphi, t_2)$ is exactly one of the three terms constituting $T_8^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$.

Applying Lemma 3.1 and using the fact that the support of $\varphi$ is compact gives the existence of some
$C'_1 > 0$ only depending on $V_1, V_2, T$ and $\varphi$ such that $|\bar{\varphi}_{\mathcal{M}}(i, g(i, y, \delta t)) - \bar{\varphi}_{\mathcal{M}}(i, y)| \le (\delta t + 2|\mathcal{M}|)C'_1$.
Let $\theta_1$ be given (for instance $\theta_1 = 1$), we then get that there exists $\delta t_1 > 0$ such that, for $\delta t \le \delta t_1$ and
$|M| \le \theta_1 \delta t$, we have

$$ \left| T_{15}^{M, \delta t}(\varphi, t_1) - \sum_{i \in E} \int_{\mathbb{R}^d} \varphi(i, x) P_{t_1}^{M, \delta t}(i, x) dx \right| \le \varepsilon. $$

We have $T_{14}^{M, \delta t}(\varphi, t_1, t_2) = T_{17}^{M, \delta t}(\varphi, t_1, t_2) + T_{18}^{M, \delta t}(\varphi, t_1, t_2)$, where we define

$$ T_{17}^{M, \delta t}(\varphi, t_1, t_2) = \frac{1}{\delta t} \sum_{i \in E} \int_{(n_1+1)\delta t}^{n_2\delta t} \int_{\mathbb{R}^d} P_s^{M, \delta t}(i, x)[\varphi(i, g(i, x, \delta t)) - \varphi(i, x)]dx ds $$

and

$$ T_{18}^{M, \delta t}(\varphi, t_1, t_2) = \frac{1}{\delta t} \sum_{i \in E} \int_{(n_1+1)\delta t}^{n_2\delta t} \int_{\mathbb{R}^d} P_s^{M, \delta t}(i, x)[\bar{\varphi}_M(i, g(i, x, \delta t)) - \varphi(i, g(i, x, \delta t))]dx ds. $$

Defining $T_{19}^{M, \delta t}(\varphi, t'_1, t'_2)$ by

$$ T_{19}^{M, \delta t}(\varphi, t'_1, t'_2) = \sum_{i \in E} \int_{t'_1}^{t'_2} \int_{\mathbb{R}^d} P_s^{M, \delta t}(i, x) \mathbf{v}(i, x) \cdot \nabla \varphi(i, x) dx ds $$

for all $t'_1, t'_2 \in \mathbb{R}^+$ such that $0 \le t'_1 \le t'_2 \le T$, we have $T_{19}^{M, \delta t}(\varphi, t_1, t_2) = T_{19}^{M, \delta t}(\varphi, t_1, (n_1+1)\delta t) + T_{19}^{M, \delta t}(\varphi, (n_1+1)\delta t, n_2\delta t) + T_{19}^{M, \delta t}(\varphi, n_2\delta t, t_2)$. We note that there exists $\delta t_2 > 0$, only depending on $\varphi, V_1, V_2, T$ and a bound of $\mathbf{v}$ on supp($\varphi$), such that, for $\delta t \le \delta t_2$, $|T_{19}^{M, \delta t}(\varphi, t_1, (n_1+1)\delta t)| \le \varepsilon$ and $|T_{19}^{M, \delta t}(\varphi, n_2\delta t, t_2)| \le \varepsilon$. We can write, since $\partial_\tau \varphi(i, g(i, x, \tau)) = \mathbf{v}(i, g(i, x, \tau)) \cdot \nabla \varphi(i, g(i, x, \tau))$,

$$ \varphi(i, g(i, x, \delta t)) - \varphi(i, x) = \int_{0}^{\delta t} \mathbf{v}(i, g(i, x, \tau)) \cdot \nabla \varphi(i, g(i, x, \tau)) d\tau. $$

Thanks to the regularity properties of $\varphi$, this gives the existence of some function $h: \mathbb{R}_+ \to \mathbb{R}$ with $h(x) \to 0$ as $x \to 0$, which only depends on $\varphi, V_1, V_2, T$ and on the properties of $\mathbf{v}$ on supp($\varphi$), such that

$$ \left| \frac{\varphi(i, g(i, x, \delta t)) - \varphi(i, x)}{\delta t} - \mathbf{v}(i, x) \cdot \nabla \varphi(i, x) \right| \le h(\delta t) \quad \forall x \in \mathbb{R}^d. $$

We thus obtain that there exists $\delta t_3 > 0$ such that, for $\delta t \le \delta t_3$, we get

$$ |T_{17}^{M, \delta t}(\varphi, t_1, t_2) - T_{19}^{M, \delta t}(\varphi, (n_1+1)\delta t, n_2\delta t)| \le \varepsilon. $$

We now observe that, thanks to (3.36),

$$ |T_{18}^{M, \delta t}(\varphi, t_1, t_2)| \le C_{\varphi} T \frac{|M|}{\delta t}. $$

We can then choose $\theta_2 = \varepsilon$, which gives, for $|M| \le \delta t \theta_2$, $|T_{18}^{M, \delta t}(\varphi, t_1, t_2)| \le C_{\varphi} T \varepsilon$. Gathering the above results, it suffices to choose $\varepsilon$, such that $4\varepsilon + C_{\varphi} T \varepsilon = \varepsilon_0$, $\delta t_0 = \min(\delta t_1, \delta t_2, \delta t_3)$ and $\theta_0 = \min(\theta_1, \theta_2)$, to conclude the proof of the lemma. $\square$
**REMARK 3.2** The hypothesis $|\mathcal{M}|/\delta t \to 0$ is sufficient to guarantee the existence of some $\theta > 0$ such that $|\mathcal{M}| \le \theta\delta t$. This hypothesis is only imposed in order to ensure that $T_{18}^{\mathcal{M},\delta t}(\varphi, t_1, t_2) \to 0$, and was not needed elsewhere. In some special cases, however, the term $T_{18}^{\mathcal{M},\delta t}(\varphi, t_1, t_2)$ vanishes so that the condition $|\mathcal{M}|/\delta t \to 0$ is not required any more. Consider, for example, the case $d = 1$, $\mathcal{M} = \{I_k = [k\delta t, (k+1)\delta t] \text{ with } k \in \mathbb{Z}\}$ and $g(i, x, s) = x + m_i s$ with $m_i \in \mathbb{Z}$.

In the case where $\rho_0(i, x) = u_0(i, x)dx$, with $u_0 \in L^2(\mathbb{R}^d)^E$ and $\mu(i, j, x)(dy) = \delta_x$, it is possible to prove a weak bounded variation inequality, namely, the existence of some $C > 0$ such that

$$ \sum_{i \in E} \sum_{n=0}^{N} \sum_{K, L \in \mathcal{M}} m_{KL}^{(i)} (p_n^{(i,L)} - p_{n+1}^{(i,K)})^2 \le C, \quad N\delta t \le T < (N+1)\delta t, $$

and an $L^2(\mathbb{R}^d \times \mathbb{R}_+)$ estimate

$$ \sum_{i \in E} \sum_{K \in \mathcal{M}} m_K (p_{n+1}^{(i,K)})^2 \le C \quad \forall n \in \mathbb{N} \text{ s.t. } n\delta t \le T. $$

The weak limit is then a function of $L^2(\mathbb{R}^d \times \mathbb{R}_+)^E$, and the convergence can be proved under the hypothesis $|\mathcal{M}| \le \theta\delta t$. We do not give the details of this result here since it does not lead to a convergence proof for more general $\rho_0$ and $\mu$.

### 3.3 Continuity and compactness

We can now state a lemma which ensures the continuity of the approximate solution with respect to the time variable. Let $\varphi \in C_b(\mathbb{R}^d)^E$ be a given function. For a given admissible mesh in the sense given in Section 2, and for a given $\delta t > 0$, using all the notations of Section 2, we set

$$ S_{\varphi}(\mathcal{M}, \delta t, s, t) = \sum_{i \in E} \int_{\mathbb{R}^d} \varphi(i, x) P_t^{\mathcal{M}, \delta t}(i, x)dx \\ - \sum_{i \in E} \int_{\mathbb{R}^d} \varphi(i, x) P_s^{\mathcal{M}, \delta t}(i, x)dx \quad \forall s, t \in \mathbb{R}_+. \quad (3.37) $$

The statement is then the following.

**LEMMA 3.8** Under assumptions (H), let $\varphi \in C_c^1(\mathbb{R}^d)^E$ and let $T > 0$ be given. Then, there exists $C_S$, which only depends on $T$, $\varphi$, $V_1$, $V_2$, $A$, Card$E$, $\rho_0$ and $\mu$, such that, for each $\varepsilon_0 > 0$, there exist $\delta t_0$ and $\theta_0$ which only depends on $\varepsilon_0$, $T$, $\varphi$, $v$, $V_1$, $V_2$, $A$, Card$E$, $\rho_0$, and $\mu$ such that for all $\delta t \in [0, \delta t_0]$, for all admissible meshes $\mathcal{M}$ in the sense given in Section 2, such that $|\mathcal{M}| \le \theta_0\delta t$, with $P^{\mathcal{M},\delta t}$ given by (2.4)–(2.11) and $S_{\varphi}(\mathcal{M}, \delta t, s, t)$ defined by (3.37), we have that

$$ |S_{\varphi}(\mathcal{M}, \delta t, s, t)| \le C_S |t-s| + \varepsilon_0, \quad (3.38) $$

which means precisely that $\limsup_{\delta t \to 0, |\mathcal{M}|/\delta t \to 0} |S_{\varphi}(\mathcal{M}, \delta t, s, t)| \le C_S |t-s|$.

**REMARK 3.3** We could complete the proof of the above lemma assuming only $|\mathcal{M}| \le \theta\delta t$ instead of $|\mathcal{M}|/\delta t \to 0$, which would only lead to a different choice for $C_S$.

*Proof.* Let $\varepsilon_0 > 0$ be given. Let $\theta_0$ be given by Lemma 3.7, and let $\delta t_0$ be the minimum between the values $\delta t_0$ given by Lemmas 3.5, 3.6 and 3.7, setting in $\theta = \theta_0$ in Lemmas 3.5 and 3.6 let $\delta t \in (0, \delta t_0)$ and $\mathcal{M}$ be an admissible mesh in the sense given in Section 2, such that $|\mathcal{M}| \le \theta_0\delta t$, with $P^{\mathcal{M},\delta t}$ given by (2.4)–(2.11).
Starting from (3.31), we get

$$
\begin{aligned}
& \left| \sum_{i \in E} \int_{\mathbb{R}^d} \varphi(i, x) P_{t_1}^{\mathcal{M}, \delta t}(i, x) dx - \sum_{i \in E} \int_{\mathbb{R}^d} \varphi(i, x) P_{t_2}^{\mathcal{M}, \delta t}(i, x) dx \right| \\
& \leq \left| -T_4^{\mathcal{M}, \delta t}(\varphi, t_1, t_2) + T_6^{\mathcal{M}, \delta t}(\varphi, t_1, t_2) - T_8^{\mathcal{M}, \delta t}(\varphi, t_1, t_2) \right| \\
& \qquad + \left| \sum_{i \in E} \int_{t_1}^{t_2} \int_{\mathbb{R}^d} \left( \sum_{j \in E} a(i, j, x) \int \varphi(j, y) \mu(i, j, x) (dy) \right) \cdot \right. \\
& \qquad \left. - b(i, x) \varphi(i, x) + \mathbf{v}(i, x) \cdot \nabla \varphi(i, x) \right| P_t^{\mathcal{M}, \delta t}(i, x) dx \, dt
\end{aligned}
$$

Using (3.30) and then Lemmas 3.6 and 3.7, we have

$$
\begin{aligned}
& | - T_4^{\mathcal{M}, \delta t}(\varphi, t_1, t_2) + T_6^{\mathcal{M}, \delta t}(\varphi, t_1, t_2) - T_8^{\mathcal{M}, \delta t}(\varphi, t_1, t_2)| \\
&= |T_3^{\mathcal{M}, \delta t}(\varphi, t_1, t_2) - T_4^{\mathcal{M}, \delta t}(\varphi, t_1, t_2) + T_6^{\mathcal{M}, \delta t}(\varphi, t_1, t_2) \\
& \qquad - T_5^{\mathcal{M}, \delta t}(\varphi, t_1, t_2) - T_8^{\mathcal{M}, \delta t}(\varphi, t_1, t_2) + T_7^{\mathcal{M}, \delta t}(\varphi, t_1, t_2)| \\
&\le 3\epsilon_0.
\end{aligned}
$$

Let $V$ be a bound of $\mathbf{v}$ on $\text{supp}(\varphi)$. Setting

$$
C_S = 2B \sup_{i \in E} \|\varphi(i, \cdot)\|_{L^{\infty}(\mathbb{R}^d)} + V \sup_{i \in E} \|\nabla\varphi(i, \cdot)\|_{L^{\infty}(\mathbb{R}^d)}
$$

provides the result by substituting $\epsilon_0/3$ by $\epsilon_0$ at the beginning of the proof.

We can now conclude with the convergence theorem.

**THEOREM 3.1** Under assumptions (H), let ($M_m, \delta t_m)_{m \in \mathbb{N}$ be a sequence such that, for all $m \in \mathbb{N}$, $M_m$ is an admissible mesh in the sense given in Section 2 and such that $\delta t_m \to 0$ and $|M_m| / \delta t_m \to 0$ as $m \to \infty$. Then, for all $t \in \mathbb{R}_+$, the tight sequence of probability measures ($P_t^{M_m, \delta t_m}(\cdot, x) dx)_{m \in \mathbb{N}$ weakly converges to a probability measure on $E \times \mathbb{R}^d$, denoted by $\rho(t)(\cdot, dx)$. Moreover, this family $(\rho(t)(\cdot, dx))_{t \in \mathbb{R}_+}$ is the unique weak solution of the Chapman-Kolmogorov equations (1.2).

*Proof.* Let $(t_p)_{p \in \mathbb{N}}$ be a sequence of real numbers, dense in $\mathbb{R}_+$ and let $(M_{\sigma_p(k)}, \delta t_{\sigma_p(k)})_{k \in \mathbb{N}}$ be a sub-sequence of $(M_m, \delta t_m)_{m \in \mathbb{N}}$ such that $(P_{t_p}^{M_{\sigma_p(k)}, \delta t_{\sigma_p(k)}}(\cdot, x)dx)_{k \in \mathbb{N}}$ is weakly convergent to a probability measure (thanks to the tightness property Lemma 3.4 and to the Prohorov theorem). Using a diagonal method, we can choose a sub-sequence $(M_\sigma(k), \delta t_\sigma(k))_{k \in \mathbb{N}}$ of $(M_m, \delta t_m)_{m \in \mathbb{N}}$, such that $(P_{t_p}^{M_\sigma(k), \delta t_\sigma(k)}(\cdot, x)dx)_{k \in \mathbb{N}}$ is weakly convergent for all $p \in \mathbb{N}$. Such a sub-sequence is again denoted by $(M_m, \delta t_m)_{m \in \mathbb{N}}$ in the following and $(P_{t_p}^{M_m, \delta t_m}(\cdot, x)dx)_{m \in \mathbb{N}}$ is now weakly convergent to a probability distribution $\rho(t_p)(\cdot, dx)$ for all $p \in \mathbb{N}$.

Let us now prove that, for all $\psi \in C_c^1(\mathbb{R}^d)^E$, the sequence $(\sum_{i \in E} \int_{\mathbb{R}^d} \psi(i, x) P_t^{M_m, \delta t_m}(i, x) dx)_{m \in \mathbb{N}}$ is a Cauchy sequence for all $t > 0$.
Let $k, m \in \mathbb{N}$. For all $p \in \mathbb{N}$, we may write, using the notation (3.7),

$$ \left| \sum_{i \in E} \int_{\mathbb{R}^d} \psi(i, x) P_t^{M_k, \delta t_k}(i, x) dx - \sum_{i \in E} \int_{\mathbb{R}^d} \psi(i, x) P_t^{M_m, \delta t_m}(i, x) dx \right| \\ \le |S_\psi(M_k, \delta t_k, t, t_p)| + |S_\psi(M_m, \delta t_m, t, t_p)| \\ + \left| \sum_{i \in E} \int_{\mathbb{R}^d} \psi(i, x) P_{t_p}^{M_k, \delta t_k}(i, x) dx - \sum_{i \in E} \int_{\mathbb{R}^d} \psi(i, x) P_{t_p}^{M_m, \delta t_m}(i, x) dx \right| \quad (3.39) $$

Let $T = t + 1$ and let $\varepsilon \in (0, 1)$ be given. We first choose $p \in \mathbb{N}$ such that $|t - t_p| \le \varepsilon$. Since the sequence $(P_{t_p}^{M_m, \delta t_m}(\cdot, x)dx)_{m \in \mathbb{N}}$ is weakly convergent, there exists $m_0$ such that, for all $m, k \ge m_0$,

$$ \left| \sum_{i \in E} \int_{\mathbb{R}^d} \psi(i, x) P_{t_p}^{M_k, \delta t_k}(i, x) dx - \sum_{i \in E} \int_{\mathbb{R}^d} \psi(i, x) P_{t_p}^{M_m, \delta t_m}(i, x) dx \right| \le \varepsilon. $$

We can now choose $m_1 \ge 0$ such that, for all $m \ge m_1$, $\delta t_m \le \delta t_0$ and $|M_m| \le \theta_0 \delta t_m$, where $\delta t_0$ and $\theta_0$ are given by Lemma 3.8 for the value $\varepsilon$. We thus have, for $m, k \ge m_1$,

$$ |S_\psi(M_k, \delta t_k, t, t_p)| + |S_\psi(M_m, \delta t_m, t, t_p)| \le 2(C_S|t - t_p| + \varepsilon) \le (2C_S + 2)\varepsilon. $$

This completes the proof that $(\sum_{i \in E} \int_{\mathbb{R}^d} \psi(i, x) P_t^{M_m, \delta t_m}(i, x) dx)_{m \in \mathbb{N}}$ is a Cauchy sequence for all $t > 0$. We easily obtain the same property for all $\psi \in C_c(\mathbb{R}^d)^E$ by density of $C_c^1(\mathbb{R}^d)^E$ and then for all $\psi \in C_b(\mathbb{R}^d)^E$, thanks to the tightness Lemma 3.5. Therefore, there exists a probability measure $\rho(t)(\cdot, dx)$ such that the sequence $(P_t^{M_m, \delta t_m}(\cdot, x)dx)_{m \in \mathbb{N}}$ is weakly convergent to $\rho(t)(\cdot, dx)$.

Moreover, passing to the limit in (3.38), we get that, for all $\varphi \in C_c^1(\mathbb{R}^d)^E$,

$$ \left| \sum_{i \in E} \int_{\mathbb{R}^d} \varphi(i, x) \rho(t)(i, dx) - \sum_{i \in E} \int_{\mathbb{R}^d} \varphi(i, x) \rho(s)(i, dx) \right| \le C_S |s - t| \quad \forall s, t \in \mathbb{R}_+, $$

which gives the continuity of the mapping $t \mapsto \int \sum_{i \in E} \varphi(i, x) \rho(t)(i, dx)$. The continuity property is easily extended by density to all $\varphi \in C_c(\mathbb{R}^d)^E$ and to all $\varphi \in C_b(\mathbb{R}^d)^E$ due to tightness. Thanks to Lemmas 3.6 and 3.7 and the dominated convergence theorem, we have the following properties for all $\varphi \in C_c^1(\mathbb{R}^d)^E$ and $t_1, t_2 \ge 0$:

$$ \lim_{m \to \infty} T_3^{M_m, \delta t_m}(\varphi, t_1, t_2) = \lim_{m \to \infty} T_4^{M_m, \delta t_m}(\varphi, t_1, t_2) = \sum_{i \in E} \int_{t_1}^{t_2} \int_{\mathbb{R}^d} b(i, x) \varphi(i, x) \rho(t)(i, dx) dt, $$

$$ \lim_{m \to \infty} T_5^{M_m, \delta t_m}(\varphi, t_1, t_2) = \lim_{m \to \infty} T_6^{M_m, \delta t_m}(\varphi, t_1, t_2) \\ = \sum_{i \in E} \int_{t_1}^{t_2} \int_{\mathbb{R}^d} \sum_{j \in E} a(i, j, x) \left( \int \varphi(j, y) \mu(i, j, x)(dy) \right) \rho(t)(i, dx) dt, $$

$$ \lim_{m \to \infty} T_7^{M_m, \delta t_m}(\varphi, t_1, t_2) = \lim_{m \to \infty} T_8^{M_m, \delta t_m}(\varphi, t_1, t_2) \\ = - \sum_{i \in E} \left[ \int_{t_1}^{t_2} \int_{\mathbb{R}^d} \mathbf{v}(i, x) \cdot \nabla \varphi(i, x) \rho(t)(i, dx) dt \right. \\ \left. + \int_{\mathbb{R}^d} \varphi(i, x) \rho(t_1)(i, dx) - \int_{\mathbb{R}^d} \varphi(i, x) \rho(t_2)(i, dx) \right]. $$
We now pass to the limit in (3.30). Since it is clear that $\rho(0)(i, dx) = \rho_0(i, dx)$ thanks to (2.4), we deduce that $(\rho(t)(i, dx))_{t \in \mathbb{R}_+}$ is the solution to the Chapman–Kolmogorov equations (1.2). Thanks to the uniqueness result proven in Cocozza-Thivent et al., we deduce in a classical way that the whole initial sequence converges in the same sense. $\square$

**REMARK 3.4 (ALTERNATE ASSUMPTIONS)** In assumptions (H), transition rates are assumed to be bounded, which is not always the case in practice (consider, for example, Weibull distribution with shape parameter greater than one). In case of unbounded continuous transition rates, the results are still valid, however, under additional assumptions on $\mu(i, j, x)(dy)$ and on the initial distribution $\rho_0(\cdot, dx)$ which allow one to control the expansion of the support of $\rho(t)(\cdot, dx)$ (or alternatively of $P_t^{M_n, \delta_{tM}}(i, x))$. Indeed, assume for example that the support of $\rho_0(\cdot, dx)$ is included in some $\bar{B}(0, M_0)(\text{the closure of } B(0, M_0))$ and that for all $i, j \in E$, the support of $\mu(i, j, x)(dy)$ is included in $\bar{B}(0, |x|)$, it is then easy to show that, for all $T > 0$, all $0 < t < T$, the support of $\rho(t)(\cdot, dx)$ is bounded and is included in some $\bar{B}(0, M)$ where $M$ depends on $M_0, V_1, V_2$ and $T$. As transition rates are bounded on $\bar{B}(0, M)$, we first derive from Cocozza-Thivent et al. the uniqueness for the solution of (1.2) and second, similarly as in the present paper, the weak convergence of $(P_t^{M_m, \delta_{tM}}(\cdot, x)dx)_{m \in \mathbb{N}}$ to $\rho(t)(\cdot, dx)$ for all $0 < t < T$.

## 4. Numerical example

We present the results obtained using the finite-volume scheme (2.4)–(2.10) on a benchmark taken from Labeau (1996). The states of the system are defined by $E = \{1, 2, 3\}$, and the environment variables are described by $x \in \mathbb{R}^2$. The only possible transitions between states are $1 \to 2 \to 3$, which means that $a(i, j, x) \neq 0$ only for $(i, j) = (1, 2)$ or $(i, j) = (2, 3)$. We assume that $a(1, 2, x) = \lambda_1$ and $a(2, 3, x) = \lambda_2$ for all $x = (x^{(1)}, x^{(2)}) \in \mathbb{R}^2$, where $\lambda_1 = 0.5$ and $\lambda_2 = 0.3$ and that $\mu(i, j, x)(dy) = \delta_x(y)$ (this means that the environmental variables are not modified by transitions). The velocity field is assumed to satisfy $v(i, x) = A_i x$, with $A_1 = \begin{pmatrix} 1.5 & 0 \\ 1 & 0.5 \end{pmatrix}$, $A_2 = \begin{pmatrix} 1 & 0 \\ 0.75 & 0.25 \end{pmatrix}$ and $A_3 = \begin{pmatrix} 0.75 & 0 \\ 0 & 0.5 \end{pmatrix}$. The initial data $\rho_0$ is given by the Dirac measure $\delta_{x_0}$ with $x_0 = (1, 1)$ in the state $i = 1$, which means that $\sum_{i \in E} \int_{\mathbb{R}^2} \varphi(i, x) \rho_0(i, dx) = \varphi(1, x_0)$, for all $\varphi \in C_b^1(\mathbb{R}^2, \mathbb{R})^E$.

An analytical solution can be obtained in this case Labeau (1996), which satisfies that $\rho(t)(1, \cdot)$ is a weighted Dirac measure located at the point $g(1, x_0, t)$ of $\mathbb{R}^2$, that $\rho(t)(2, \cdot)$ is supported by a line of $\mathbb{R}^2$ and that $\rho(t)(3, \cdot)$ admits a density with respect to the Lebesgue measure, denoted $p(t)$, supported by a triangle.

We consider meshes on $\mathbb{R}^2$ consisting of squares with side $|M|$. In order to confine the computations as much as possible to the geometrical domain which supports $\rho(t)$, we performed consecutive simulations with the values $|M| = 5 \times 10^{-3}$, $|M| = 3.13 \times 10^{-4}$ and $|M| = 10^{-5}$, each simulation being used to restrict the computational domain of the following one. Figure 1 presents the consecutive numerical supports of the computational domain. In this figure, the triangular line represents the analytical contour of the support of $\rho(0.1)(3, \cdot)$. We then compare, in Fig. 2, the numerical solution (top) with the analytical one (bottom) for the density $p(0.1)$ of $\rho(0.1)(3, \cdot)$, with $|M| = 10^{-5}$. We observe good agreement between the approximate and the analytical solution, with smoothing due to numerical diffusion. Nevertheless, this solution appears to be sufficiently accurate from the point of view of the densities $p_x(t)$ and $p_y(t)$ of the marginal distributions in the state $i = 3$, defined by $p_x(t)(x^{(1)}) = \int_{\mathbb{R}} p(t)(x^{(1)}, x^{(2)})dx^{(2)}$, for a.e. $x^{(1)} \in \mathbb{R}$, and $p_y(t)(x^{(2)}) = \int_{\mathbb{R}} p(t)(x^{(1)}, x^{(2)})dx^{(1)}$, for a.e. $x^{(2)} \in \mathbb{R}$. We compare in Figs 3 and 4, respectively, the analytical and numerical results for $p_x(t)$ and $p_y(t)$: these figures show an excellent agreement.
<figure><img src="image_2.png" /><figcaption>FIG. 1. Consecutive numerical supports of the computational domain.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Density $p(0.1)$. Top: numerical approximation with $|M| = 10^{-5}$, bottom: analytical solution.</figcaption></figure>

If we compare the use of the finite-volume method for the approximation of the marginal distributions, and Monte Carlo simulations, we can make the following observations:

* For both methods, it is possible to adjust the precision, by increasing the number of histories for Monte Carlo simulations and by refining the mesh for the finite-volume method. In the above example (see Labeau (1996)), the computing time is here about three quarters of an hour on a PC, which is longer than those given in Labeau (1996). The precision in Labeau (1996) is, however, lower (this is shown by a comparison between the marginal distributions).

* The finite-volume method has been shown to be efficient in some cases where some other methods have been unsuccessful (see Cocozza-Thivent & Eymard and Cocozza-Thivent & Eymard (2003)), with relatively short computing times on 1D problems. However, one can wonder if these computing times will remain acceptable if the dimension of the problem increases. On the contrary, it seems that Monte Carlo simulation is the method of choice for higher dimensional problems.
<figure><img src="image_4.png" /><figcaption>FIG. 3. Analytical and numerical function $p_x(0.1)(x)$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4. Analytical and numerical function $p_y(0.1)(y)$.</figcaption></figure>

* Finally, an advantage of the finite-volume methods is that it can be extended to the computation of some stationary cases, with a slight modification of the scheme. In this case, a linear system must then be solved. One can expect the development of some finite-volume schemes which are able to handle long-term reliability studies, by mixing a transient model with the techniques used for the stationary problem.

## REFERENCES

* ABIA, L. M., ANGULO, O. & LÓPEZ-MARCOS, J. C. (2004) Size-structured population dynamics models and their numerical solutions. *Discrete Contin. Dyn. Syst. Ser. B*, **4**, 1203–1222.

* ACKLEH, A. S., BANKS, H. T. & DENG, K. (2002) A finite difference approximation for a coupled system of nonlinear size-structured populations. *Nonlinear Anal.*, **50**, 727–748.
ACKLEH, A. S. & ITO, K. (1998) An implicit finite difference scheme for the nonlinear size-structured population model. *Numer. Funct. Anal. Optim.*, **8**, 865–884.

BELL, G. I. & ANDERSON, E. C. (1967) Cell growth and division I. A mathematical model with applications to cell volume distributions in mammalian suspension cultures. *Biophys. J.*, **7**, 329–351.

BILLINGSLEY, P. (1968) *Convergence of Probability Measures*, Wiley Series in Probability and Mathematical Statistics. New York: John Wiley & sons.

COCOZZA-THIVENT, C. & EYMARD, R. (2003) Marginal distributions of a semi-Markov process and their computations. *Ninth ISSAT International Conference on Reliability and Quality in Design* (H. Pham & S. Yamada eds). New Brunswick, NJ: International Society of Science and Applied Technologies.

COCOZZA-THIVENT, C. & EYMARD, R. (2004) Approximation of the marginal distributions of a semi-Markov process using a finite volume scheme *ESAIM: M2AN*, **38**, 853–875.

COCOZZA-THIVENT, C. & EYMARD, R. (to appear) Numerical computation of the marginal distributions of a semi-Markov process (H. Pharm ed.). *Reliab. Model. Anal. Optim*. (to appear).

COCOZZA-THIVENT, C. EYMARD, R., MERCIER, S. & ROUSSIGNOL, M. (to appear) Characterization of the marginal distributions of Markov processes used in dynamic reliability. *J. Appl. Math. Stoch. Anal.* (to appear).

DAVIS, M. H. A. (1984) Piecewise-deterministic Markov processes: a general class of non-diffusion stochastic models. *J. R. Stat. Soc. Ser. B*, **46**, 353–388.

DAVIS, M. H. A. (1993) *Markov Models and Optimization*. London: Chapman & Hall.

EYMARD, R., GALLOUËT, T. & HERBIN, R. (2000) *Finite Volume Methods, Handbook of Numerical Analysis* (P. G. Ciarlet & J. L. Lions eds), VII. North-Holland, Amsterdam: Elsevier, pp. 713–1020.

FELLER, W. (1966) *An Introduction to Probability Theory and its Applications*, vol. II. New York: Wiley.

LABEAU, P. E. (1996) A Monte Carlo estimation of the marginal distributions in a problem of probabilistic dynamics. *Reliab. Eng. Sys. Saf.*, **52**, 65–75.

MISCHLER, S. PERTHAME, B. & RYZHIK, L. (2002) Stability in a nonlinear population maturation model. *Math. Models Methods Appl. Sci.*, **12**, 1–22.

SINKO, J. W. & STREIFER, W. (1967) A new model for age-size structure of a population. *Ecology*, **48**, 910–918.