# Approximation of the marginal distributions of a semi-Markov process using a finite volume scheme

Christiane Cocozza-Thivent, Robert Eymard

► To cite this version:

Christiane Cocozza-Thivent, Robert Eymard. Approximation of the marginal distributions of a semi-Markov process using a finite volume scheme. ESAIM: Mathematical Modelling and Numerical Analysis, 2004, 38 (5), pp.853-875. 10.1051/m2an:2004043. hal-00693113

HAL Id: hal-00693113

https://hal.science/hal-00693113

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
# APPROXIMATION OF THE MARGINAL DISTRIBUTIONS OF A SEMI-MARKOV PROCESS USING A FINITE VOLUME SCHEME

CHRISTIANE COCOZZA-THIVENT<fnref n="1" /> AND ROBERT EYMARD<fnref n="1" />

**Abstract.** In the reliability theory, the availability of a component, characterized by non constant failure and repair rates, is obtained, at a given time, thanks to the computation of the marginal distributions of a semi-Markov process. These measures are shown to satisfy classical transport equations, the approximation of which can be done thanks to a finite volume method. Within a uniqueness result for the continuous solution, the convergence of the numerical scheme is then proven in the weak measure sense, and some numerical applications, which show the efficiency and the accuracy of the method, are given.

**Mathematics Subject Classification.** 60K15, 60K20, 65C20, 65M60.

## INTRODUCTION

The reliability theory is devoted to the study of the probability that a system be available at a given time, when the failure or the repair of each component of this system is modeled using stochastic processes. The general mathematical framework of such problems is, in the general case, highly complex. In this paper, we focus on the more simple study of the availability of a single component, the state of which belongs to a finite set, the transitions of state being determined by probabilistic laws. Let us first consider the simple example of a component which has been created before the time $t = 0$, and which can either be working, either be in repair. Hence the set of the possible states is defined by $E = \{0, 1\}$, 1 for the working state, 0 for the repair state.

Let us denote the failure rate (resp. the repair rate) of the component by $b(1, \cdot)$ (resp. $b(0, \cdot)$). This means that the probability for the component to fail (resp. to be repaired) in the time interval $(t, t+dt)$, knowing that the elapsed time in the working (resp. failure) state is $t$, is equal to $b(1, t)dt$ (resp. $b(0, t)dt$). We assume that as soon as the component fails, the repair begins immediately and that, after being repaired, the component is equivalent to a new one and is immediately working.

Considering this component at a given time $t \ge 0$, let us assume that the probability, the component be in state 1 (resp. 0) and the elapsed time in this working (resp. failure) state be in the time interval $(x, x + dx)$, is given by the measure $\bar{u}_t(1, \cdot)dx$ (resp. $\bar{u}_t(0, \cdot)dx$). The functions $\bar{u}_t(i, \cdot) \in L^1(\mathbb{R}_+)$, for $(i, t) \in E \times \mathbb{R}_+$, are then (non strictly) positive, and must verify, for $t \in \mathbb{R}_+$, that $\sum_{i \in E} \int_{\mathbb{R}_+} \bar{u}_t(i, x)dx = 1$ (thanks to the law of

<footnote>Keywords and phrases. Renewal equation, semi-Markov process, convergence of a finite volume scheme.</footnote>

<footnote n="1">¹ Laboratoire d'Analyse et de Mathématiques Appliquées (UMR 8050 CNRS), Université de Marne la Vallée, Cité Descartes, 5 boulevard Descartes, Champs sur Marne, 77454 Marne-la-Vallée Cedex 2, France.</footnote>

<footnote>e-mail: cocozza@univ-mlv.fr, eymard@univ-mlv.fr</footnote>
large numbers, these functions could as well be interpreted as providing the distribution of the ages in each of the two states, 1 and 0, for a population of similar components, all created before the time $t = 0$).

It can then be shown (see Sect. 1 for a detailed proof) that all the above hypotheses give that the function $(t, i, x) \mapsto \bar{u}_t(i, x)$ is the weak solution of the following equations

$$ \frac{\partial \bar{u}_t}{\partial t}(i, x) + \frac{\partial \bar{u}_t}{\partial x}(i, x) = -b(i, x) u_t(i, x), \quad \text{for } (i, t, x) \in E \times \mathbb{R}_+ \times \mathbb{R}_+, \quad (1) $$

with an initial condition

$$ \bar{u}_0(i, x) = \bar{u}_{\text{ini}}(i, x), \quad \text{for } (i, x) \in E \times \mathbb{R}_+, \quad (2) $$

and a boundary condition

$$ \bar{u}_t(i, 0) = \int_{\mathbb{R}_+} b(1 - i, x) \bar{u}_t(1 - i, x) dx, \quad \text{for } (i, t) \in E \times \mathbb{R}_+. \quad (3) $$

Equations (1)–(3) constitute a system of linear hyperbolic equations on $\mathbb{R}_+ \times \mathbb{R}_+$, only coupled here by the boundary conditions. Such a model belongs to the class of the Markov renewal models which are generalizations of the renewal model (see for example [2] or [3]). When the state space $E$ has only one element, we obtain the class of renewal processes which have been extensively studied from the mathematical point of view (see for example [10] or [12] and references therein). Note that similar models are widely used in the framework of population dynamics (the so-called McKendrick-Von Foerster model is a one-state version of (1)–(3), see for example [14] and references therein). In the particular case where $E$ has two elements (as in the example considered in this introduction), one says that this model is an alternative renewal model. We consider in this paper an extension of Model (1)–(3), used for computing the marginal distributions of the stochastic process constituted by the pair $(\eta_t, X_t)$, where $\eta_t \in E$ denotes the state at time $t \in \mathbb{R}_+$ of a given component (we now consider the case of any finite set $E$ with, in practice, at least two elements) and $X_t \in \mathbb{R}_+$ is the elapsed time since this component is in the state $\eta_t$. A given initial probability measure for $(\eta_0, X_0)$ is then assumed.

We thus show in Section 1 that, under some hypotheses, this stochastic process is a Markov process, the marginal distributions of which, denoted by $\rho_t(i, dx)$ for $(i, t) \in E \times \mathbb{R}_+$, are measures, solution of a system of linear hyperbolic equations coupled by the boundary conditions, generalizing (1)–(3). We provide the mathematical background showing that this problem is well-posed. In particular, we prove, for the sake of completeness, a uniqueness theorem under the exact hypotheses used in this paper.

Note that numerical methods are needed to approximate these marginal distributions $\rho_t(i, dx)$, in order to evaluate the probability $A_i(t)$ that a component be in state $i \in E$ at time $t \in \mathbb{R}_+$, given by $A_i(t) = \int_{\mathbb{R}_+} \rho_t(i, dx)$. Matrix formula can be used for the computation of the Laplace transform of these marginal distributions [5], but the inversion of the Laplace transform is known to induce some stability problems. Otherwise, the phase method [1] can be used in order to boil down to a Markov process, but then the phase approximations cannot be automated. More recently, in the particular case of an alternative renewal process, computations using the resolution of Volterra equations have been proposed (see [13]) but some problems arise when the mean durations in each of the two states are contrasted.

We present in Section 2 the use of a finite volume method to find a direct approximation of these marginal distributions. Since we approximate measures instead of functions, the convergence of the method is proven thanks to compactness estimates in measure spaces, using the uniqueness result proven in Section 1.

Finally, we present some numerical results in Section 3. These results show that the finite volume method considered here is accurate and efficient in many cases.

# 1. MARGINAL DISTRIBUTIONS OF A SEMI-MARKOV PROCESS

Let us assume that the probabilistic model for the state of a component (assuming that the set $E$ of all the possible states is finite), is given by a Markov renewal process $(Y_n, T_n)_{n \in \mathbb{N}}$, where for all $n \in \mathbb{N}$, $Y_n \in E$
denotes the state of the component for any time $t$ such that $T_n \le t < T_{n+1}$. We assume that $T_0 \le 0 < T_1$ and that the semi-Markov kernel of this Markov process has a density $q$ with respect to the Lebesgue measure. These hypotheses mean that, for all $n \in \mathbb{N}$, for all $i_0, i_1, \dots, i_{n-1}, i, j \in E$, for all real values such that $s_0 \le 0 < s_1 < \dots < s_n$ and for all bounded measurable function $f$ defined on $\mathbb{R}_+$, we have

$$ 
\begin{aligned}
& \mathbb{E} \left( \begin{array}{l} 1_{\{Y_{n+1}=j\}} f(T_{n+1} - T_n) / \\ T_0 = s_0, Y_0 = i_0, T_1 = s_1, Y_1 = i_1, \dots, T_{n-1} = s_{n-1}, Y_{n-1} = i_{n-1}, T_n = s_n, Y_n = i \end{array} \right) \\
& = \mathbb{E}(1_{\{Y_{n+1}=j\}} f(T_{n+1} - T_n) / Y_n = i) \\
& = \int_{\mathbb{R}_+} f(t) q(i, j, t) \, dt.
\end{aligned}
$$

We consider the semi-Markov process $(\eta_t)_{t \in \mathbb{R}_+}$, defined by

$$ \eta_t = Y_n \quad \text{for all } n \in \mathbb{N} \text{ and for all } t \in \mathbb{R}_+ \text{ such that } T_n \le t < T_{n+1}. \quad (4) $$

Intuitively speaking, such a process $(\eta_t)_{t \in \mathbb{R}_+}$ forgets its past at times $T_n$ and satisfies

$$ 
\begin{aligned}
& \mathbb{P}(\eta_{T_{n+1}} = j, T_{n+1} - T_n \in [t, t + dt] / \eta_{T_n} = i) \\
& = \mathbb{P}(\eta_{T_1} = j, T_1 - T_0 \in [t, t + dt] / \eta_0 = i) \\
& = q(i, j, t) \, dt.
\end{aligned}
$$

We then define the transition rate $a$ by

$$ a(i, j, t) = \begin{cases} 0 & \text{if } \mathbb{P}(T_1 - T_0 > t / \eta_0 = i) = \int_t^{+\infty} \sum_{k \in E} q(i, k, u) \, du = 0, \\ \frac{q(i, j, t)}{\mathbb{P}(T_1 - T_0 > t / \eta_0 = i)} & \text{otherwise.} \end{cases} \quad (5) $$

The function $a$ intuitively satisfies the following formula

$$ \mathbb{P}(\eta_{T_1} = j, T_1 - T_0 \in [t, t + dt] / T_1 - T_0 > t, \eta_0 = i) = a(i, j, t) \, dt. $$

Defining the function $b$ by

$$ b(i, x) = \sum_{j \in E} a(i, j, x), \quad \forall (i, x) \in E \times \mathbb{R}_+, \quad (6) $$

it can be shown that

$$ \mathbb{P}(T_1 - T_0 > t / \eta_0 = i) = \exp \left( - \int_0^t b(i, u) \, du \right). \quad (7) $$

The above equation means that $b(i, \cdot)$ is the hazard rate of $T_1 - T_0$ knowing $\{\eta_0 = i\}$.

**Remark 1.1** (The case of an alternative renewal process). In order to detail the probabilistic framework in the particular case of the example given in the introduction of this paper, let us consider a component which can be in two different states: a working state denoted by 1 and a failure state denoted by 0. We suppose that the successive working periods and the successive failure periods are stochastically independent, that the working periods have the same probability density function $q(1, 0, \cdot)$ and that the failure periods have the same probability density function $q(0, 1, \cdot)$. Let $T_0 \le 0 < T_1 < \dots < T_n < \dots$ be the successive times at which the
component changes of state. Then $(T_n)_{n \ge 0}$ is an alternative renewal process and the process $\eta_t$ giving the state of the component at time $t$ is a semi-Markov process taking its values in $E = \{0, 1\}$. We have

$$b(i, t) = a(i, 1 - i, t) = \frac{q(i, 1 - i, t)}{\int_{t}^{+\infty} q(i, 1 - i, s) \, ds}, \quad \forall i = 0, 1, \forall t \in \mathbb{R}_{+}, \qquad (8)$$

or, equivalently,

$$q(i, 1 - i, t) = b(i, t) \exp \left( - \int_{0}^{t} b(i, s) \, ds \right), \quad \forall i = 0, 1, \forall t \in \mathbb{R}_{+}, \qquad (9)$$

where $a(1, 0, \cdot)$ and $a(0, 1, \cdot)$ are respectively the failure rate and the repair rate of the component.

For $k \in \mathbb{N}$, let us denote by $C^k(E \times \mathbb{R}_+)$ the space $C^k(\mathbb{R}_+)^E$. Let $X_t$ be the elapsed time in the current state at time $t$, given by

$$X_t = t - T_n, \text{ for all } t \in [T_n, T_{n+1}), \text{ for all } n \in \mathbb{N}. \qquad (10)$$

Then the process $(\eta_t, X_t)_{t \in \mathbb{R}_+}$ is a Markov process taking its values in $E \times \mathbb{R}_+$. Let $P_t$ be its semi-group, defined by

$$(P_t \xi)(i, x) = \mathbb{E}(\xi(\eta_t, X_t) / (\eta_0, X_0) = (i, x)), \quad \forall \xi \in C^0(E \times \mathbb{R}_+). \qquad (11)$$

We then have the following proposition.

**Proposition 1.2.** Let $(Y_n, T_n)_{n \ge 0}$, with $T_0 \le 0 < T_1 < \dots < T_n \dots$, be a Markov renewal process taking its values in a finite set $E$, such that its semi-Markov kernel has a density $q$ with respect to the Lebesgue measure. Let $(\eta_t)_{t \ge 0}$, $a$, $(X_t)_{t \ge 0}$ and $P_t$ be respectively defined by (4), (5), (10) and (11). We assume that $a \in C^0(E \times E \times \mathbb{R}_+)$ and that there exists $R_0 > 0$ such that $\mathbb{P}(X_0 \le R_0) = 1$. Then the following properties hold.

(1) For all $\xi \in C^0(E \times \mathbb{R}_+)$, $P_t \xi \in C^0(E \times \mathbb{R}_+)$.

(2) Let $L: C^1(E \times \mathbb{R}_+) \to C^0(E \times \mathbb{R}_+)$ be the operator defined, for all $\xi \in C^1(E \times \mathbb{R}_+)$ by

$$(i, x) \mapsto L\xi(i, x) = \sum_{j \in E} a(i, j, x)(\xi(j, 0) - \xi(i, x)) + \frac{\partial \xi}{\partial x}(i, x).$$

Then we get

$$P_t \xi = \xi + \int_0^t P_s L \xi \, ds, \quad \forall \xi \in C^1(E \times \mathbb{R}_+),$$

which leads, for all $\xi \in C^1(E \times \mathbb{R}_+)$, to

$$\mathbb{E}(\xi(\eta_t, X_t)) = \mathbb{E}(\xi(\eta_0, X_0)) + \int_0^t \mathbb{E}(L\xi(\eta_s, X_s)) \, ds.$$
(3) Therefore, the probability distribution $\rho_t$ of $(\eta_t, X_t)$ satisfies

$$ 
\begin{aligned}
\int \xi d\rho_t &= \sum_{i \in E} \int_{\mathbb{R}_+} \xi(i, x) \rho_t(i, dx) \\
&= \sum_{i \in E} \int_{\mathbb{R}_+} \xi(i, x) \rho_{\text{ini}}(i, dx) + \sum_{i \in E} \xi(i, 0) \int_0^t \sum_{j \in E} \int_{\mathbb{R}_+} a(j, i, x) \rho_s(j, dx) ds \\
&\quad - \sum_{i \in E} \int_0^t \int_{\mathbb{R}_+} b(i, x) \xi(i, x) \rho_s(i, dx) ds \\
&\quad + \sum_{i \in E} \int_0^t \int_{\mathbb{R}_+} \frac{\partial \xi}{\partial x}(i, x) \rho_s(i, dx) ds, \quad \forall \xi \in C^1(E \times \mathbb{R}_+),
\end{aligned} 
\quad (12) $$

with $\rho_{\text{ini}}(\cdot, dx) = \rho_0(\cdot, dx)$.

The proof of Proposition 1.2 is straightforward, following the ideas of [6, 7].

For $k \in \mathbb{N}$, let us denote, in the same manner as above, by $C^k(E \times \mathbb{R}_+ \times \mathbb{R}_+)$ the space $(C^k(\mathbb{R}_+ \times \mathbb{R}_+))^E$. Let $C_{ct}^k(E \times \mathbb{R}_+ \times \mathbb{R}_+)$ be the set of functions $\varphi$ of $C^k(E \times \mathbb{R}_+ \times \mathbb{R}_+)$ with a compact support in time $(i.e. \exists T > 0 : \forall t > T, \forall i \in E, \forall x \in \mathbb{R}_+, \varphi(i, x, t) = 0)$. The following result holds.

**Proposition 1.3.** Let $E$ be a finite set, let $a \in C^0(E \times E \times \mathbb{R}_+)$ be given and $b$ be defined by (6), and let $\rho_{\text{ini}}$ be a probability measure on $E \times \mathbb{R}_+$, such that there exists $R_0 > 0$ with $\sum_{i \in E} \int_{[0, R_0]} \rho_{\text{ini}}(i, dx) = 1$. Let us assume that there exists, for all $t \in \mathbb{R}_+$, a probability measure $\rho_t$ on $E \times \mathbb{R}_+$ such that (12) is satisfied. Then

$$ 
\begin{aligned}
0 = \sum_{i \in E} \int_{\mathbb{R}_+ \times \mathbb{R}_+} \left( \frac{\partial \varphi}{\partial t}(i, x, t) + \frac{\partial \varphi}{\partial x}(i, x, t) - b(i, x) \varphi(i, x, t) \right) \rho_t(i, dx) dt \\
+ \sum_{i \in E} \int_{\mathbb{R}_+} \varphi(i, x, 0) \rho_{\text{ini}}(i, dx) + \sum_{i, j \in E} \int_{\mathbb{R}_+ \times \mathbb{R}_+} a(j, i, x) \varphi(i, 0, t) \rho_t(j, dx) dt, \\
\forall \varphi \in C_{ct}^1(E \times \mathbb{R}_+ \times \mathbb{R}_+).
\end{aligned} 
\quad (13) $$

Moreover, for all $\xi \in C^0(E \times \mathbb{R})$, the function $t \mapsto \sum_{i \in E} \int_{\mathbb{R}_+} \xi(i, x) \rho_t(i, dx)$ is continuous, and $\rho_0(\cdot, dx) = \rho_{\text{ini}}(\cdot, dx)$ holds.

*Proof.* We proceed in the same way as in [7]. We first remark that the continuity of $t \mapsto \sum_{i \in E} \int_{\mathbb{R}_+} \xi(i, x) \rho_t(i, dx)$ and $\rho_0(\cdot, dx) = \rho_{\text{ini}}(\cdot, dx)$ are consequence of (12), thanks to regularizations in $C^1(E \times \mathbb{R})$ of $\xi \in C^0(E \times \mathbb{R})$. We then apply (12) to the function $\xi : (i, x) \mapsto \frac{\partial \varphi}{\partial t}(i, x, t)$, and we integrate the result on $t \in \mathbb{R}_+$. We then get (13) thanks to Fubini's theorem. $\square$

Then the following uniqueness result holds.

**Theorem 1.4.** Let $E$ be a finite set, let $a \in C^0(E \times E \times \mathbb{R}_+)$ be given and $b$ be defined by (6), and let $\rho_{\text{ini}}$ be a probability measure on $E \times \mathbb{R}_+$, such that there exists $R_0 > 0$ with $\sum_{i \in E} \int_{[0, R_0]} \rho_{\text{ini}}(i, dx) = 1$. Then there exists at most one measure $\mu$ on $E \times \mathbb{R}_+ \times \mathbb{R}_+$ such that, for all $T > 0$, there exists $R_T > 0$ with $\sum_{i \in E} \int_0^T \int_{[0, R_T]} \mu(i, dx, dt) = \sum_{i \in E} \int_0^T \int_{\mathbb{R}_+} \mu(i, dx, dt) < +\infty$ and which satisfies

$$ 
\begin{aligned}
0 = \sum_{i \in E} \int_{\mathbb{R}_+ \times \mathbb{R}_+} \left( \frac{\partial \varphi}{\partial t}(i, x, t) + \frac{\partial \varphi}{\partial x}(i, x, t) - b(i, x) \varphi(i, x, t) \right) \mu(i, dx, dt) \\
+ \sum_{i \in E} \int_{\mathbb{R}_+} \varphi(i, x, 0) \rho_{\text{ini}}(i, dx) + \sum_{i, j \in E} \int_{\mathbb{R}_+ \times \mathbb{R}_+} a(j, i, x) \varphi(i, 0, t) \mu(j, dx, dt), \\
\forall \varphi \in C_{ct}^1(E \times \mathbb{R}_+ \times \mathbb{R}_+).
\end{aligned} 
\quad (14) $$
Moreover, if, for all $t \in \mathbb{R}_+$, there exists a measure $\rho_t$ on $E \times \mathbb{R}_+$ such that for all $\xi \in C^0(E \times \mathbb{R}_+)$, the function $t \mapsto \sum_{i \in E} \int_{\mathbb{R}_+} \xi(i, x)\rho_t(i, dx)$ is continuous and $\mu(i, dx, dt) = \rho_t(i, dx)dt$, then the measure $\rho_t$ is unique too and satisfies (13).

*Proof.* We adapt the methods of [7] to the particular case considered here. Let us suppose that there exist two Radon measures satisfying the hypotheses of the above theorem. The measure $\mu$, now defined as the difference between these two solutions, thus verifies

$$0 = \sum_{i \in E} \iint_{\mathbb{R}_+ \times \mathbb{R}_+} \left( \frac{\partial \varphi}{\partial t}(i, x, t) + \frac{\partial \varphi}{\partial x}(i, x, t) - b(i, x) \varphi(i, x, t) + \sum_{j \in E} a(i, j, x) \varphi(j, 0, t) \right) \mu(i, dx, dt), \\ \forall \varphi \in C_{ct}^1(E \times \mathbb{R}_+ \times \mathbb{R}_+). \quad (15)$$

We have to prove that $\mu$ is equal to 0, which is equivalent to the relation

$$\sum_{i \in E} \iint_{\mathbb{R}_+ \times \mathbb{R}_+} \psi(i, x, t) \mu(i, dx, dt) = 0, \quad \forall \psi \in C_{ct}^0(E \times \mathbb{R}_+ \times \mathbb{R}_+). \quad (16)$$

For a given $\psi \in C_{ct}^0(E \times \mathbb{R}_+ \times \mathbb{R}_+)$, let $T > 0$ be such that $\psi(\cdot, \cdot, t) = 0$ for all $t \ge T$ and let $R_T > 0$ be such that $\sum_{i \in E} \int_0^T \int_{[0, R_T]} \mu(i, dx, dt) = \sum_{i \in E} \int_0^T \int_{\mathbb{R}} \mu(i, dx, dt) < \infty$. We prolong by continuity the function $a$ on $E \times E \times \mathbb{R}$, we define $b$ as in (6) by $b(i, x) = \sum_{j \in E} a(i, j, x)$, for all $(i, x) \in E \times \mathbb{R}$ and we prolong $\psi$ by continuity on $E \times \mathbb{R} \times \mathbb{R}_+$. Let us define the function $A: [0, T] \times C^0(E \times [-T, R_T]) \to C^0(E \times [-T, R_T])$, $(s, \xi) \mapsto A(s, \xi)$ by $(i, y) \mapsto (A(s, \xi))(i, y) = b(i, y + s) \xi(i, y) - \sum_{j \in E} a(i, j, y + s) \xi(j, -s) + \psi(i, y + s, s)$. Since the function $A$ is continuous with respect to $s$ and Lipschitz-continuous with respect to $\xi$, we deduce from the Cauchy-Lipschitz theorem the existence and the uniqueness of a function $\tilde{\varphi} \in C^0(E \times [-T, R_T] \times [0, T])$, such that $\frac{\partial \tilde{\varphi}}{\partial s} \in C^0(E \times [-T, R_T] \times [0, T])$ and

$$\begin{gathered} \tilde{\varphi}(i, y, T) = 0, \quad \forall (i, y) \in E \times [-T, R_T], \\ \frac{\partial \tilde{\varphi}}{\partial s}(i, y, s) = b(i, y + s) \tilde{\varphi}(i, y, s) - \sum_{j \in E} a(i, j, y + s) \tilde{\varphi}(j, -s, s) + \psi(i, y + s, s), \quad (17) \\ \forall (i, y, s) \in E \times [-T, R_T] \times [0, T]. \end{gathered}$$

Prolonging $\tilde{\varphi}$ in order to get that $\tilde{\varphi} \in C_{ct}^0(E \times \mathbb{R} \times \mathbb{R}_+)$ and $\frac{\partial \tilde{\varphi}}{\partial s} \in C_{ct}^0(E \times \mathbb{R} \times \mathbb{R}_+)$ (by continuity, we can take $\tilde{\varphi}(\cdot, \cdot, s) = 0$ for all $s \ge T$), let us prove that

$$0 = \sum_{i \in E} \iint_{\mathbb{R}_+ \times [0, T]} \left( \frac{\partial \tilde{\varphi}}{\partial s}(i, x - t, t) - b(i, x) \tilde{\varphi}(i, x - t, t) + \sum_{j \in E} a(i, j, x) \tilde{\varphi}(j, -t, t) \right) \mu(i, dx, dt), \quad (18)$$

which, thanks to the definition (17) of $\tilde{\varphi}$, gives (16).

Let $\beta_0$ be a real nonnegative indefinitely differentiable function with a compact support in $(-1, 1)$, such that $\int_{\mathbb{R}} \beta_0(x) dx = 1$. We consider, for all $n \in \mathbb{N}_*$, the function $\tilde{\varphi}_n \in C_{ct}^1(E \times \mathbb{R} \times \mathbb{R}_+)$ defined by $(i, y, s) \mapsto \tilde{\varphi}_n(i, y, s) = \int_{\mathbb{R}} \tilde{\varphi}(i, x, s) \beta_n(y - x) dx$, where $\beta_n(x) = n\beta_0(nx)$. We then get that $\tilde{\varphi}_n$ and $\frac{\partial \tilde{\varphi}_n}{\partial s}$ respectively converge to $\tilde{\varphi}$ and $\frac{\partial \tilde{\varphi}}{\partial s}$ in $C^0(E \times [-T, R_T] \times [0, T])$, using the uniform continuity in a compact set of any continuous function. We then apply (15) to the element $\varphi_n \in C_{ct}^1(E \times \mathbb{R}_+ \times \mathbb{R}_+)$ such that $\varphi_n(i, x, t) = \tilde{\varphi}_n(i, x - t, t)$ for all $(i, x, t) \in E \times \mathbb{R}_+ \times \mathbb{R}_+$. Since we have $\frac{\partial \varphi_n}{\partial t}(i, x, t) = -\frac{\partial \tilde{\varphi}_n}{\partial y}(i, x - t, t) + \frac{\partial \tilde{\varphi}_n}{\partial s}(i, x - t, t)$ and
$\frac{\partial \varphi_n}{\partial x}(i, x, t) = \frac{\partial \tilde{\varphi}_n}{\partial y}(i, x - t, t)$, this produces

$$0 = \sum_{i \in E} \iint_{\mathbb{R}_+ \times \mathbb{R}_+} \left( \frac{\partial \tilde{\varphi}_n}{\partial s}(i, x-t, t) - b(i, x) \tilde{\varphi}_n(i, x-t, t) + \sum_{j \in E} a(i, j, x) \tilde{\varphi}_n(j, -t, t) \right) \mu(i, dx, dt). \quad (19)$$

Passing to the limit $n \to \infty$ in (19), we thus get that (18) holds, which concludes the first part of the proof of the theorem.

Let now $\mu$ be a measure satisfying the hypotheses of the above theorem. The uniqueness, for all $t \in \mathbb{R}_+$, of a measure $\rho_t$ on $E \times \mathbb{R}_+$ such that for all $\xi \in C^0(E \times \mathbb{R})$, the function $t \mapsto \sum_{i \in E} \int_{\mathbb{R}} \xi(i, x) \rho_t(i, dx)$ is continuous and $\mu(i, dx, dt) = \rho_t(i, dx) dt$, immediately results from this continuity property since we can write

$$\sum_{i \in E} \int_{\mathbb{R}} \xi(i, x) \rho_t(i, dx) = \lim_{\tau \to 0} \sum_{i \in E} \frac{1}{\tau} \int_{\mathbb{R} \times (t, t+\tau)} \xi(i, x) \mu(i, dx, ds). \quad \square$$

Hence equation (12) appears as a weak formulation of a partial differential equation problem, the solution to which is a measure. A functional interpretation of this equation holds in the particular case where there exist some functions $\bar{u}_{\text{ini}}$ and $\bar{u}$ such that $\rho_{\text{ini}}(i, dx) = \bar{u}_{\text{ini}}(i, x)dx$ and $\rho_t(i, dx) = \bar{u}_t(i, x)dx$: in such a case, we can then identify equation (12) with a weak formulation of the system of linear scalar hyperbolic equations

$$\frac{\partial \bar{u}_t}{\partial t}(i, x) + \frac{\partial \bar{u}_t}{\partial x}(i, x) = -b(i, x) \bar{u}_t(i, x), \quad \forall (i, t, x) \in E \times \mathbb{R}_+ \times \mathbb{R}_+, \quad (20)$$

with the initial conditions

$$\bar{u}_0(i, x) = \bar{u}_{\text{ini}}(i, x), \quad \forall (i, x) \in E \times \mathbb{R}_+, \quad (21)$$

and the boundary conditions

$$\bar{u}_t(i, 0) = \sum_{j \in E} \int_{\mathbb{R}_+} a(j, i, x) \bar{u}_t(j, x) dx, \quad \forall (i, t) \in E \times \mathbb{R}_+. \quad (22)$$

(this system (20)–(22) is the generalization of Model (1)–(3) to the case of any finite set $E$ of all possible states). Then equation (12) can be obtained by multiplying (20) by any function $\xi \in C^1(E \times \mathbb{R}_+)$, and integrating by parts with respect to $x$. This interpretation is used below in order to define a numerical approximation of the solution to this equation.

## 2. THE NUMERICAL SCHEME AND ITS CONVERGENCE

### 2.1. The numerical scheme

In order to compute an approximation of the probability distribution $\rho_t$ of $(\eta_t, X_t)$ for $t \in \mathbb{R}_+$, we use a finite volume scheme (see for example [11] for a general introduction on such numerical methods). Let $h > 0$ be given and let us divide $\mathbb{R}_+$ into intervals of length $h$. The principle of the method is to approximate the measure $\rho_t(i, dx)$ by the measure $\bar{u}_t^h(i, x) dx$, where the function $\bar{u}$ is constant on each square $[mh, (m+1)h] \times [nh, (n+1)h]$ for all integers $n$ and $m$, i.e.

$$\bar{u}_t^h(i, x) = u_n^h(i, m), \quad \forall (x, t) \in [mh, (m+1)h] \times [nh, (n+1)h], \quad \forall (m, n) \in \mathbb{N} \times \mathbb{N}.$$

The finite volume scheme is inspired by the functional interpretation (20)–(22). An explicit upstream weighted finite volume scheme, in which the time step is taken equal to the maximum value authorized by a CFL condition,
can be written as

$$h \frac{u_{n+1}^h(i, m) - u_n^h(i, m)}{h} + u_n^h(i, m) - u_n^h(i, m-1) = -h b(i, mh) u_{n+1}^h(i, m) \text{ for all } m \ge 1,$$

or equivalently

$$u_{n+1}^h(i, m) = \frac{u_n^h(i, m-1)}{1 + h b(i, mh)} \quad \text{for all } m \ge 1, \qquad (23)$$

and

$$u_{n+1}^h(i, 0) = \sum_{j \in E} \sum_{m \ge 1} h a(j, i, mh) u_{n+1}^h(j, m), \qquad (24)$$

where we recall that $b(i, mh) = \sum_{j \in E} a(i, j, mh)$.

The initial condition is given by the natural discretization of measure $\rho_{\text{ini}}$:

$$u_0^h(i, m) = \frac{1}{h} \int_{[mh, (m+1)h]} \rho_{\text{ini}}(i, dx) \quad \text{for all } m \ge 0. \qquad (25)$$

**Remark 2.1** (Implicit schemes and CFL). We consider in this paper only explicit schemes under the condition CFL = 1, mainly because we want to avoid the smearing of Dirac masses during the transport step. Nevertheless, it can be necessary, for precision purposes, to choose different space steps with respect to the state (see the last numerical example presented below). In this case, it is all the same possible to generalize Scheme (23)–(25) while avoiding the smearing of Dirac masses in all the states. Indeed, let us assume that the time step is given by $h > 0$ and that, for each state $i \in E$, the space step is given by $N_i h$ with $N_i \in \mathbb{N}^*$: consider the scheme given by the initialization

$$u_0^h(i, m) = \frac{1}{N_i h} \int_{[mN_i h, (m+1)N_i h]} \rho_{\text{ini}}(i, dx) \quad \text{for all } m \ge 0, \qquad (26)$$

and, for all $n \in \mathbb{N}$ and $i \in E$,

$$ \left. \begin{array}{l} \text{if there exists } k \in \mathbb{N} \text{ such that } n+1 = kN_i, \text{ then} \\ \tilde{u}_{n+1}^h(i, 0) = 0 \text{ and } \tilde{u}_{n+1}^h(i, m) = \frac{u_n^h(i, m-1)}{1 + h b(i, mN_i h)} \text{ for all } m \ge 1, \\ \text{else} \\ \tilde{u}_{n+1}^h(i, m) = \frac{u_n^h(i, m)}{1 + h b(i, mN_i h)} \text{ for all } m \ge 0. \end{array} \right\} \qquad (27) $$

Then we set, for all $n \in \mathbb{N}$ and $i \in E$,

$$ \begin{aligned} u_{n+1}^h(i, m) &= \tilde{u}_{n+1}^h(i, m) \text{ for all } m \ge 1, \\ u_{n+1}^h(i, 0) &= \tilde{u}_{n+1}^h(i, 0) + \frac{1}{N_i} \sum_{j \in E} \sum_{m \ge 0} N_j h a(j, i, mN_j h) \tilde{u}_{n+1}^h(j, m). \end{aligned} \qquad (28) $$

Then Scheme (26)–(28) satisfies the mass conservation, i.e.

$$\sum_{i \in E} \sum_{m \in \mathbb{N}} N_i h u_{n+1}^h(i, m) = \sum_{i \in E} \sum_{m \in \mathbb{N}} N_i h u_0^h(i, m) = 1.$$

Moreover, it keeps respecting the non-diffusive effect resulting from the condition CFL= 1, even for all the states $i \in E$ such that $N_i > 1$. The convergence analysis of Scheme (26)–(28) can be done, following the method
presented in this paper for Scheme (23)–(25). It is applied in [9] to different examples, where keeping the same space step for all states would either produce inaccurate results, either lead to inadmissible computing times. In [8], we consider more general cases where the velocity can depend on the state, in multi-dimensional spaces. Note that in such cases, we cannot avoid the smearing effect of CFL< 1 and the discretization must be refined at the location where Dirac masses are transported. In such a case, implicit schemes must then be completely excluded for accuracy, but also for computing time reasons.

Some simplification of the expressions can be drawn from letting $m \in \mathbb{Z}$ instead of $m \in \mathbb{N}$. The numerical scheme (23)–(25) can then be written in the following way

$$
\begin{align*}
& a^h(j, i, m) = a(j, i, mh), \quad \forall i, j \in E, \forall m \in \mathbb{Z}, m \ge 1, \\
& a^h(j, i, m) = 0, \quad \forall i, j \in E, \forall m \in \mathbb{Z}, m < 1, \\
& b^h(i, m) = \sum_{j \in E} a^h(i, j, m), \quad \forall i \in E, \forall m \in \mathbb{Z}, \\
& u_0^h(i, m) = \frac{1}{h} \int_{[mh, (m+1)h]} \rho_{\text{ini}}(i, dx), \quad \forall i \in E, \forall m \in \mathbb{N}, \\
& u_0^h(i, m) = 0, \quad \forall i \in E, \forall m \in \mathbb{Z} \setminus \mathbb{N}, \\
& (1 + h b^h(i, m)) u_{n+1}^h(i, m) = u_n^h(i, m-1) + U_{n+1}^h(i, m), \quad \forall i \in E, \forall (m, n) \in \mathbb{Z} \times \mathbb{N}, \\
& U_{n+1}^h(i, m) = 0, \quad \forall i \in E, \forall m \in \mathbb{Z} \setminus \{0\}, \forall n \in \mathbb{N}, \\
& U_{n+1}^h(i, 0) = \sum_{j \in E} \sum_{m \in \mathbb{Z}} h a^h(j, i, m) u_{n+1}^h(j, m), \quad \forall i \in E, \forall n \in \mathbb{N}, \\
& \bar{u}_t^h(i, x) = u_n^h(i, m), \quad \forall (x, t) \in [mh, (m+1)h) \times [nh, (n+1)h), \\
& \qquad \forall i \in E, \forall (m, n) \in \mathbb{Z} \times \mathbb{N}. \tag{29}
\end{align*}
$$

In the following statements, we refer to Scheme (23)–(25) using the label Scheme (29), which gathers all the notations which are needed. We now state, in the following proposition, the property of mass conservation of this scheme.

**Proposition 2.2.** Let $E$ be a finite set, and let $a \in C_b^0(E \times E \times \mathbb{R}_+)$. Let $\rho_{\text{ini}}$ be a probability measure on the set $E \times \mathbb{R}_+$, such that there exists $R_0 > 0$ with $\sum_{i \in E} \int_{[0, R_0]} \rho_{\text{ini}}(i, dx) = 1$. Let $h > 0$ be given, and let $(u_n^h(i, m))_{(i,m,n) \in E \times \mathbb{Z} \times \mathbb{N}}$ and $\bar{u}^h : E \times \mathbb{R} \times \mathbb{R}_+ \to \mathbb{R}$ be given by Scheme (29). Then

$$
u_n^h(i, m) \ge 0, \quad \forall i \in E, \forall n \in \mathbb{N}, \forall m \in \mathbb{Z}, \tag{30}
$$

$$
u_n^h(i, m) = 0, \quad \forall i \in E, \forall n \in \mathbb{N}, \forall m \in \mathbb{Z}, m < 0, \qquad (31)
$$

$$
h \sum_{i \in E} \sum_{m \in \mathbb{Z}} u_n^h(i, m) = 1, \quad \forall n \in \mathbb{N}, \tag{32}
$$

and

$$
\sum_{i \in E} \int_{\mathbb{R}} \bar{u}_t^h(i, x) \, dx = 1, \quad \forall t \in \mathbb{R}_+. \qquad (33)
$$

Moreover, denoting, for any positive real $x$, by $[x]$ the greatest integer lower or equal to $x$, then

$$
u_n^h(i, m) = 0, \forall i \in E, \forall n \in \mathbb{N}, \forall m \in \mathbb{N}, m \gg [R_0/h] + n. \quad (34)
$$
Proof. Properties (30)–(32) hold for $n = 0$ since $\rho_{\text{ini}}$ is a probability measure with support on $E \times [0, R_0]$. Assume that these properties hold for a given $n \in \mathbb{N}$. We get from Scheme (29) that Properties (30) and (31) are verified for $m + 1$. Moreover

$$h \sum_{i \in E} \sum_{m \in \mathbb{Z}, m \ge 1} (1 + h b^h(i, m)) u_{n+1}^h(i, m) = h \sum_{i \in E} \sum_{m \in \mathbb{Z}, m \ge 0} u_n^h(i, m) = 1.$$

It delivers, since we can write

$$u_{n+1}^{h}(i, 0) = U_{n+1}^{h}(i, 0) = \sum_{j \in E} \sum_{m \in \mathbb{Z}} h a^{h}(j, i, m) u_{n+1}^{h}(j, m), \quad \forall i \in E,$$

that

$$h \sum_{i \in E} u_{n+1}^{h}(i, 0) = h \sum_{i \in E} \sum_{m \in \mathbb{Z}, m \ge 1} h b^{h}(i, m) u_{n+1}^{h}(i, m),$$

which gives

$$h \sum_{i \in E} \sum_{m \in \mathbb{Z}} u_{n+1}^{h}(i, m) = h \sum_{i \in E} \sum_{m \in \mathbb{Z}, m \ge 1} (1 + h b^{h}(i, m)) u_{n+1}^{h}(i, m) = 1.$$

This concludes the proof of (32), which immediately produces (33). The proof of (34) then comes from a proof by induction. □

## 2.2. Convergence of the numerical scheme

Let us begin with proving a property of continuity with respect to the time of the measures $\bar{u}_t^h$ dx.

**Lemma 2.3.** Let $E$ be a finite set, and let $a \in C^0(E \times E \times \mathbb{R}_+)$. Let $\rho_{\text{ini}}$ be a probability measure on the set $E \times \mathbb{R}_+$, such that there exists $R_0 > 0$ with $\sum_{i \in E} \int_{[0, R_0]} \rho_{\text{ini}}(i, dx) = 1$. Let $h \in (0, 1)$ be given, and let $(u_n^h(i, m))_{(i,m,n) \in E \times \mathbb{Z} \times \mathbb{N}}$ and $\bar{u}^h : E \times \mathbb{R} \times \mathbb{R}_+ \to \mathbb{R}$ be given by Scheme (29). Let $\psi \in C^1(E \times \mathbb{R})$. Then, for all $T \in \mathbb{R}_+$, there exists a real $C_1(\psi, T)$, which only depends on $\psi, T, a$ and $\rho_{\text{ini}}$, such that

$$\left| \sum_{i \in E} \int_{\mathbb{R}} \psi(i, x) \bar{u}_t^h(i, x) \, dx - \sum_{i \in E} \int_{\mathbb{R}} \psi(i, x) \bar{u}_s^h(i, x) \, dx \right| \le C_1(\psi, T) (|t-s| + h), \quad \forall s, t \in [0, T]. \quad (35)$$

Proof. Let us assume the hypotheses of the above lemma. Let $R_T = R_0 + T + 1$. Let $s, t$ be given such that $0 \le s \le t \le T$. Let us define

$$\Delta^h(s, t) = \sum_{i \in E} \int_{\mathbb{R}} \psi(i, x) \bar{u}_t^h(i, x) \, dx - \sum_{i \in E} \int_{\mathbb{R}} \psi(i, x) \bar{u}_s^h(i, x) \, dx.$$

We then have

$$\Delta^h(s, t) = \sum_{i \in E} \sum_{m \in \mathbb{Z}} h \bar{\psi}(i, m) u_{[t/h]}^h(i, m) - \sum_{i \in E} \sum_{m \in \mathbb{Z}} h \bar{\psi}(i, m) u_{[s/h]}^h(i, m),$$

where we set, for all $m \in \mathbb{Z}$, $\bar{\psi}(i, m) = \frac{1}{h} \int_{mh}^{(m+1)h} \psi(i, x) \, dx$. We get

$$\Delta^h(s, t) = \sum_{n=[s/h]}^{[t/h]-1} \sum_{i \in E} \sum_{m \in \mathbb{Z}} h \bar{\psi}(i, m) [u_{n+1}^h(i, m) - u_n^h(i, m)].$$
From Scheme (29), we obtain

$$
\begin{align*}
\sum_{i \in E} \sum_{m \in \mathbb{Z}} h \bar{\psi}(i, m) [u_{n+1}^h(i, m) - u_n^h(i, m)]
&= \sum_{i \in E} \sum_{m \in \mathbb{Z}} h \bar{\psi}(i, m) \left( \begin{array}{l} -u_n^h(i, m) + u_n^h(i, m-1) \\ -h b^h(i, m) u_{n+1}^h(i, m) + U_{n+1}^h(i, m) \end{array} \right) \\
&= -h \sum_{i \in E} \sum_{m \in \mathbb{Z}} [\bar{\psi}(i, m+1) - \bar{\psi}(i, m)] u_n^h(i, m) \\
&\quad - h^2 \sum_{i \in E} \sum_{m \in \mathbb{Z}} \bar{\psi}(i, m) b^h(i, m) u_{n+1}^h(i, m) \\
&\quad + h^2 \sum_{i \in E} \bar{\psi}(i, 0) \sum_{j \in E} \sum_{m \in \mathbb{Z}} a^h(j, i, m) u_{n+1}^h(j, m).
\end{align*}
$$

Hence, we deduce

$$
\begin{align*}
\Delta^h(s,t) = & -h \sum_{i \in E} \sum_{m \in \mathbb{Z}} [\bar{\psi}(i, m+1) - \bar{\psi}(i, m)] \sum_{n=[s/h]}^{[t/h]-1} u_n^h(i, m) \\
& -h^2 \sum_{i \in E} \sum_{m \in \mathbb{Z}} \bar{\psi}(i, m) b^h(i, m) \sum_{n=[s/h]}^{[t/h]-1} u_{n+1}^h(i, m) \\
& +h^2 \sum_{i \in E} \bar{\psi}(i, 0) \sum_{j \in E} \sum_{m \in \mathbb{Z}} a^h(j, i, m) \sum_{n=[s/h]}^{[t/h]-1} u_{n+1}^h(j, m).
\end{align*}
$$

Since $\bar{\psi}(i, m+1) - \bar{\psi}(i, m) = \frac{1}{h} \int_{mh}^{(m+1)h} \int_0^h \frac{\partial \psi}{\partial x}(i, x+y) \, dy \, dx$ and since $\sum_{n=[s/h]}^{[t/h]-1} u_n^h(i, m) \neq 0$ implies $0 \le mh \le (m+2)h \le R_T$, we have

$$
|\psi(i, m+1) - \psi(i, m)| \sum_{n=[s/h]}^{[t/h]-1} u_n^h(i, m) \le \max_{i \in E} \sup_{x \in [0, T_T]} \left| \frac{\partial \psi}{\partial x}(i, x) \right| \sum_{n=[s/h]}^{[t/h]-1} u_n^h(i, m).
$$

Using similar arguments for the other terms, we deduce that

$$
\begin{align*}
|\Delta^h(s,t)| \le h^2 \max_{i \in E} \sup_{x \in [0, R_T]} \left| \frac{\partial \psi}{\partial x}(i,x) \right| & \sum_{n=[s/h]}^{[t/h]-1} \sum_{i \in E} \sum_{m \in \mathbb{Z}} u_n^h(i,m) \\
& + h^2 \max_{i \in E} \sup_{x \in [0, R_T]} |\psi(i,x)| \sup_{x \in [0, R_T]} |b(i,x)| \sum_{n=[s/h]+1}^{[t/h]} \sum_{i \in E} \sum_{m \in \mathbb{Z}} u_n^h(i,m) \\
& + h^2 \max_{i \in E} \sup_{x \in [0,1]} |\psi(i,x)| \max_{i \in E} \sup_{x \in [0, R_T]} |b(i,x)| \sum_{n=[s/h]+1}^{[t/h]} \sum_{i \in E} \sum_{m \in \mathbb{Z}} u_n^h(i,m).
\end{align*}
$$

Using (32), we obtain

$$
|\Delta^h(s,t)| \leq C_1(\psi, T) h ([t/h] - [s/h]) \leq C_1(\psi, T) (t - s + h),
$$

where $C_1(\psi, T)$ is a real which only depends on $\psi, T, a$ and $\rho_{\text{ini}}$.

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

We can now prove the following theorem.
**Theorem 2.4.** Let E be a finite set, and let $a \in C^0(E \times E \times \mathbb{R}_+)$. Let $\rho_{\text{ini}}$ be a probability measure on the set $E \times \mathbb{R}_+$, such that there exists $R_0 > 0$ with $\sum_{i \in E} \int_{[0, R_0]} \rho_{\text{ini}}(i, dx) = 1$. Let $(h_p)_{p \in \mathbb{N}}$ be a sequence of strictly positive real numbers, which converges to 0 as $p \to \infty$. Then there exists a subsequence of $(h_p)_{p \in \mathbb{N}}$, again denoted $(h_p)_{p \in \mathbb{N}}$, and, for all $t \ge 0$, there exists a probability measure $\nu_t$ on $E \times \mathbb{R}_+$ with a finite support, such that the sequence of measures $((\bar{u}_t^{h_p})(i, x) dx, i \in E))_{p \in \mathbb{N}}$ converges to $\nu_t$ as $p \to \infty$, where $\bar{u}^{h_p} : E \times \mathbb{R} \times \mathbb{R}_+ \to \mathbb{R}$ is given by Scheme (29) for $h = h_p$. The measure convergence means that for all $\xi \in C^0(E \times \mathbb{R}_+)$, we have

$$ \lim_{p \to +\infty} \sum_{i \in E} \int_{E \times \mathbb{R}_+} \xi(i, x) \bar{u}_t^{h_p}(i, x) \, dx = \int_{E \times \mathbb{R}_+} \xi \, d\nu_t. $$

Moreover, for all $\xi \in C^0(E \times \mathbb{R}_+)$, the function $t \mapsto \int_{E \times \mathbb{R}_+} \xi \, d\nu_t$ is continuous.

*Proof.* Without lack of generality, we assume that $h_p \in (0, 1)$ for all $p \in \mathbb{N}$. Let $(t_m)_{m \in \mathbb{N}}$ be a sequence of real numbers, which is dense into $\mathbb{R}_+$. For any $m \in \mathbb{N}$ and $p \in \mathbb{N}$, we deduce from (32) and (34) that the measure $\bar{u}_{t_m}^{h_p} dx$ is a probability measure on $E \times \mathbb{R}_+$ with a support included in $[0, R_0 + t_m + 1]$. Thus, using Prohorov's theorem, we deduce that there exists a subsequence of $(h_p)_{p \in \mathbb{N}}$ (depending on $m$) and a probability measure $\nu_{t_m}$ such that the sequence of measures $((\bar{u}_{t_m}^{h_p})(i, x) dx, i \in E))_{p \in \mathbb{N}}$ converges to $\nu_{t_m}$ as $p \to \infty$ in the following sense: for all bounded functions $\psi \in C^0(E \times \mathbb{R}_+)$:

$$ \lim_{p \to +\infty} \sum_{i \in E} \int_{E \times \mathbb{R}_+} \psi(i, x) \bar{u}_{t_m}^{h_p}(i, x) \, dx = \int_{E \times \mathbb{R}_+} \psi \, d\nu_{t_m}. $$

Using the diagonal method, we then can build a subsequence, again denoted $(h_p)_{p \in \mathbb{N}}$ such that for all $m \in \mathbb{N}$, the sequence $(\bar{u}_{t_m}^{h_p}(i, x) dx)_{p \in \mathbb{N}}$ converges to the measure $\nu_{t_m}$ in the sense given above. Let $t \in \mathbb{R}_+$ and let $\xi \in C^0(E \times \mathbb{R}_+)$ be given. Let us prove that $(\sum_{i \in E} \int_{\mathbb{R}} \xi(i, x) \bar{u}_t^{h_p}(i, x) dx)_{p \in \mathbb{N}}$ converges as $p \to \infty$. To that purpose, we prove that it is a Cauchy sequence. We set $T = t+1$ and $R_T = R_0 + T + 1$. Let $\varepsilon > 0$ be given. We approximate function $\xi$ by a bounded function $\psi \in C^1(E \times \mathbb{R}_+)$ such that $\sup_{x \in [0, R_T]} |\xi(x) - \psi(x)| \le \varepsilon$. Let $m \in \mathbb{N}$ be such that $|t_m - t| \le \min(\varepsilon/C_1(\psi, T), 1)$. We then obtain for all $p \in \mathbb{N}$, using Lemma 2.3,

$$ \left| \sum_{i \in E} \int_{\mathbb{R}} \psi(i, x) (\bar{u}_t^{h_p}(i, x) - \bar{u}_{t_m}^{h_p}(i, x)) dx \right| \le C_1(\psi, T) (|t - t_m| + h_p). $$

We now choose $n \in \mathbb{N}$ such that, for all $p, q \ge n$, $\left| \sum_{i \in E} \int_{\mathbb{R}} \psi(i, x) (\bar{u}_{t_m}^{h_p}(i, x) - \bar{u}_{t_m}^{h_q}(i, x)) dx \right| \le \varepsilon$ (this is possible since the sequence $(\sum_{i \in E} \int_{\mathbb{R}} \psi(i, x) \bar{u}_{t_m}^{h_p}(i, x) dx)_{p \in \mathbb{N}}$ converges) and $h_p$ and $h_q$ less or equal to $\varepsilon/C_1(\psi, T)$. We then get, for $p, q \ge n$,

$$ \begin{aligned} \left| \sum_{i \in E} \int_{\mathbb{R}} \xi(i, x) (\bar{u}_t^{h_p}(i, x) - \bar{u}_t^{h_q}(i, x)) dx \right| & \le \left| \sum_{i \in E} \int_{\mathbb{R}} \psi(i, x) (\bar{u}_t^{h_p}(i, x) - \bar{u}_t^{h_q}(i, x)) dx \right| + 2\varepsilon \\ & \le \left| \sum_{i \in E} \int_{\mathbb{R}} \psi(i, x) (\bar{u}_t^{h_p}(i, x) - \bar{u}_{t_m}^{h_p}(i, x)) dx \right| \\ & \qquad + \left| \sum_{i \in E} \int_{\mathbb{R}} \psi(i, x) (\bar{u}_{t_m}^{h_p}(i, x) - \bar{u}_{t_m}^{h_q}(i, x)) dx \right| \\ & \qquad + \left| \sum_{i \in E} \int_{\mathbb{R}} \psi(i, x) (\bar{u}_t^{h_q}(i, x) - \bar{u}_{t_m}^{h_q}(i, x)) dx \right| + 2\varepsilon \\ & \le 7\varepsilon. \end{aligned} $$
This proves that $(\sum_{i \in E} \int_{\mathbb{R}} \xi(i, x) \bar{u}_t^{h_p}(i, x) dx)_{p \in \mathbb{N}}$ is a Cauchy sequence. Let us define

$$S_t(\xi) = \lim_{p \to \infty} \sum_{i \in E} \int_{\mathbb{R}} \xi(i, x) \bar{u}_t^{h_p}(i, x) dx.$$

The operator $S_t$ is positive linear form on $C^0(E \times \mathbb{R}_+)$ and $S_t(1) = 1$. Consequently it is a probability denoted by $\nu_t$, i.e. $S_t(\xi) = \int \xi \, d\nu_t = \sum_{i \in E} \int_{\mathbb{R}_+} \xi(i, x) \nu_t(i, dx)$.

Let us prove the continuity of $t \mapsto \int \xi \, d\nu_t$, for any $\xi \in C^0(E \times \mathbb{R} \times \mathbb{R}_+)$. We again set $T = t+1$ and $R_T = R_0+T+1$. Let $\varepsilon > 0$. We again approximate $\xi$ by $\psi \in C^1(E \times \mathbb{R}_+)$ such that $\sup_{x \in [0, R_T]} |\xi(x) - \psi(x)| \le \varepsilon$. We let $h = h_p$ in (35), and we let $p \to \infty$. We thus get, for $s \in [0, T]$,

$$\left| \sum_{i \in E} \int_{\mathbb{R}} \psi(i, x) \rho_t(i, dx) - \sum_{i \in E} \int_{\mathbb{R}} \psi(i, x) \rho_s(i, dx) \right| \le C_1(\psi, T) |t-s|.$$

We then have, for all $s \in [0, T]$ such that $|t-s| \le \varepsilon/C_1(\psi, T)$,

$$\left| \sum_{i \in E} \int_{\mathbb{R}} \xi(i, x) \nu_t(i, dx) - \sum_{i \in E} \int_{\mathbb{R}} \xi(i, x) \nu_s(i, dx) \right| \le \left| \sum_{i \in E} \int_{\mathbb{R}} \psi(i, x) \nu_t(i, dx) - \sum_{i \in E} \int_{\mathbb{R}} \psi(i, x) \nu_s(i, dx) \right| + 2\varepsilon \\ \le 3\varepsilon,$$

which completes this proof of continuity.

□

We now get the following corollary.

**Corollary 2.5.** Let $E$ be a finite set, and let $a \in C^0(E \times E \times \mathbb{R}_+)$. Let $\rho_{\text{ini}}$ be a probability measure on the set $E \times \mathbb{R}_+$, such that there exists $R_0 > 0$ with $\sum_{i \in E} \int_{[0, R_0]} \rho_{\text{ini}}(i, dx) = 1$. Let $(h_p)_{p \in \mathbb{N}}$ be a sequence of strictly positive real numbers, which converges to 0 as $p \to \infty$. Let us again denote by $(h_p)_{p \in \mathbb{N}}$ the subsequence of $(h_p)_{p \in \mathbb{N}}$ and, for all $t \ge 0$, by $\nu_t$ the probability measure with a finite support given by Theorem 2.4. Then the sequence of measures $(\bar{u}_t^{h_p}(i, x) \, dx \, dt)_{p \in \mathbb{N}}$, where $\bar{u}_p^{h_p} : E \times \mathbb{R} \times \mathbb{R}_+ \to \mathbb{R}$ is given by Scheme (29) for $h = h_p$, converges to the measure $\mu$ defined by $\mu(i, dx, dt) = \nu_t(i, dx) \, dt$ as $p \to \infty$ in the following sense: for all $\xi \in C_{0,ct}(E \times \mathbb{R} \times \mathbb{R}_+)$ we have

$$\lim_{p \to \infty} \sum_{i \in E} \int_{\mathbb{R}_+} \int_{\mathbb{R}_+} \xi(i, x, t) \bar{u}_t^{h_p}(i, x) dx dt = \sum_{i \in E} \int_{\mathbb{R}_+} \int_{\mathbb{R}_+} \xi(i, x, t) \nu_t(i, dx) dt = \int_{E \times \mathbb{R}_+ \times \mathbb{R}_+} \xi \, d\mu.$$

Moreover, let $\xi \in C_{0,ct}(E \times \mathbb{R} \times \mathbb{R}_+)$ and $A$ be such that $\xi(i, x, t) = 0, \forall t > A, \forall i \in E, \forall x \in \mathbb{R}$, and let $R_A = R_0 + A + 1$. We then have

$$\lim_{p \to +\infty} h_p^2 \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=0}^{[A/h_p]-1} \xi(i, mh, nh) u_n^h(i, m) = \int_{E \times \mathbb{R}_+ \times \mathbb{R}_+} \xi \, d\mu. \quad (36)$$

*Proof.* Let $\xi \in C^0(E \times \mathbb{R} \times \mathbb{R}_+)$. We get that

$$\lim_{p \to \infty} \sum_{i \in E} \int_{\mathbb{R}_+} \int_{\mathbb{R}_+} \xi(i, x, t) \bar{u}_t^{h_p}(i, x) dx dt = \sum_{i \in E} \int_{\mathbb{R}_+} \int_{\mathbb{R}_+} \xi(i, x, t) \rho_t(i, dx) dt,$$

thanks to the Lebesgue's dominated convergence theorem.
Let us now prove (36). We can write

$$
\begin{align*}
\sum_{i \in E} \int_{\mathbb{R}_+ \times \mathbb{R}_+} \xi(i, x, t) \bar{u}_t^h(i, x) \, dx \, dt &= \sum_{i \in E} \sum_{m \ge 0} \sum_{n \ge 0} \iint_{[mh,(m+1)h \times [nh,(n+1)h[} \xi(i, x, t) u_n^h(i, m) \, dx \, dt \\
&= \sum_{i \in E} \sum_{m \ge 0} \sum_{n=0}^{[A/h]-1} u_n^h(i, m) \iint_{[mh,(m+1)h \times [nh,(n+1)h[} \xi(i, x, t) \, dx \, dt \\
&\quad + \sum_{i \in E} \sum_{m \ge 0} u_n^h(i, m) \iint_{[mh,(m+1)h \times [h[A/h],A]} \xi(i, x, t) \, dx \, dt.
\end{align*}
$$

On the one hand:

$$
\left| \sum_{i \in E} \sum_{m \ge 0} u_n^h(i, m) \iint_{[mh,(m+1)h \times [h[A/h],A]} \xi(i, x, t) \, dx \, dt \right| \le \sum_{i \in E} \sum_{m \ge 0} u_n^h(i, m) \|\xi\|_{L^\infty(E \times \mathbb{R}_+ \times \mathbb{R}_+)} h(A-h[A/h]) \\
= \|\xi\|_{L^\infty(E \times \mathbb{R}_+ \times \mathbb{R}_+)} (A-h[A/h])
$$

and the right hand side of the above expression tends to 0 when $h$ tends to 0.

On the other hand, we have

$$
\begin{align*}
\delta_{\xi}(h) = \left| h^2 \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=0}^{[A/h]-1} \xi(i, mh, nh) u_n^h(i, m) \right. \\
\left. - \sum_{i \in E} \sum_{m \ge 0} \sum_{n=0}^{[A/h]-1} u_n^h(i, m) \iint_{[mh,(m+1)h \times [nh,(n+1)h[} \xi(i, x, t) \, dx \, dt \right| \\
\le \sum_{i \in E} \sum_{m \ge 0} \sum_{n=0}^{[A/h]-1} u_n^h(i, m) \iint_{[mh,(m+1)h \times [nh,(n+1)h[} |\xi(i, mh, nh) - \xi(i, x, t)| \, dx \, dt.
\end{align*}
$$

Using (34) we obtain

$$
\delta'_{\xi}(h) \le \sum_{i \in E} \sum_{m=0}^{[(R_0+A)/h]-1} \sum_{n=0}^{[A/h]-1} u_n^h(i, m) \iint_{[mh,(m+1)h \times [nh,(n+1)h[} |\xi(i,mh,nh) - \xi(i,x,t)| \, dx \, dt.
$$

Let us define

$$
\epsilon_{A,\xi}(h) = \sup\{|\xi(i, x_1, t_1) - \xi(i, x_2, t_2)| : i \in E, 0 \le x_1, x_2 \le R_0 + A, 0 \le t_1, t_2 \le T, \\
|x_1 - x_2| \le h, |t_1 - t_2| \le h\}.
$$

Then we have

$$
0 \le \delta_{\xi}(h) \le \sum_{i \in E} \sum_{m=0}^{[(R_0+A)/h]-1} \sum_{n=0}^{[A/h]-1} \epsilon_{A,\xi}(h) h2 u_n^h(i,m) \le \epsilon_{A,\xi}(h) [A/h] \le A \epsilon_{A,\xi}(h)
$$

and the right hand side of the above inequalities tends to 0 when $h$ tends to $O$, which ends the proof.

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

Let us now prove that the measure $\mu(i, dx, dt) = \nu_t(i, dx) dt$ satisfies (14).
**Proposition 2.6.** Let E be a finite set, and let $a \in C^0(E \times E \times \mathbb{R}_+)$. Let $\rho_{\text{ini}}$ be a probability measure on the set $E \times \mathbb{R}_+$, such that there exists $R_0 > 0$ with $\sum_{i \in E} \int_{[0,R_0]} \rho_{\text{ini}}(i, dx) = 1$. Let $(h_p)_{p \in \mathbb{N}}$ be a sequence of strictly positive real numbers, which converges to 0 as $p \to \infty$. Let us again denote by $(h_p)_{p \in \mathbb{N}}$ the subsequence of $(h_p)_{p \in \mathbb{N}}$ and, for all $t \ge 0$, by $\nu_t$ the probability measure with a finite support given by Theorem 2.4. Then the measure $\mu(i, dx, dt) = \rho_t(i, dx) dt$ is such that for all $T > 0$, there exists $R_T > 0$ with $\sum_{i \in E} \int_0^T \int_{[0,R_T]} \mu(i, dx, dt) = \sum_{i \in E} \int_0^T \int_{\mathbb{R}_+} \mu(i, dx, dt) < \infty$ and that (14) is satisfied.

*Proof.* Let us assume the assumptions of the above proposition, and let $\varphi \in C_{ct}^1(E \times \mathbb{R} \times \mathbb{R}_+)$. Without lack of generality, we can assume that $h_p \in (0, 1)$. We consider a time $T > 1$ such that

$$ \varphi(i, x, t) = 0 \quad \forall (i, x, t) \in E \times \mathbb{R} \times [T-1, +\infty), \quad (37) $$

and we set $R_T = R_0 + T + 1$. Let us simplify the notation in all the following proof, replacing, in any equation where $p \in \mathbb{N}$ is given, $h_p$ by $h$. We also set

$$ \mu^h(i, dx, dt) = \bar{u}_t^h(i, x) dx dt, \quad \forall (i, x, t) \in E \times \mathbb{R} \times \mathbb{R}_+. \quad (38) $$

Let us notice that (33) gives

$$ \int_{E \times \mathbb{R}} \times [0, T] d\mu = \int_{E \times \mathbb{R}_+} \times [0, T] d\mu = T. \quad (39) $$

Let $\epsilon \in (0, 1)$ and $h \in (0, 1 - \epsilon)$. Let us first prove that

$$ \int_{E \times \mathbb{R} \times [0, T]} d\mu^h = \int_{E \times \mathbb{R}_+ \times [0, T]} d\mu^h = \int_{E \times [0, R_T - \epsilon] \times [0, T]} d\mu^h, \quad (40) $$

$$ \int_{E \times \mathbb{R}_+ \times [0, T]} d\mu = \int_{E \times [0, R_T] \times [0, T]} d\mu. $$

We have

$$ 0 \le \int_{E \times [R_T - \epsilon, +\infty \times [0, T]} d\mu^h \le \sum_{m \ge \frac{R_T - \epsilon}{h}} \sum_{n=0}^{[T/h]} h^2 u_n^h(i, m). $$

Then Formula (34) gives $u_n^h(i, m) = 0$ for $0 \le n \le [T/h]$, and $m > [R_0/h] + [T/h]$. It is easy to check that $[(R_T - \epsilon)/h] < [R_0/h] + [T/h]$. Therefore we get

$$ \int_{E \times [R_T - \epsilon, +\infty \times [0, T]} d\mu^h = 0, $$

which proves the result for measure $\mu^h$.

Now, let $\epsilon > 0$ and let $\varphi_\epsilon \in C^0(E \times \mathbb{R}_+ \times \mathbb{R}_+)$ a function which satisfies: $0 \le \varphi_\epsilon \le 1$, $\varphi_\epsilon = 1$ on $E \times [R_T, +\infty \times [0, T]$ and $\varphi_\epsilon = 0$ on $E \times [R_T - \epsilon, +\infty \times [0, T + \epsilon]$. We obtain

$$ 0 \le \int_{E \times [R_T, +\infty \times [0, T]} d\mu \le \int \varphi_\epsilon d\mu = \lim_{h \to 0} \int \varphi_\epsilon d\mu^h. $$

For $h < 1 - \epsilon$, we have

$$ 0 \le \int \varphi_\epsilon d\mu^h \le \int_{E \times [R_T - \epsilon, +\infty \times [T, T + \epsilon]} d\mu^h \le \int_T^{T+\epsilon} \left( \sum_{i \in E} \int_{R_+} \bar{u}_t^h(i, x) dx \right) dt \le \epsilon. $$
Let us now prove that measure $\mu$ satisfies (14). We get, from Scheme (29),

$$0 = u_{n+1}^h(i, m) - u_n^h(i, m) + u_n^h(i, m) - u_n^h(i, m-1) \\ + h b^h(i, m) u_{n+1}^h(i, m) - U_{n+1}^h(i, m).$$

Multiplying the above equation by $-h \varphi(i, mh, nh)$ and summing on $i \in E, m \in \mathbb{Z}$ and $0 \le n \le [T/h] - 1$, we obtain

$$T_1^h + T_2^h + T_3^h + T_4^h = 0,$$

with

$$T_1^h = - \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=0}^{[T/h]-1} (u_{n+1}^h(i, m) - u_n^h(i, m)) h \varphi(i, mh, nh),$$

$$T_2^h = - \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=0}^{[T/h]-1} (u_n^h(i, m) - u_n^h(i, m-1)) h \varphi(i, mh, nh),$$

$$T_3^h = - \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=0}^{[T/h]-1} h^2 b^h(i, m) u_{n+1}^h(i, m) \varphi(i, mh, nh),$$

$$T_4^h = \sum_{i \in E} \sum_{n=0}^{[T/h]-1} h U_{n+1}^h(i, 0) \varphi(i, 0, nh).$$

### Step 1. Study of $T_1^h$

We have

$$T_1^h = -h \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=1}^{[T/h]} u_n^h(i, m) \varphi(i, mh, (n-1)h) \\ + h \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=0}^{[T/h]-1} u_n^h(i, m) \varphi(i, mh, nh) \\ = h \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=1}^{[T/h]-1} u_n^h(i, m) (\varphi(i, mh, nh) - \varphi(i, mh, (n-1)h)) \\ + h \sum_{i \in E} \sum_{m \in \mathbb{Z}} u_0^h(i, m) \varphi(i, mh, 0),$$

since, if $h$ is small enough, i.e. if $p$ is large enough, we have $\varphi(i, mh, ([T/h] - 1)h) = 0$. Let us define $T_1^{p}$ and $T_1^{pp}$ by

$$T_1^{p} = h \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=1}^{[T/h]-1} u_n^h(i, m) (\varphi(i, mh, nh) - \varphi(i, mh, (n-1)h)), \\ T_1^{pp} = h \sum_{i \in E} \sum_{m \in \mathbb{Z}} u_0^h(i, m) \varphi(i, mh, 0).$$
We also define $\psi^h$ by

$$ \begin{aligned} \psi^h(i, m, n) &= \frac{\varphi(i, mh, nh) - \varphi(i, mh, (n-1)h)}{h} \quad \text{if } n \neq 0, \\ \psi^h(i, m, 0) &= \frac{\partial \varphi}{\partial t}(i, mh, 0). \end{aligned} $$

and $\bar{\psi}^h$ by

$$ \bar{\psi}^h(i, x, t) = \psi^h(i, m, n) \text{ if } (x, t) \in [mh, (m+1)h] \times [nh, (n+1)h], \forall m \in \mathbb{Z}, 0 \le n \le [T/h] - 1. $$

For $n \ge 1$, we notice that there exists $t_n^h \in [(n-1)h, nh]$ such that

$$ \psi^h(i, m, n) = \frac{\partial \varphi}{\partial t}(i, mh, t_n^h). $$

Since the function $\frac{\partial \varphi}{\partial t}$ is continuous, it is uniformly continuous on $[0, R_T] \times [0, T]$, and therefore

$$ \lim_{p \to \infty} \left\| \frac{\partial \varphi}{\partial t} - \bar{\psi}^{h_p} \right\|_{L^{\infty}(E \times [0, R_T] \times [0, T])} = 0. \quad (41) $$

We have

$$ \begin{aligned} \int_{E \times \mathbb{R} \times \mathbb{R}_+} \bar{\psi}^h \, d\mu^h &= \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=0}^{[T/h]-1} h^2 \psi^h(i, m, n) u_n^h(i, m) + \int_{E \times \mathbb{R} \times [h[T/h], T]} \bar{\psi}^h \, d\mu^h \\ &= T_1^{p} + h^2 \sum_{i \in E} \sum_{m \in \mathbb{Z}} u_0^h(i, m) \psi^h(i, m, 0) + \int_{E \times \mathbb{R} \times [h[T/h], T]} \bar{\psi}^h \, d\mu^h. \end{aligned} \quad (42) $$

On the one hand, using (37), (40) and (39), we get

$$ \begin{aligned} \left| \int_{E \times \mathbb{R} \times \mathbb{R}_+} \bar{\psi}^h \, d\mu^h - \int_{E \times \mathbb{R} \times \mathbb{R}_+} \frac{\partial \varphi}{\partial t} \, d\mu \right| &\le \left| \int_{E \times \mathbb{R} \times [0, T)} \left( \bar{\psi}^h - \frac{\partial \varphi}{\partial t} \right) \, d\mu^h \right| \\ &\quad + \left| \int_{E \times \mathbb{R} \times \mathbb{R}_+} \frac{\partial \varphi}{\partial t} \, d\mu^h - \int_{E \times \mathbb{R} \times \mathbb{R}_+} \frac{\partial \varphi}{\partial t} \, d\mu \right| \\ &\le T \left\| \frac{\partial \varphi}{\partial t} - \bar{\psi}^h \right\|_{L^{\infty}(E \times [0, R_T] \times [0, T])} \\ &\quad + \left| \int_{E \times \mathbb{R} \times \mathbb{R}_+} \frac{\partial \varphi}{\partial t} \, d\mu^h - \int_{E \times \mathbb{R} \times \mathbb{R}_+} \frac{\partial \varphi}{\partial t} \, d\mu \right|. \end{aligned} $$

Since $\frac{\partial \varphi}{\partial t} \in C0(E \times \mathbb{R} \times \mathbb{R}_+)$, Corollary 2.5 and (41) imply that

$$ \lim_{p \to \infty} \int_{E \times \mathbb{R} \times [0, T)} \bar{\psi}^{h_p} \, d\mu^{h_p} = \int_{E \times \mathbb{R} \times [0, T)} \frac{\partial \varphi}{\partial t} \, d\mu. \quad (43) $$
On the other hand, using (32) and (39) we obtain

$$
\begin{aligned}
& \left| h^2 \sum_{i \in E} \sum_{m \in \mathbb{Z}} u_0^h(i, m) \psi^h(i, m, 0) + \int_{E \times \mathbb{R} \times [h[T/h], T]} \bar{\psi}^h \, d\mu^h \right| \\
& \leq h \|\bar{\psi}^h\|_{L^\infty(E \times [0, R_T] \times [0, T])} h \sum_{i \in E} \sum_{m \in \mathbb{Z}} u_0^h(i, m) + h \|\bar{\psi}^h\|_{L^\infty(E \times [0, R_T] \times [0, T])} \\
& = 2h \|\bar{\psi}^h\|_{L^\infty(E \times [0, R_T] \times [0, T])}.
\end{aligned}
$$

But (41) implies that $\|\bar{\psi}^h\|_{L^\infty(E \times [0, R_T] \times [0, T])}$ is bounded. Consequently we obtain

$$
\lim_{p \to \infty} \left( h_p^2 \sum_{i \in E} \sum_{m \in \mathbb{Z}} u_0^{h_p}(i, m) \psi^{h_p}(i, m, 0) + \int_{E \times \mathbb{R} \times [h_p[T/h_p], T]} \bar{\psi}^{h_p} \, d\mu^{h_p} \right) = 0. \quad (44)
$$

Equations (42)-(44) yield

$$
\lim_{p \to \infty} T_1^{p} = \int \frac{\partial \varphi}{\partial t} \, d\mu = \int_{E \times \mathbb{R}_+ \times \mathbb{R}_+} \frac{\partial \varphi}{\partial t} \, d\mu . \quad (45)
$$

To study $T_1^{p''}$, let us notice that if $f$ is a real-valued continuous function, we have

$$
\lim_{p \to \infty} h_p \sum_{m \in \mathbb{Z}} f(mh_p) u_0^{h_p}(i, m) = \int_{\mathbb{R}} f(x) \rho_{\text{ini}}(i, dx).
$$

Indeed, we can write

$$
\begin{aligned}
\left| \int_{\mathbb{R}} f(x) \rho_{\text{ini}}(i, dx) - h \sum_{m \in \mathbb{Z}} f(mh) u_0^h(i, m) \right| &= \left| \int_{\mathbb{R}} f(x) \rho_{\text{ini}}(i, dx) - \sum_{m \in \mathbb{Z}} \int_{[mh,(m+1)h]} f(mh) \rho_{\text{ini}}(i, dx) \right| \\
&\le \sum_{m \in \mathbb{Z}} \int_{[mh,(m+1)h]} |f(x) - f(mh)| \rho_{\text{ini}}(i, dx) \\
&\le \sup_{\substack{x_1, x_2 \in [0, R_0] \\ |x_1 - x_2| \le h}} |f(x_1) - f(x_2)| \sum_{m \in \mathbb{Z}} \int_{[mh,(m+1)h]} \rho_{\text{ini}}(i, dx) \\
&= \sup_{\substack{x_1, x_2 \in [0, R_0] \\ |x_1 - x_2| \le h}} |f(x_1) - f(x_2)|,
\end{aligned}
$$

and the uniform continuity of $f$ on $[0, R_0]$ implies the convergence to 0 when $p$ tends to $\infty$. Therefore we get

$$
\lim_{p \to \infty} T_1^{p''} = \sum_{i \in E} \int_{\mathbb{R}} \varphi(i, x, 0) \rho_{\text{ini}}(i, dx).
$$

Gathering the previous results, we have proven that

$$
\lim_{p \to \infty} T_1^h = \int_{E \times \mathbb{R}_+ \times \mathbb{R}_+} \frac{\partial \varphi}{\partial t} \, d\mu + \sum_{i \in E} \int_{\mathbb{R}_+} \varphi(i, x, 0) \rho_{\text{ini}}(i, dx). \quad (46)
$$
### Step 2. Study of $T_2^h$

We have

$$ 
\begin{aligned}
T_2^h &= -h \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=0}^{[T/h]-1} u_n^h(i, m) \varphi(i, mh, nh) \\
&\quad + h \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=0}^{[T/h]-1} u_n^h(i, m) \varphi(i, (m+1)h, nh) \\
&= h \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=0}^{[T/h]-1} u_n^h(i, m) (\varphi(i, (m+1)h, nh) - \varphi(i, mh, nh)).
\end{aligned}
$$

Using the same methods as in step 1, we show that

$$ \lim_{p \to \infty} T_2^h = \int_{E \times \mathbb{R}_+ \times [0, T)} \frac{\partial \varphi}{\partial x} d\mu. $$

### Step 3. Study of $T_3^h$

We have already seen that, if $h$ is small enough, $\varphi(i, mh, ([T/h] - 1)h) = 0$, consequently we get

$$ T_3^h = -h^2 \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=1}^{[T/h]-1} b^h(i, m) u_n^h(i, m) \varphi(i, mh, (n-1)h). $$

We notice that

$$ \lim_{h \to 0} T_3^h = \lim_{h \to 0} \tilde{T}_3^h $$

with

$$ \tilde{T}_3^h = -h^2 \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=1}^{[T/h]-1} b^h(i, m) u_n^h(i, m) \varphi(i, mh, nh), $$

since, using (32) and (34) we obtain

$$ 
\begin{aligned}
& h^2 \sum_{i \in E} \sum_{m \in \mathbb{Z}} \sum_{n=1}^{[T/h]-1} b^h(i, m) u_n^h(i, m) |\varphi(i, mh, nh) - \varphi(i, mh, (n-1)h)| \\
&= h^2 \sum_{i \in E} \sum_{\substack{m \in \mathbb{Z}, \\ mh \le R_T}} \sum_{n=1}^{[T/h]-1} b^h(i, m) u_n^h(i, m) |\varphi(i, mh, nh) - \varphi(i, mh, (n-1)h)| \\
&\le \sup_{x \in [0, R_T]} \sup_{|t_1 - t_2| \le h} |\varphi(i, x, t_1) - \varphi(i, x, t_2)| \|b\|_{L^\infty(E \times [0, R_T])} h^2 \sum_{i \in E} \sum_{n=0}^{[T/h]-1} \sum_{m \in \mathbb{Z}} u_n^h(i, m) \\
&\le \sup_{x \in [0, R_T]} \sup_{|t_1 - t_2| \le h} |\varphi(i, x, t_1) - \varphi(i, x, t_2)| \|b\|_{L^\infty(E \times [0, R_T])} T,
\end{aligned}
$$

and the uniform continuity of $\varphi$ on $[0, R_T] \times [0, T]$ gives the result.
In fact, the following relation holds

$$
\begin{aligned}
\tilde{T}_3^h = & -h^2 \sum_{i \in E} \sum_{m \ge 0} \sum_{n=0}^{[T/h]-1} b^h(i, m) u_n^h(i, m) \varphi(i, mh, nh) \\
& + h^2 \sum_{i \in E} \sum_{m \in \mathbb{Z}} b^h(i, m) u_0^h(i, m) \varphi(i, mh, 0).
\end{aligned}
$$

On the one hand, we have, using (36)

$$
\lim_{p \to \infty} h_p^2 \sum_{i \in E} \sum_{m \ge 0} \sum_{n=0}^{[T/h_p]-1} b^{h_p}(i,m) \varphi(i,mh_p,nh_p) u_n^{h_p}(i,m) = \int_{E \times \mathbb{R}_+ \times \mathbb{R}_+} b \varphi \, d\mu.
$$

On the other hand, thanks to (32) and (34), we have

$$
\begin{aligned}
& |h^2 \sum_{i \in E} \sum_{m \in \mathbb{Z}} b^h(i, m) u_0^h(i, m) \varphi(i, mh, 0)| \\
& \le \|b\|_{L^\infty([0, R_T])} \|\varphi\|_{L^\infty(E \times [0, R_T] \times [0, T])} h^2 \sum_{i \in E} \sum_{m \in \mathbb{Z}} u_0^h(i, m) \\
& = \|b\|_{L^\infty([0, R_T])} \|\varphi\|_{L^\infty(E \times [0, R_T] \times [0, T])} h,
\end{aligned}
$$

therefore this term tends to 0 when p tends to $\infty$.

We have therefore shown that

$$
\lim_{p \to \infty} T_3^h = - \int_{E \times \mathbb{R}_+ \times \mathbb{R}_+} b \varphi \, d\mu
$$

**Step 4.** Study of $T_4^h$

We have

$$
T_4^h = \sum_{i \in E} \sum_{n=0}^{[T/h]-1} h \sum_{j \in E} \sum_{m \in \mathbb{Z}} h a^h(j, i, m) u_{n+1}^h(j, m) \varphi(i, 0, nh).
$$

and therefore as previously

$$
\lim_{h \to 0} T_4^h = \sum_{i \in E} \sum_{j \in E} \iint_{\mathbb{R} \times [0, T)} a(i, j, x) \varphi(j, 0, t) \mu(i, dx, dt)
$$

**Step 5.** Synthesis

Since

$$
0 = T_{1}^{h} + T_{2}^{h} + T_{3}^{h} + T_{4}^{h},
$$
we have

$$
\begin{aligned}
0 &= \lim_{h \to 0} (T_1^h + T_2^h + T_3^h + T_4^h) \\
&= \sum_{i \in E} \iint_{\mathbb{R}_+ \times \mathbb{R}_+} \frac{\partial \varphi}{\partial t}(i, x, t) \mu(i, dx, dt) + \sum_{i \in E} \int_{\mathbb{R}_+} \varphi(i, x, 0) \rho_{\text{ini}}(i, dx) \\
&\quad + \sum_{i \in E} \iint_{\mathbb{R}_+ \times \mathbb{R}_+} \frac{\partial \varphi}{\partial x}(i, x, t) \mu(i, dx, dt) \\
&\quad - \sum_{i \in E} \iint_{\mathbb{R}_+ \times \mathbb{R}_+} b(i, x) \varphi(i, x, t) \mu(i, dx, dt) \\
&\quad + \sum_{i \in E} \sum_{j \in E} \iint_{\mathbb{R}_+ \times \mathbb{R}_+} a(i, j, x) \varphi(j, 0, t) \mu(i, dx, dt),
\end{aligned}
$$

and this is the requested result.

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

**Theorem 2.7.** Let E be a finite set, and let $a \in C^0(E \times E \times \mathbb{R}_+)$. Let $\rho_{\text{ini}}$ be a probability measure on the set $E \times \mathbb{R}_+$, such that there exists $R_0 > 0$ with $\sum_{i \in E} \int_{[0,R_0]} \rho_{\text{ini}}(i, dx) = 1$. Let $h > 0$ be a strictly positive real number. Then, for all $t \in \mathbb{R}_+$, the probability measure $\bar{u}_t^h(i, x)$ $dx$ converges to $\rho_t(i, dx)$, the unique solution of (12) as $h$ tends to 0.

*Proof.* The above proposition is a direct consequence of Theorem 1.4, Proposition 2.6 and Theorem 2.4.

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

### 3. NUMERICAL APPLICATIONS

We now present some classical examples in reliability theory, in the framework of Remark 1.1: we recall that $E = \{0,1\}$ and that $q$ and $b$ are linked by relations (8) and (9). The mean value of the duration in the state $i$, for $i = 0,1$, is then given by $m_i = \int_{\mathbb{R}_+} tq(i, 1-i,t)dt$, the standard deviation is given by $\sigma_i^2 = \int_{\mathbb{R}_+} t^2q(i, 1-i,t)dt - m_i^2$ and the coefficient of variation $C_i^v$ is given by $C_i^v = \sigma_i/m_i$ (see for example [3]). We thus approximate the marginal distributions $\rho_t(i, dx)$ using the numerical scheme which is described in this paper. The availability $A_1(t)$ is then defined by $A_1(t) = \mathbb{P}(\eta_t = 1) = \int_{\mathbb{R}_+} \rho_t(1, dx)$ and the unavailability is defined by $A_0(t) = \int_{\mathbb{R}_+} \rho_t(0, dx) = 1 - A_1(t)$. It is possible to obtain an indicator of the error committed on $A_0(t)$, using the fact that the asymptotic unavailability is given by

$$ \lim_{t \to \infty} A_0(t) = \frac{m_0}{m_1 + m_0}. \qquad (47) $$

We examine three numerical examples of the calculation of the unavailability $A_0(t) = \mathbb{P}(\eta_t = 0)$ of the component as a function of $t$. In the following table, we precise for $i = 0,1$ the laws of duration $q(i, 1-i, \cdot)$ kept for both states of the three examples. We also give the mean values $m_i$ and the coefficients of variation $C_i^v$.

<table><thead><tr><th></th><th>q(1, 0, ·)</th><th>m<sub>1</sub></th><th>C<sub>1</sub><sup>v</sup> = σ<sub>1</sub>/m<sub>1</sub></th><th>q(0, 1, ·)</th><th>m<sub>0</sub></th><th>C<sub>0</sub><sup>v</sup> = σ<sub>0</sub>/m<sub>0</sub></th></tr></thead><tbody><tr><td>Ex.1</td><td>gamma</td><td>10.</td><td>0.32</td><td>log-normal</td><td>1.65</td><td>1.31</td></tr><tr><td>Ex.2</td><td>log-normal</td><td>4.48</td><td>1.31</td><td>Weibull</td><td>0.9</td><td>0.28</td></tr><tr><td>Ex.3</td><td>Weibull</td><td>886</td><td>0.52</td><td>Weibull</td><td>0.9</td><td>0.28</td></tr></tbody></table>

The Weibull distribution, with parameters $\alpha$ and $\beta$, is given by

$$ q(t) = b(t) \exp \left( - \int_{0}^{t} b(s) ds \right), \quad b(t) = \frac{\beta}{\alpha^{\beta}} t^{\beta-1}, \quad m = \alpha \Gamma \left( 1 + \frac{1}{\beta} \right), \quad \sigma^2 = \alpha^2 \left( \Gamma \left( 1 + \frac{2}{\beta} \right) - \Gamma^2 \left( 1 + \frac{1}{\beta} \right) \right), $$
<figure><img src="image_6.png" /><figcaption>FIGURE 1. Example 1. Left: numerical value of the unavailability $A_0(t)$ as a function of the time, right: approximation of $\bar{u}_{80/3}(i, x)$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 2. Numerical values of the unavailability $A_0(t)$ as a function of the time, Example 2.</figcaption></figure>

the log-normal distribution with parameters $\alpha$ and $\beta$, is given by

$$q(t) = \frac{1}{\sqrt{2\pi}\beta t} \exp\left(\frac{-(\log t - \alpha)^2}{2\beta^2}\right), \quad b(t) = \frac{q(t)}{\int_t^{+\infty} q(s) \, ds}, \quad m = e^{\alpha+\beta^2/2}, \quad \sigma^2 = e^{2\alpha}(e^{2\beta^2} - e^{\beta^2}),$$

and the gamma distribution, with parameters $\alpha$ and $\beta$, is given by

$$q(t) = \frac{1}{\Gamma(\alpha) \beta^{\alpha}} t^{\alpha-1} e^{-t/\beta}, \quad b(t) = \frac{q(t)}{\int_t^{+\infty} q(s) \, ds}, \quad m = \alpha\beta, \quad \sigma^2 = \alpha\beta^2.$$

In Figure 1, the left picture shows, in the case of Example 1, the unavailability as a function of the time, the horizontal dashed line representing the asymptotic value. The right picture represents the curves $\bar{u}_{80/3}(i, x)$ for $i = 0$ and $i = 1$, as a function of $x$, where $\bar{u}_{80/3}(i, x) \, dx$ is an approximation of $\rho_{80/3}(i, \, dx)$ for $t = 80/3$. The peak, precisely obtained for $x = 80/3$, corresponds to the integration in the control volume of a Dirac mass at this point, whose weight is the probability that the component is working during the time interval $[0, t]$. In this case, we have taken $h = 0.01$.

Figure 2 represents the unavailability as a function of the time for Example 2, and the horizontal dashed line again represents the asymptotic value. The precision of the calculation can be evaluated by the distance between the numerical results and the asymptotic value. The results are accurate enough, with $h = 0.01$, which leads to 4000 control volumes.

In the case of Example 3, the results are not precise enough, using Scheme (23)-(25) with $h = 1$. and 2000 control volumes (left picture of Fig. 3). Note that other results of the literature show much less accuracy in this case: this lack of accuracy is due to the high ratio between the mean durations of a working period and of a
<figure><img src="image_8.png" /><figcaption>FIGURE 3. Numerical values of the unavailability $A_0(t)$ as a function of the time (left: Example 3 with Scheme (23)–(25) and $h = 1$; right: Example 3 with Scheme (26)–(28) and $h = 0.01$, $N_1 = 100$, $N_0 = 1$.)</figcaption></figure>

repair period, although these values are usual for reliable materials. Using Scheme (26)–(28), we are then able to obtain much more accurate results (right picture of Fig. 3) with $h = 0.01$ and again 2000 control volumes, using the values $N_1 = 100$ and $N_0 = 1$ (a number of 200 000 control volumes would have been needed for the same calculation with Scheme (23)–(25) and $h = 0.01$).

Other numerical experiments can be found in [4,9].

REFERENCES

[1] S. Asmussen, Fitting phase-type distributions via the EM algorithm. Scand. J. Statist. **23** (1996) 419–441.

[2] E. Cinlar, *Introduction to Stochastic processes*. Prentice-Hall (1975).

[3] C. Cocozza-Thivent, *Processus stochastiques et fiabilité des systèmes*. Springer, Paris, Collection Mathématiques & Applications **28** (1997).

[4] C. Cocozza-Thivent and R. Eymard, Marginal distributions of a semi-Markov process and their computations, *Ninth ISSAT International Conference on Reliability and Quality in Design*, International Society of Science and Applied Technologies, H. Pham and S. Yamada Eds. (2003).

[5] C. Cocozza-Thivent and M. Roussignol, Semi-Markov process for reliability studies. *ESAIM: PS 1* (1997) 207–223.

[6] C. Cocozza-Thivent and M. Roussignol, A general framework for some asymptotic reliability formulas. *Adv. Appl. Prob.* **32** (2000) 446–467.

[7] C. Cocozza-Thivent, R. Eymard, S. Mercier and M. Roussignol, On the marginal distributions of Markov processes used in dynamic reliability, *Prépublications du Laboratoire d'Analyse et de Mathématiques Appliquées UMR CNRS 8050*, 2/2003 (January 2003).

[8] C. Cocozza-Thivent, R. Eymard and S. Mercier, A numerical scheme to solve integro-differential equations in the dynamic reliability field, *PSAM7-ESREL'04*, Berlin (June 2004).

[9] C. Cocozza-Thivent, R. Eymard and S. Mercier, Méthodologie et algorithmes pour la quantification de petits systèmes redondants, *Proceedings of the Conference λ/µ 14*, Bourges, France (October 2004).

[10] D.R. Cox, *Renewal Theory*. Chapman and Hall, London (1982).

[11] R. Eymard, T. Gallouët and R. Herbin, Finite Volume Methods, *Handbook of Numerical Analysis*, P.G. Ciarlet and J.L. Lions Eds., **VII** (2000) 723–1020.

[12] W. Feller, *An Introduction to Probability Theory and its Applications*. Volume II, Wiley (1966).

[13] A. Fritz, P. Pozsgai and B. Bertsche, Notes on the Analytic Description and Numerical Calculation of the Time Dependent Availability, *MMR'2000: Second International Conference on Mathematical Methods in Reliability*, Bordeaux, France, July 4–7 (2000) 413–416.

[14] S. Mischler, B. Perthame and L. Ryzhik, Stability in a nonlinear population maturation model. *Math. Models Met. App. Sci.* **12** (2002) 1–22.