C Cocozza-Thivent

R Eymard

S Mercier
email: sophie.mercier@univ-mlv.fr

A finite-volume scheme for dynamic reliability models

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

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 R d . We finally show on a numerical example the accuracy and efficiency of the approximation method.

Introduction

We consider a system which is described, at each time t ∈ R + , by its physical state I t ∈ E, where E is a finite set, and by environmental variables X t ∈ 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

dX t dt = v(i, X t ) ∀ t ∈ R + s.t., I t = i ∀ i ∈ E, (1.1)
where v is a mapping from E × R d to R d . Stochastic transitions from state (i, x) ∈ E × R d to ( j, y) ∈ E × R d are defined, thanks to a rate a(i, j, x) and a probability measure µ(i, j, x)(dy). We denote by ρ(t)(i, dx) the probability distribution of (I t , X t ), for all t ∈ R + . For such processes (I t , X t ) t∈R + (called 'piecewise deterministic Markov processes' by [START_REF] Davis | Piecewise-deterministic Markov processes: a general class of non-diffusion stochastic models[END_REF][START_REF] Davis | Markov Models and Optimization[END_REF]), the family of marginal distributions (ρ(t)(i, dx)) t∈R + is solution of the Chapman-Kolmogorov equations

t 0 i∈E R d j∈E a(i, j, x) R d
ϕ( j, y)µ(i, j, x)(dy) -ϕ(i, x) ρ(s)(i, dx)ds

+ t 0 i∈E R d v(i, x) • ∇ϕ(i, x)ρ(s)(i, dx)ds - i∈E R d ϕ(i, x)ρ(t)(i, dx) + i∈E R d ϕ(i, x)ρ 0 (i, dx) = 0 ∀ i ∈ E, ∀ t ∈ R + , ∀ ϕ ∈ C 1 c (R d ) E , (1.2)
• the function v(i, •) is sub-linear, i.e. there are some V 1 > 0 and V 2 > 0 such that

∀ i ∈ E, ∀ x ∈ R d , |v(i, x)| V 1 |x| + V 2 .

The function µ is such that

• for all i, j ∈ E and x ∈ R d , µ(i, j, x) ∈ P(R d ), the set of probability measures on R d ,

• for all i, j ∈ E and all ψ ∈ C b (R d , R), the function x → ψ(y)µ(i, j, x)(dy) is continuous from R d to R.

5. We assume that ρ 0 ∈ P(E × R d ), where P(E × R d ) is the set of probability measures on E × R d .

In Cocozza-Thivent et al., we have proved that under (H), there exists one and only one solution ρ(•) to equation (1.2) , namely, a single function ρ(•): R + → P(E × R d ) such that the function t → i∈E R d ψ(i, x)ρ(t) (i, dx) is Lebesgue-integrable and bounded on [0, T ], for all T ∈ R + and for all ψ ∈ C b (R d ) E and such that (1.2) is fulfilled. Moreover, the unique solution ρ(t) is actually such that the function t

→ i∈E R d ψ(i, x) ρ(t)(i, dx) is continuous from R + to R for all ψ ∈ C b (R d ) E .
REMARK 1.1 Hypothesis (H-4) includes the important practical cases µ(i, j, x)(dy) = δ 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 µ(i, j, x)(dy) = δ 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 ∈ L 1 (R d ) E and u ∈ C(R + , L 1 (R d )) E such that ρ 0 (i, dx) = u 0 (i, x)dx and ρ(t)(i, dx) = u(i, x, t)dx, and assuming that µ(i, j, x)(dy) = M(i, j, x, y)dy for all i, j ∈ E and x, y ∈ 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 ∂ t u (i, x, t) + div(u(i, x, t) ( j, i, y)u( j, y, t) M( j, i, y, x) dyu(i, x, t)b(i, x), for a.e. x ∈ R d , ∀ t ∈ R + , ∀ i ∈ E,

v(i, x)) = j∈E R d a
(1.3) u(i, x, 0) = u 0 (i, x), for a.e. x ∈ R d , ∀ i ∈ E.

Let us now take the more particular case, where d = 1; there exist regular functions u 0 ∈ L 1 (R + ) E and u ∈ C(R + , L 1 (R + )) E such that ρ 0 (i, dx) = u 0 (i, x)dx and ρ(t)(i, dx) = u(i, x, t)dx, assuming that µ(i, j, x)(dy) = δ 0 (dy) for all i, j ∈ E and x ∈ R + (see Remark 1.1) and v(i, •) 0 (in order to ensure the existence of a solution, the support of which remains a subset of R + ). Then the term involving µ in (1.2) can be seen as resulting from the weak formulation of a boundary condition on {0} × R + (see [START_REF] Cocozza-Thivent | Approximation of the marginal distributions of a semi-Markov process using a finite volume scheme[END_REF]. 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} × R + :

∂ t u(i, x, t) + div(u(i, x, t)v(i, x)) = -u(i, x, t)b(i, x), for a.e.

x ∈ R + , ∀ t ∈ R + , ∀ i ∈ E, v(i, 0) u(i, 0, t) = j∈E R d a( j, i, y)u( j, y, t)dy, ∀ t ∈ R + , ∀ i ∈ E, (1.4) u(i, x, 0) = u 0 (i, x), for a.e. x ∈ R d , ∀ i ∈ E.
Hence, in the particular case

E = {1}, System (1.4) is completely identical, dropping the condition b(1, •) = a(1, 1, •)
, to a population dynamics model called the McKendrick-Von Foerster model, improved by [START_REF] Sinko | A new model for age-size structure of a population[END_REF] and [START_REF] Bell | Cell growth and division I. A mathematical model with applications to cell volume distributions in mammalian suspension cultures[END_REF]. In this model, u(1, x, t) represents the density of population of age x at the time t, 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, andreferences therein, andsee 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 ρ(t)(•, dx) of the process (I t , X t ). Actually, as can be seen in the literature (see, e.g. [START_REF] Labeau | A Monte Carlo estimation of the marginal distributions in a problem of probabilistic dynamics[END_REF], 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. [START_REF] Labeau | A Monte Carlo estimation of the marginal distributions in a problem of probabilistic dynamics[END_REF], and other papers by P. E. Labeau).

We construct here an approximation of ρ(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 [START_REF] Cocozza-Thivent | Approximation of the marginal distributions of a semi-Markov process using a finite volume scheme[END_REF] for the 'semi-Markov' case). Note that similar numerical approximations for (1.4) have been studied in the 1D case (see [START_REF] Ackleh | A finite difference approximation for a coupled system of nonlinear size-structured populations[END_REF], Ackleh & Ito (1997) and [START_REF] Abia | Size-structured population dynamics models and their numerical solutions[END_REF], 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. [START_REF] Eymard | Finite Volume Methods, Handbook of Numerical Analysis[END_REF]. 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 × R d ), and not in

(L 1 (R d ) ∩ L ∞ (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 ∞ (R d ) must be replaced by that of a family of probability measures in P(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 µ(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 δt → 0 and |M|/δt → 0, where |M| is the space step and δt is the time step. Such a condition instead of a simple bound on |M|/δ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 ∞ framework, where the C F L condition for the convergence of the explicit scheme is expressed as a bound on δ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).

The numerical scheme

We first denote by g:

E × R d × R → R d the solution of ∂g ∂t (i, x, t) = v(i, g(i, x, t)) ∀ (i, x, t) ∈ E × R d × R, (2.1) with g(i, x, 0) = x ∀ (i, x) ∈ E × R d . (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 ∈ C 1 (R; C(R d )) E , where C 1 (R; C(R d ))
is the set of continuously differentiable functions with respect to time, valued in C(R d ), 2. the function x → g(i, x, t) is locally Lipschitz continuous with respect to x on R d , for all t > 0, all i ∈ E, 3. the following property holds:

g(i, g(i, x, s), t) = g(i, x, s + t) ∀ (i, x, s, t) ∈ E × R d × R × R.
(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 ∀ (i, y) ∈ E × R d and ∂ ∂s J (i, y, s) = div v(i, g(i, y, s)) J (i, y, s) ∀ (i, y, s) ∈ E × R d × R,
which leads to

J (i, y, s) = exp s 0 div v(i, g(i, y, τ ))dτ ∀ (i, y, s) ∈ E × R d × R.
We then introduce the notion of admissible mesh M of R d , i.e. a family of measurable subsets of

R d such that 1. K ∈M K = R d , 2. ∀ (K , L) ∈ M 2 , K = L ⇒ K ∩ L = ∅ (M is a partition of R d ), 3. m K := K dx > 0 for all K ∈ M. 4. sup K ∈M diam(K ) < +∞ where diam(K ) = sup {(x,y)∈K 2 } |x -y|. We then set |M| := sup K ∈M diam(K ).
Let a real number δt > 0 and an admissible mesh 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 (2.4) and in order to determine the discrete rates of state transition, we set m L a

(i,K ) 0 by m K p (i,K ) 0 = K ρ 0 (i, dx) ∀ K ∈ M ∀ i ∈ E,
(i j) K L = L a( j, i, y) K µ( j, i, y)(dx)dy ∀ (K , L) ∈ M 2 , ∀ (i, j) ∈ E 2 , (2.5) m K b (i) K = K b(i, x)dx ∀ K ∈ M, ∀ i ∈ E.
We define the following coefficients for the convection part of the scheme, using the notations given at the beginning of this section:

m (i) K L = {x∈K /g(i,x,-δt)∈L} J (i, x, -δt)dx = {y∈L/g(i,y,δt)∈K } dy, (2.6)
which thus satisfy

K ∈M m (i) K L = m L ∀ L ∈ M.
(2.7)

For n ∈ N fixed, assuming p (i,K ) n to be constructed (all i, K ), let us construct p

(i,K )
n+1 . 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 (i,K ) n ∀ x ∈ K , ∀ K ∈ M, ∀ i ∈ E, (2.8)
we consider the weak solution u of the partial differential equation

∂ t u(x, t) + div( u(x, t)v(i, x)) = 0, for a.e. (x, t) ∈ R d × (0, δt), u(x, 0) = P n (i, x), for a.e. x ∈ R d .
This function u is then such that u(x, t) = J (i, x, -t)P n (i, g(i, x, -t)), for a.e. x ∈ R d and all t ∈ (0, δt). We then define m K p (i,K ) n := K u(x, δt)dx, which thus satisfies, thanks to the above definitions,

m K p (i,K ) n = L∈M m (i) K L p (i,L) n ∀ K ∈ M, ∀ i ∈ E.
(2.9)

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

m K p (i,K ) n+1 = 1 δtb (i) K + 1 m K p (i,K ) n + δt j∈E L∈M a (i j) K L δtb ( j) L + 1 m L p ( j,L) n ∀ K ∈ M, ∀ i ∈ E. (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 righthand 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 M,δt

t (•, x)dx of the measure ρ(t)(•, dx) such that P M,δt t (•, x) is constant on each {i} × [nδt, (n + 1)δt[ × K : P M,δt t (i, x) := p (i,K ) n = P n (i, x) ∀ i ∈ E, ∀ x ∈ K , ∀ t ∈ [nδt, (n + 1)δt[, ∀ K ∈ M, ∀ n ∈ N.
(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

(i,K ) n
) K ∈M has to be finite at each time step n (n ∈ N). This can be ensured e.g. by the following assumptions: ρ 0 has compact support and inf K ∈M m K > 0, which implies that the cardinality of {K ∈ M: K ⊂ B(0, R)} is finite for all R > 0, where B(0, R) = {x ∈ 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 δt, although it would be easy to generalise our results to the case of a variable time step.

Convergence study

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| |x| + V 2 V 1 (exp(V 1 t) -1) ∀ x ∈ R d , ∀ t ∈ R + , ∀ i ∈ E. (3.1)
Proof. Let t 0 and x ∈ R d be given. We can write CONSERVATION) Under Hypotheses (H), let M be an admissible mesh on R d in the sense given in Section 2 and let δt > 0. Then there exists one and only one family of non-negative real values ( p

|g(i, x, t) -x| = t 0 v(i, g(i, x, s))ds t 0 (V 1 (|g(i, x, s) -x| + |x|) + V 2 )ds = (V 1 |x| + V 2 )t + V 1 t 0 |g(i, x, s) -x|ds. We set h(t) = |g(i, x, t) -x|. We then have h(t) (V 1 |x| + V 2 )t + V 1 t 0 h(τ )dτ . Defining the continuous function h by h(t) = (V 1 |x| + V 2 )t + V 1 t 0 h(τ )dτ , we then get, using Gronwall's lemma, that h(t) -h(t) 0 for all t 0. Since h(t) = V 1 |x|+V 2 V 1 (exp(V 1 t) -1), we get (3.1). LEMMA 3.2 (MASS
(i,K ) n ) n∈N,i∈E,K ∈M which is the solution of (2.4)-(2.10). Moreover, this family is such that i∈E K ∈M m K p (i,K ) n = 1 ∀ n ∈ N. (3.2)
Proof. The existence and uniqueness of the solution ( p

(i,K ) n ) n∈N,i∈E,K ∈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 i∈E K ∈M m K p

(i,K ) n = 1. Then, P n , p (•,•) n and consequently p (•,•)
n+1 are uniquely defined by (2.8)-(2.10). Moreover, we have

i∈E K ∈M m K p (i,K ) n+1 = i∈E K ∈M 1 δtb (i) K + 1 m K p (i,K ) n + δt j∈E L∈M m L p ( j,L) n δtb ( j) L + 1 i∈E K ∈M a (i j) K L , with i∈E K ∈M a (i j) K L = 1 m L i∈E L a( j, i, y) K ∈M K µ( j, i, y)(dx)dy = b ( j) L .
This leads to

i∈E K ∈M m K p (i,K ) n+1 = i∈E K ∈M 1 δtb (i) K + 1 m K p (i,K ) n + δt j∈E L∈M b ( j) L m L p ( j,L) n δtb ( j) L + 1 = i∈E K ∈M m K p (i,K ) n = i∈E K ∈M L∈M m (i) K L p (i,L) n = i∈E L∈M m L p (i,L) n , (3.3)
which completes the proof.

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 ∈ E, the family {µ(i, j, x)(dy): |x| R} is weakly compact, which is a direct consequence of (H -4). The family {µ(i, j, x)(dy)} x∈B(0,R) is then tight (see, e.g., [START_REF] Billingsley | Convergence of Probability Measures[END_REF]), which allows us to introduce the function f i, j : R + × R + → R + , given by

∀ R ∈ R + , ∀ ε ∈ R + , f i, j (R, ε) = inf r R, sup x∈B(0,R) R d \B(0,r ) µ(i, j, x)(dy) ε . (3.4)
We get the following result.

LEMMA 3.3 (PROPAGATION OF THE PROBABILITY MASS) Under Hypotheses (H), let ε > 0, T > 0, θ > 0 and R 0 > 0 be given. Let the sequences (R (k) ) k∈N and (R (k) ) k∈N be defined by

R (0) = R 0 , (3.5) R (k) = exp(V 1 T )R (k) + V 2 + θ exp(V 1 T ) V 1 (exp(V 1 T ) -1) ∀ k ∈ N, (3.6) R (k+1) = sup i, j∈E f i, j (R (k) + θ T, ε) ∀ k ∈ N. (3.7) For all δt ∈ (0, T ), defining N ∈ N such that N δt T < (N + 1)δt, let the sequence (R (k) n ) n=0,...,N be defined for all k ∈ N by R (k) 0 = R (k) , R (k) n+1 = (R (k) n + θδt) exp(V 1 δt) + V 2 V 1 (exp(V 1 δt) -1) ∀ n = 0, . . . , N -1. (3.8)
Then the following properties hold:

R (k) n = exp(V 1 nδt)R (k) + V 2 V 1 + θδt exp(V 1 δt) exp(V 1 δt) -1 (exp(V 1 nδt) -1) ∀ n = 0, . . . , N , ∀ k ∈ N, (3.9) R (k) n < R (k) n+1 , R (k) R (k) n R (k) < R (k+1) R (k+1) n ∀ n = 0, . . . , N , ∀ k ∈ N, (3.10) sup x∈B(0,R (k-1) n+1 +θδt), i, j∈E R d \B(0,R (k) n+1 ) µ(i, j, x)(dy) ε ∀ n = 0, . . . , N -1, ∀ k ∈ N , (3.11) ∀ k ∈ N, ∀ n = 0, . . . , N -1 ∀ i ∈ E, ∀ x ∈ R d , |x| < R (k) n + θδt ⇒ |g(i, x, δt)| < R (k) n+1 .
(3.12)

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 (k) . Then, using that f i, j (R, ε) is nondecreasing withrespect to R, the relations (3.10) and (3.7) yield

(V 1 δt) -1 V 1 δt and therefore θδt exp(V 1 δt)-1 θ V 1 and of f i, j (R (k) + θ T, ε) R (k) + θ T > R
f i, j (R (k-1) n+1 + θδt, ε) f i, j ( R(k-1) + θ T, ε) R (k) < R (k) n+1 ∀ n = 0, . . . , N -1, ∀ k ∈ N , ∀ i, j ∈ E
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 δt > 0 be given and let M be an admissible mesh on R d in the sense given in Section 2. Let θ > |M|/δt and T > δt.

Let ( p (i,K ) n ) i∈E,n∈N,K ∈M be the solution of (2.4)-(2.10).
Then, for all ε 0 > 0, there exists R > 0 only depending on ε 0 , T , V 1 , V 2 , A, CardE, ρ 0 , µ and θ, such that

K ∈M K ∩B(0,R) =∅ i∈E m K p (i,K ) n 1 -ε 0 ∀ n ∈ N s.t. nδt T, (3.13)
whence the tightness of the family of probability measures

(P M,δt t (•, x)dx) {(t,M,δt):|M|/δt<θ ;0<t,δt<T } on E × R d .
Proof. Let ε > 0 be given (this real number will be chosen as a function of ε 0 at the end of the proof). We denote by R 0 > 0 a real number such that i∈E R d \B(0,R 0 ) ρ 0 (i, dx) ε.

(3.14)

Hence, R 0 only depends on ε and ρ 0 . We then denote by (R (k) ) k∈N and (R (k) ) k∈N the sequences defined for these values ε > 0 and R 0 > 0 by (3.5), (3.6) and (3.7), which therefore only depend on ε, ρ 0 , T , V 1 , V 2 , µ and θ. We finally denote by N ∈ N the integer, such that N δt T < (N + 1)δt, and by (R (k) n ) n=0,...,N the sequence defined for all k ∈ N by (3.8). For all n = 0, . . . , N and k ∈ N, we denote by

M (k) n = {K ∈ M, K ∩ B(0, R (k)
n ) = ∅}, and we define

q (k) n = K ∈M (k) n i∈E m K p (i,K ) n ∀ n = 0, . . . , N , ∀ k ∈ N.
which leads, using (3.14), to q

(k) 0 ε ∀ k ∈ N. (3.15)
The principle of this proof is to find some k 0 ∈ N and some C > 0, independent of M, δt and ε, such that q

(k 0 ) n
Cε, for all n ∈ N such that nδt T , which will be sufficient to conclude, choosing ε = ε 0 /C and R R (k 0 ) for all n ∈ N such that nδt T . Let k ∈ N and n = 0, . . . , N -1 be given.

We have

q (k) n+1 = K ∈M (k) n+1 i∈E ⎛ ⎝ 1 δtb (i) K + 1 m K p (i,K ) n + δt j∈E L∈M a (i j) K L δtb ( j) L + 1 m L p ( j,L) n ⎞ ⎠ .
We thus get

q (k) n+1 = K ∈M (k) n+1 i∈E 1 δtb (i) K + 1 m K p (i,K ) n + δt L∈M (k-1) n+1 j∈E K ∈M (k) n+1 i∈E a (i j) K L δtb ( j) L + 1 m L p ( j,L) n + δt L∈M\M (k-1) n+1 j∈E K ∈M (k) n+1 i∈E a (i j) K L δtb ( j) L + 1 m L p ( j,L) n . (3.16)
We first look for a bound for the third term of the right-hand side of (3.16). Let j ∈ E and L ∈ M\M

(k-1)
n+1 be given. We have

i∈E K ∈M (k) n+1 a (i j) K L = i∈E 1 m L L a( j, i, y) K ∈M (k) n+1 K µ( j, i, y)(dx)dy. (3.17) Since we have L ∩ B(0, R (k-1) n+1 ) = ∅ and diam(L) |M| < θδt, thus we have y ∈ B(0, R (k-1) n+1 + θδt). Since K ⊂ R d \B(0, R (k) n+1 ) we can then write, using (3.11), ∀ y ∈ L K ∈M (k) n+1 K µ( j, i, y)(dx) R d \B(0,R (k) n+1 ) µ( j, i, y)(dx) ε.
(3.18) Thus, we get from (3.17) and (3.18) (recall that B = A CardE)

∀ n = 0, . . . , N -1, ∀ k ∈ N , ∀ j ∈ E, ∀ L ∈ M \ M (k-1) n+1 , i∈E K ∈M (k) n+1 a (i j) K L Bε. (3.19)
Since we get a bound of the second term at the right-hand side of (3.16), using

K ∈M (k) n+1 i∈E a (i j) K L b ( j) L ∀ j ∈ E, ∀ L ∈ M,
we thus get from (3.16) and (3.19)

q (k) n+1 K ∈M (k) n+1 i∈E 1 δtb (i) K + 1 m K p (i,K ) n + δt L∈M (k-1) n+1 j∈E b ( j) L δtb ( j) L + 1 m L p ( j,L) n + δt Bε, using L∈M j∈E m L p ( j,L) n = 1 (see (3.3)), which yields q (k) n+1 K ∈M (k) n+1 i∈E m K p (i,K ) n + δt L∈M (k-1) n+1 \M (k) n+1 j∈E b ( j) L δtb ( j) L + 1 m L p ( j,L) n + δt Bε.
Therefore, since b

q (k) n+1 K ∈M (k) n+1 i∈E m K p (i,K ) n + δt B δt B + 1 L∈M (k-1) n+1 \M (k) n+1 j∈E m L p n ( j, L) + δt Bε. Since R (k-1) n+1 R (k) n+1 implies M (k) n+1 ⊂ M (k-1) n+1 , we get q (k) n+1 1 δt B + 1 K ∈M (k) n+1 i∈E m K p (i,K ) n + δt B δt B + 1 L∈M (k-1) n+1 j∈E m L p ( j,L) n + δt Bε.
(3.20)

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

∀ k ∈ N, ∀ n = 0, . . . , N -1, ∀ i ∈ E, ∀ K ∈ M (k) n+1 , ∀ L ∈ M, m (i) K L = 0 ⇒ L ∈ M (k) n .
This implies using (2.9) and (2.7)

i∈E K ∈M (k) n+1 m K p (i,K ) n q (k) n ∀ n = 0, . . . , N , ∀ k ∈ N.
Hence, we get from (3.20)

q (k) n+1 1 δt B + 1 q (k) n + δt B δt B + 1 q (k-1) n + δt Bε. Since q (k) 0
ε for all k ∈ N and q (0) n 1 for all n ∈ N, we then easily prove by induction on n and on k that, for all k, n ∈ N, q

(k) n ε + q(k) n + nδt Bε, where ( q(k) n ) n,k∈N is defined by ( q(0) n = 1 ∀ n ∈ N) and ( q(k) 0 = 0 ∀ k ∈ N ) and ( q(k) n+1 = q(k) n + δt B q(k-1) n ∀ n = 0, . . . , N , ∀ k ∈ N ).
Then an easy inductive argument shows that, for all n = 0, . . . , N , we have

q(k) n = n! k!(n-k)! (δt B
) k for all k = 0, . . . , n (with the usual convention 0! = 1) and q(k) n = 0 for all k > n. This gives

q(k) n n k k! (δt B) k (BT ) k k! ∀ n ∈ N s.t. nδt T, ∀ k ∈ N.
We thus obtain

q (k) n ε + (BT ) k k! + BεT ∀ n ∈ N s.t. nδt T, ∀ k ∈ N.
We then choose k 0 ∈ N such that (BT ) k 0 k 0 ! ε. Thus, k 0 only depends on ε, B (with B = A CardE) and T . We then get

q (k 0 ) n ε(2 + BT ) ∀ n ∈ N s.t. nδt T.
The conclusion of the above lemma follows by setting R = R (k 0 ) (defined by (3.6)) and ε = ε 0 /(2 + BT ).

Convergence lemmas

Let f ∈ C b (R d ) E be a given function, and let t 1 , t 2 ∈ R + be such that 0 t 1 t 2 . For a given admissible mesh in the sense given in Section 2, and for a given δt > 0, using all the notations of Section 2, we define the expressions T M,δt

1 ( f, t 1 , t 2 ) and T M,δt 2 ( f, t 1 , t 2 ) by T M,δt 1 ( f, t 1 , t 2 ) = δt i∈E n 2 -1 n=n 1 K ∈M m K f (i,K ) p (i,K ) n δtb (i) K + 1 , (3.21) T M,δt 2 ( f, t 1 , t 2 ) = i∈E [t 1 ,t 2 ] R d f (i, x)P M,δt t (i, x)dx dt, (3.22)
where

f (i,K ) = 1 m K K f (i, x)dx ∀ i ∈ E, ∀ K ∈ M, (3.23)
represents the mean value of f on K and where we denote by n 1 , n 2 ∈ N the integers such that n 1 δt t 1 < (n 1 + 1)δt and n 2 δt t 2 < (n 2 + 1)δt.

REMARK 3.1 Throughout this paper, we use the convention that, if n 1 n 2 , all the sums n 2 -1 n=n 1 (•) are set equal to zero.

Similarly, for ϕ ∈ C 1 c (R d ) E , we define the expressions T M,δt 3 (ϕ, t 1 , t 2 ), T M,δt 4 (ϕ, t 1 , t 2 ), T M,δt 5 (ϕ, t 1 , t 2 ), T M,δt 6 (ϕ, t 1 , t 2 ), T M,δt 7 (ϕ, t 1 , t 2 ) and T M,δt 8 (ϕ, t 1 , t 2 ) by T M,δt 3 (ϕ, t 1 , t 2 ) = i∈E n 2 -1 n=n 1 K ∈M m K p (i,K ) n δtb (i) K δtb (i) K + 1 ϕ (i,K ) , (3.24) T M,δt 4 (ϕ, t 1 , t 2 ) = i∈E t 2 t 1 R d b(i, x)ϕ(i, x)P M,δt t (i, x)dx dt, (3.25) T M,δt 5 (ϕ, t 1 , t 2 ) = i∈E n 2 -1 n=n 1 K ∈M δt j∈E L∈M a (i j) K L m L p ( j,L) n δtb ( j) L + 1 ϕ (i,K ) , (3.26) T M,δt 6 (ϕ, t 1 , t 2 ) = i∈E t 2 t 1 R d j∈E a(i, j, x) ϕ( j, y)µ(i, j, x)(dy) P M,δt t (i, x)dx dt, (3.27) T M,δt 7 (ϕ, t 1 , t 2 ) = i∈E n 2 -1 n=n 1 K ∈M [m K p (i,K ) n+1 -m K p (i,K ) n ]ϕ (i,K ) , (3.28) T M,δt 8 (ϕ, t 1 , t 2 ) = i∈E - t 2 t 1 R d v(i, x) • ∇ϕ(i, x)P M,δt t (i, x)dx dt - R d ϕ(i, x)P M,δt t 1 (i, x)dx + R d ϕ(i, x)P M,δt t 2 (i, x)dx . (3.29)
Note that considering Scheme (2.10) multiplied by ϕ (i,K ) and summed over i ∈ E, K ∈ M and n = n 1 , . . . , n 2 -1, we exactly get

-T M,δt 3 (ϕ, t 1 , t 2 ) + T M,δt 5 (ϕ, t 1 , t 2 ) -T M,δt 7 (ϕ, t 1 , t 2 ) = 0. (3.30)
On the other hand, using the above definitions, we can write

-T M,δt 4 (ϕ, t 1 , t 2 ) + T M,δt 6 (ϕ, t 1 , t 2 ) -T M,δt 8 (ϕ, t 1 , t 2 ) = i∈E ⎡ ⎣ t 2 t 1 R d ⎡ ⎣ j∈E a(i, j, x) ϕ( j, y)µ(i, j, x)(dy) -b(i, x)ϕ(i, x) + v(i, x) • ∇ϕ(i, x) ⎤ ⎦ P M,δt t (i, x)dx dt + R d ϕ(i, x)P M,δt t 1 (i, x)dx - R d ϕ(i, x)P M,δt t 2 (i, x)dx ⎤ ⎦ . (3.31)
The following results will lead to the convergence of the scheme, by proving that the differences between T

> 0, f ∈ C b (R d ) E , T > 0 and t 1 , t 2 ∈ R + be such that 0 t 1 t 2 T.
Then, for each ε 0 > 0 , there exists δt 0 , which only depends on ε 0 , T , f , V 1 , V 2 , A, CardE, ρ 0 , µ and θ , such that for all δt ∈ [0, δt 0 ], for all admissible meshes M in the sense given in Section 2, such that |M| θδt, with P M,δt given by (2.4-2.11) and T M,δt 1 ( f, t 1 , t 2 ) and T M,δt 2 ( f, t 1 , t 2 ) defined, respectively, by (3.21) and (3.22), we have

|T M,δt 1 ( f, t 1 , t 2 ) -T M,δt 2 ( f, t 1 , t 2 )| ε 0 , (3.32)
which means precisely that lim δt→0,|M| θδt (T M,δt

1 ( f, t 1 , t 2 ) -T M,δt 2 ( f, t 1 , t 2 )) = 0.
Proof. Let ε > 0. For a given admissible mesh in the sense given in Section 2, and for a given δt > 0, using all the notations of Section 2, we define the expression T M,δt

9 ( f, t 1 , t 2 ) by T M,δt 9 ( f, t 1 , t 2 ) = δt i∈E n 2 -1 n=n 1 K ∈M m K f (i,K ) p (i,K ) n .
We then have

|T M,δt 1 ( f, t 1 , t 2 ) -T M,δt 9 ( f, t 1 , t 2 )| δt Bδt i∈E n 2 -1 n=n 1 K ∈M m K | f (i,K ) | p (i,K ) n δt BT sup i∈E f (i, •) L ∞ (R d ) ,
and therefore there exists δt 1 > 0 only depending on ε, f , B and T such that if δt δt 1 , then

|T M,δt 1 ( f, t 1 , t 2 ) -T M,δt 9 ( f, t 1 , t 2 )| ε.
13

We now consider T M,δt 10 ( f, t 1 , t 2 ) defined by

T M,δt 10 ( f, t 1 , t 2 ) = i∈E n 2 δt n 1 δt R d f (i, x)P M,δt t (i, x)dx dt = δt i∈E n 2 -1 n=n 1 K ∈M m K f (i,K ) p (i,K ) n .
We then get that

|T M,δt 10 ( f, t 1 , t 2 ) -T M,δt 2 ( f, t 1 , t 2 )| i∈E R d (n 1 +1)δt n 1 δt | f (i, x)|P M,δt t (i, x)dt + (n 2 +1)δt n 2 δt | f (i, x)|P M,δt t (i, x)dt dx,
and therefore |T M,δt 10

( f, t 1 , t 2 )-T M,δt 2 ( f, t 1 , t 2 )| 2δt sup i∈E f (i, •) L ∞ (R d ) .
Therefore, there exists δt 2 , only depending on ε and f , such that |T M,δt 10

( f, t 1 , t 2 ) -T M,δt 2 ( f, t 1 , t 2 )| ε. Setting T M,δt 11 ( f, t 1 , t 2 ) = T M,δt 9 ( f, t 1 , t 2 ) -T M,δt 10 ( f, t 1 , t 2 ), we get T M,δt 11 ( f, t 1 , t 2 ) = δt i∈E n 2 -1 n=n 1 K ∈M p (i,K ) n L∈M m (i) L K ( f (i,L) -f (i,K ) ).
We then get that, for a given R > 0, T M,δt

11 ( f, t 1 , t 2 ) = T M,δt,R 12 ( f, t 1 , t 2 ) + T M,δt,R 13 ( f, t 1 , t 2 ), with T M,δt,R 12 ( f, t 1 , t 2 ) = δt i∈E n 2 -1 n=n 1 K ∈M K ∩B(0,R)=∅ p (i,K ) n L∈M m (i) L K ( f (i,L) -f (i,K ) )
and

T M,δt,R 13 ( f, t 1 , t 2 ) = δt i∈E n 2 -1 n=n 1 K ∈M K ∩B(0,R) =∅ p (i,K ) n L∈M m (i) L K ( f (i,L) -f (i,K ) ). Since |T M,δt,R 12 ( f, t 1 , t 2 )( f, δt, M, R)| 2 sup i∈E f (i, •) L ∞ (R d ) δt i∈E n 2 -1 n=n 1 K ∈M K ∩B(0,R)=∅ p (i,K ) n m K ,
thanks to Lemma 3.4, letting δt T and therefore |M| θ T , it is possible to choose an R which only depends on ε,

T , V 1 , V 2 , A, CardE, ρ 0 , µ, f and θ such that |T M,δt,R 12 ( f, t 1 , t 2 )| εT . Now let K ∈ M be such that K ∩ B(0, R) = ∅ and L such that m (i) L K = 0. Thanks to Lemma 3.1, we can choose R R + θ T R + diam(K ), only depending on T , V 1 , V 2 , θ and R such that L ⊂ B(0, R ) and K ⊂ B(0, R ). Since f is continuous, it is uniformly continuous on B(0, R ). Let α ∈ (0, 1) be such that | f (i, x) -f (i, y)| ε for all (i, x, y) ∈ E × B(0, R ) × B(0, R ) with |x -y| α.
Thanks again to Lemma 3.1, let δt 3 > 0, which only depends on α, θ , V 1 , V 2 and T , be such that, for all δt δt 3 , we have |g(i, z, δt) -z| α/3 for all z ∈ B(0, R ).

As m (i)

L K = 0, there is a z ∈ K such that g(i, z, δt) ∈ L. For δt min(δt 3 , α 3θ ), we deduce that, for all x ∈ K and y ∈ L,

|x -y| |x -z| + |z -g(i, z, δt)| + |g(i, z, δt) -y| α and consequently | f (i,L) -f (i,K ) | ε. This proves that |T M,δt,R

13

( f, t 1 , t 2 )| εT . Choosing at the beginning of the proof ε such that ε + ε + εT + εT = ε 0 , and then taking δt 0 = min(δt 1 , δt 2 , T, δt 3 , α 3θ ) completes the proof. We now deduce the following result.

LEMMA 3.6 Under assumptions (H), let θ > 0, ϕ ∈ C 1 c (R d ) E , T > 0 and t 1 , t 2 ∈ R + be such that 0 t 1 t 2 T .
Then, for each ε 0 > 0, there exists δt 0 , which only depends on ε 0 , T , ϕ, V 1 , V 2 , A, CardE, ρ 0 , µ and θ , such that for all δt ∈ [0, δt 0 ], for all admissible meshes M in the sense given in Section 2, such that |M| θδt, with P M,δt be given by (2.4)-(2.11) and T M,δt

3 (ϕ, t 1 , t 2 ), T M,δt 4 (ϕ, t 1 , t 2 ), T M,δt 5 (ϕ,
|T M,δt 3 (ϕ, t 1 , t 2 ) -T M,δt 4 (ϕ, t 1 , t 2 )| ε 0 , (3.33)
and

|T M,δt 5 (ϕ, t 1 , t 2 ) -T M,δt 6 (ϕ, t 1 , t 2 )| ε 0 , (3.34)
which means precisely that lim δt→0,|M| θδt (T M,δt

3 ( f, t 1 , t 2 ) -T M,δt 4 ( f, t 1 , t 2 )) = 0 and lim δt→0,|M| θδt (T M,δt 5 ( f, t 1 , t 2 ) -T M,δt 6 ( f, t 1 , t 2 )) = 0.
Proof. We consider the function f defined by

f (i, x) = b(i, x)ϕ(i, x) ∀ (i, x) ∈ E × R d . Setting C ϕ = sup i∈E ∇ϕ(i, •) L ∞ (R d ) , we then have |b (i) K ϕ (i,K ) -bϕ (i,K ) | C ϕ |M|B, which yields that |T M,δt 3 (ϕ, t 1 , t 2 ) -T M,δt 1 ( f, t 1 , t 2 )| C ϕ |M|BT. Noting that T M,δt 2 ( f, t 1 , t 2 ) = T M,δt 4
(ϕ, 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) = j∈E a(i, j, x) ϕ( j, y)µ(i, j, x)(dy) ∀ (i, x) ∈ E × R d .
Note that, although the function ϕ 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 ∈ C b (R d ) E . We then have

m K f (i,K ) - j∈E L∈M m K a ( ji) L K ϕ ( j,L) = {x∈K } j∈E a(i, j, x) L∈M {y∈L} ϕ( j, y)µ(i, j, x)(dy) dx - j∈E L∈M {x∈K } a(i, j, x) {y∈L} µ(i, j, x)(dy)dx 1 m L {z∈L} ϕ( j, z)dz j∈E L∈M {x∈K } {y∈L} a(i, j, x) ϕ( j, y) - 1 m L {z∈L} ϕ( j, z)dz µ(i, j, x)(dy)dx C ϕ |M|Bm K .
We thus get that

|T M,δt 5 (ϕ, t 1 , t 2 ) -T M,δt 1 ( f, t 1 , t 2 )| C ϕ |M|BT.
Noting again that T M,δt

2 ( f, t 1 , t 2 ) = T M,δt 6
(ϕ, t 1 , t 2 ), the conclusion follows from the application of Lemma 3.5 to f . We now study the difference between T M,δt 7 (ϕ, t 1 , t 2 ) and T M,δt 8 (ϕ, t 1 , t 2 ) defined, respectively, by (3.28) and (3.29).

LEMMA 3.7 Under assumptions (H), let ϕ

∈ C 1 c (R d ) E , T > 0 and t 1 , t 2 ∈ R + be such that 0 t 1 t 2 T.
Then, for each ε 0 > 0, there exist δt 0 and θ 0 which only depends on ε 0 , T , ϕ, v, V 1 , V 2 , A, CardE, ρ 0 , and µ such that for all δt ∈ [0, δt 0 ], for all admissible meshes M in the sense given in Section 2, such that |M| θ 0 δt, with P M,δt given by (2.4)-(2.11) and T M,δt 7 (ϕ, t 1 , t 2 ) and

T M,δt 8 ( f, t 1 , t 2 ) defined by (3.28)-(3.29), we have that |T M,δt 7 (ϕ, t 1 , t 2 ) -T M,δt 8 (ϕ, t 1 , t 2 )| ε 0 , (3.35)
which means precisely that lim δt→0,|M|/δt→0 (T M,δt

7 ( f, t 1 , t 2 ) -T M,δt 8 ( f, t 1 , t 2 )) = 0.
Proof. Let ε > 0 be given. Let M be an admissible mesh in the sense given in Section 2, and let δt > 0 be given. Using the notations from Section 2, we define the function ϕ M by ϕ M (i, x) = ϕ (i,K ) = 1 m K K ϕ(i, y)dy for all x ∈ K and all K ∈ M. We also note that, setting

C ϕ = sup i∈E ∇ϕ(i, •) L ∞ (R d ) , we have |ϕ M (i, x) -ϕ(i, x)| C ϕ |M| ∀x ∈ R d , ∀ i ∈ E. (3.36)
We get from (2.6), (2.8), (2.9), (2.11) and (3.28)

T M,δt 7 (ϕ, t 1 , t 2 ) = i∈E K ∈M n 2 -1 n=n 1 K [P n+1 (i, x)ϕ(i, x) -P n (i, g(i, x, -δt))J (i, x, -δt)ϕ M (i, x)]dx ds = 1 δt i∈E n 2 -1 n=n 1 (n+1)δt nδt R d [P M,δt s+δt (i, x)ϕ(i, x) -P M,δt s (i, g(i, x, -δt))J (i, x, -δt)ϕ M (i, x)]dx ds = 1 δt i∈E R d (n 2 +1)δt (n 1 +1)δt P M,δt t (i, y)ϕ(i, y)dt - n 2 δt n 1 δt P M,δt s (i, y)ϕ M (i, g(i, y, δt))ds dy.
This leads to T M,δt

7 (ϕ, t 1 , t 2 ) = -T M,δt 14 (ϕ, t 1 , t 2 ) -T M,δt 15 (ϕ, t 1 ) + T M,δt 16
(ϕ, t 2 ), where we define

T M,δt 14 (ϕ, t 1 , t 2 ) = 1 δt i∈E n 2 δt (n 1 +1)δt R d P M,δt s (i, x)[ϕ M (i, g(i, x, δt)) -ϕ(i, x)]dx ds, T M,δt 15 (ϕ, t 1 ) = 1 δt i∈E (n 1 +1)δt n 1 δt R d P M,δt s (i, y)ϕ M (i, g(i, y, δt))dy ds and T M,δt 16 (ϕ, t 2 ) = 1 δt i∈E (n 2 +1)δt n 2 δt R d P M,δt s (i, y)ϕ(i, y)dy ds.
We first remark that

T M,δt 15 (ϕ, t 1 ) = i∈E R d P M,δt t 1 (i, y)ϕ M (i, g(i, y, δt))dy,
and that

T M,δt 16 (ϕ, t 2 ) = i∈E R d P M,δt t 2 (i, y)ϕ(i, y)dy.
Hence, T M,δt 16 (ϕ, t 2 ) is exactly one of the three terms constituting T M,δt 8 (ϕ, t 1 , t 2 ). Applying Lemma 3.1 and using the fact that the support of ϕ is compact gives the existence of some

C 1 > 0 only depending on V 1 , V 2 , T and ϕ such that |ϕ M (i, g(i, y, δt)) -ϕ M (i, y)| (δt + 2|M|)C 1 .
Let θ 1 be given (for instance θ 1 = 1), we then get that there exists δt 1 > 0 such that, for δt δt 1 and |M| θ 1 δt, we have

T M,δt 15 (ϕ, t 1 ) - i∈E R d ϕ(i, x)P M,δt t 1 (i, x)dx ε.
We have T M,δt

14 (ϕ, t 1 , t 2 ) = T M,δt 17 (ϕ, t 1 , t 2 ) + T M,δt 18 (ϕ, t 1 , t 2 )
, where we define

T M,δt 17 (ϕ, t 1 , t 2 ) = 1 δt i∈E n 2 δt (n 1 +1)δt R d P M,δt s (i, x)[ϕ(i, g(i, x, δt)) -ϕ(i, x)]dx ds and T M,δt 18 (ϕ, t 1 , t 2 ) = 1 δt i∈E n 2 δt (n 1 +1)δt R d P M,δt s (i, x)[ϕ M (i, g(i, x, δt)) -ϕ(i, g(i, x, δt))]dx ds. Defining T M,δt 19 (ϕ, t 1 , t 2 ) by T M,δt 19 (ϕ, t 1 , t 2 ) = i∈E t 2 t 1 R d P M,δt s (i, x)v(i, x) • ∇ϕ(i, x)dx ds for all t 1 , t 2 ∈ R + such that 0 t 1 t 2 T , we have T M,δt 19 (ϕ, t 1 , t 2 ) = T M,δt 19 (ϕ, t 1 , (n 1 + 1)δt) + T M,δt 19 (ϕ, (n 1 + 1)δt, n 2 δt) + T M,δt 19
(ϕ, n 2 δt, t 2 ). We note that there exists δt 2 > 0, only depending on ϕ, V 1 , V 2 , T and a bound of v on supp(ϕ), such that, for δt δt 2 , |T M,δt

19 (ϕ, t 1 , (n 1 + 1)δt)| ε and |T M,δt 19 (ϕ, n 2 δt, t 2 )| ε. We can write, since ∂ τ ϕ(i, g(i, x, τ )) = v(i, g(i, x, τ )) • ∇ϕ(i, g(i, x, τ )), ϕ(i, g(i, x, δt)) -ϕ(i, x) = δt 0 v(i, g(i, x, τ )) • ∇ϕ(i, g(i, x, τ ))dτ.
Thanks to the regularity properties of ϕ, this gives the existence of some function h: R + → R with h(x) → 0 as x → 0, which only depends on ϕ, V 1 , V 2 , T and on the properties of v on supp(ϕ), such that ϕ(i,

g(i, x, δt)) -ϕ(i, x) δt -v(i, x) • ∇ϕ(i, x) h(δt) ∀ x ∈ R d .
We thus obtain that there exists δt 3 > 0 such that, for δt δt 3 , we get

|T M,δt 17 (ϕ, t 1 , t 2 ) -T M,δt 19 (ϕ, (n 1 + 1)δt, n 2 δt)| ε.
We now observe that, thanks to (3.36),

|T M,δt 18 (ϕ, t 1 , t 2 )| C ϕ T |M| δt .
We can then choose θ 2 = ε, which gives, for |M| δtθ 2 , |T M,δt 18 (ϕ, t 1 , t 2 )| C ϕ T ε. Gathering the above results, it suffices to choose ε, such that 4ε + C ϕ T ε = ε 0 , δt 0 = min(δt 1 , δt 2 , δt 3 ) and θ 0 = min(θ 1 , θ 2 ), to conclude the proof of the lemma. REMARK 3.2 The hypothesis |M|/δt -→ 0 is sufficient to guarantee the existence of some θ > 0 such that |M| θδt. This hypothesis is only imposed in order to ensure that T M,δt 18 (ϕ, t 1 , t 2 ) -→ 0, and was not needed elsewhere. In some special cases, however, the term T M,δt 18 (ϕ, t 1 , t 2 ) vanishes so that the condition |M|/δt -→ 0 is not required any more. Consider, for example, the case d = 1,

M = {I k = [kδt,(k + 1)δt[ with k ∈ Z} and g(i, x, s) = x + m i s with m i ∈ Z.
In the case where ρ 0 (i, x) = u 0 (i, x)dx, with u 0 ∈ L 2 (R d ) E and µ(i, j, x)(dy) = δ x , it is possible to prove a weak bounded variation inequality, namely, the existence of some C > 0 such that

i∈E N n=0 K ,L∈M m (i) K L ( p (i,L) n -p (i,K ) n+1 ) 2 C, N δt T < (N + 1)δt, and an L 2 (R d × R + ) estimate i∈E K ∈M m K ( p (i,K ) n+1 ) 2 C ∀ n ∈ N s.t. nδt T.
The weak limit is then a function of L 2 (R d × R + ) E , and the convergence can be proved under the hypothesis |M| θδt. We do not give the details of this result here since it does not lead to a convergence proof for more general ρ 0 and µ.

Continuity and compactness

We can now state a lemma which ensures the continuity of the approximate solution with respect to the time variable. Let ϕ ∈ C b (R d ) E be a given function. For a given admissible mesh in the sense given in Section 2, and for a given δt > 0, using all the notations of Section 2, we set

S ϕ (M, δt, s, t) = i∈E R d ϕ(i, x)P M,δt t (i, x)dx - i∈E R d ϕ(i, x)P M,δt s (i, x)dx ∀ s, t ∈ R + .
(3.37)

The statement is then the following.

LEMMA 3.8 Under assumptions (H), let ϕ ∈ C 1 c (R d ) E and let T > 0 be given. Then, there exists C S , which only depends on T , ϕ, V 1 , V 2 , A, CardE, ρ 0 and µ, such that, for each ε 0 > 0, there exist δt 0 and θ 0 which only depends on ε 0 , T , ϕ, v, V 1 , V 2 , A, CardE, ρ 0 , and µ such that for all δt ∈ [0, δt 0 ], for all admissible meshes M in the sense given in Section 2, such that |M| θ 0 δt, with P M,δt given by (2.4)-(2.11) and S ϕ (M, δt, s, t) defined by (3.37), we have that Proof. Let ε 0 > 0 be given. Let θ 0 be given by Lemma 3.7, and let δt 0 be the minimum between the values δt 0 given by Lemmas 3.5, 3.6 and 3.7, setting in θ = θ 0 in Lemmas 3.5 and 3.6 let δt ∈ (0, δt 0 ) and M be an admissible mesh in the sense given in Section 2, such that |M| θ 0 δt, with P M,δt given by (2.4)-(2.11).

|S ϕ (M, δt, s, t)| C S |t -s| + ε 0 , ( 3
Starting from (3.31), we get

i∈E R d ϕ(i, x)P M,δt t 1 (i, x)dx - i∈E R d ϕ(i, x)P M,δt t 2 (i, x)dx | -T M,δt 4 (ϕ, t 1 , t 2 ) + T M,δt 6 (ϕ, t 1 , t 2 ) -T M,δt 8 (ϕ, t 1 , t 2 )| + i∈E t 2 t 1 R d j∈E a(i, j, x) ϕ( j, y)µ(i, j, x)(dy) -b(i, x)ϕ(i, x) + v(i, x) • ∇ϕ(i, x) P M,δt t (i, x)dx dt .
Using (3.30) and then Lemmas 3.6 and 3.7, we have

| -T M,δt 4 (ϕ, t 1 , t 2 ) + T M,δt 6 (ϕ, t 1 , t 2 ) -T M,δt 8 (ϕ, t 1 , t 2 )| = |T M,δt 3 (ϕ, t 1 , t 2 ) -T M,δt 4 (ϕ, t 1 , t 2 ) + T M,δt 6 (ϕ, t 1 , t 2 ) -T M,δt 5 (ϕ, t 1 , t 2 ) -T M,δt 8 (ϕ, t 1 , t 2 ) + T M,δt 7 (ϕ, t 1 , t 2 )| 3ε 0 .
Let V be a bound of v on supp(ϕ). Setting

C S = 2B sup i∈E ϕ(i, •) L ∞ (R d ) + V sup i∈E ∇ϕ(i, •) L ∞ (R d )
provides the result by substituting ε 0 /3 by ε 0 at the beginning of the proof.

We can now conclude with the convergence theorem.

THEOREM 3.1 Under assumptions (H), let (M m , δt m ) m∈N be a sequence such that, for all m ∈ N, M m is an admissible mesh in the sense given in Section 2 and such that δt m → 0 and |M m | δt m → 0 as m → ∞. Then, for all t ∈ R + , the tight sequence of probability measures (P M m ,δt m t (•, x) dx) m∈N weakly converges to a probability measure on E × R d , denoted by ρ(t)(•, dx). Moreover, this family (ρ(t)(•, dx)) t∈R + is the unique weak solution of the Chapman-Kolmogorov equations (1.2).

Proof. Let (t p ) p∈N be a sequence of real numbers, dense in R + and let (M σ p (k) , δt σ p (k) ) k∈N be a sub-sequence of (M m , δt m ) m∈N such that (P

M σ p (k) ,δt σ p (k)
t p (•, x)dx) k∈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 σ (k) , δt σ (k) ) k∈N of (M m , δt m ) m∈N , such that Let us now prove that, for all

P M σ (k) ,δt σ (k)
ψ ∈ C 1 c (R d ) E , the sequence ( i∈E R d ψ(i, x)P M m ,δt m t (i, x)dx) m∈N is a Cauchy sequence for all t > 0.
Let k, m ∈ N. For all p ∈ N, we may write, using the notation (3.37),

i∈E R d ψ(i, x)P M k ,δt k t (i, x)dx - i∈E R d ψ(i, x)P M m ,δt m t (i, x)dx |S ψ (M k , δt k , t, t p )| + |S ψ (M m , δt m , t, t p )| + i∈E R d ψ(i, x)P M k ,δt k t p (i, x)dx - i∈E R d ψ(i, x)P M m ,
, i∈E R d ψ(i, x)P M k ,δt k t p (i, x)dx - i∈E R d ψ(i, x)P M m ,δt m t p (i, x)dx ε.
We can now choose m 1 0 such that, for all m m 1 , δt m δt 0 and |M m | θ 0 δt m , where δt 0 and θ 0 are given by Lemma 3.8 for the value ε. We thus have, for m, k m 1 ,

|S ψ (M k , δt k , t, t p )| + |S ψ (M m , δt m , t, t p )| 2(C S |t -t p | + ε) (2C S + 2)ε.
This completes the proof that ( i∈E R d ψ(i, x)P M m ,δt m t (i, x)dx) m∈N is a Cauchy sequence for all t > 0. We easily obtain the same property for all ψ ∈ C c (R d ) E by density of C 1 c (R d ) E and then for all ψ ∈ C b (R d ) E , thanks to the tightness Lemma 3.5. Therefore, there exists a probability measure ρ(t)(•, dx) such that the sequence (P M m ,δt m t (•, x)dx) m∈N is weakly convergent to ρ(t)(•, dx). Moreover, passing to the limit in (3.38), we get that, for all

ϕ ∈ C 1 c (R d ) E , i∈E R d ϕ(i, x)ρ(t)(i, dx) - i∈E R d ϕ(i, x)ρ(s)(i, dx) C S |s -t| ∀s, t ∈ R + ,
which gives the continuity of the mapping t → i∈E ϕ(i, x)ρ(t)(i, dx). The continuity property is easily extended by density to all ϕ ∈ C c (R d ) E and to all ϕ ∈ C b (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

ϕ ∈ C 1 c (R d ) E and t 1 , t 2 0: lim m→∞ T M m ,δt m 3 (ϕ, t 1 , t 2 ) = lim m→∞ T M m ,δt m 4 (ϕ, t 1 , t 2 ) = i∈E t 2 t 1 R d b(i, x)ϕ(i, x)ρ(t)(i, dx)dt, lim m→∞ T M m ,δt m 5 (ϕ, t 1 , t 2 ) = lim m→∞ T M m ,δt m 6 (ϕ, t 1 , t 2 ) = i∈E t 2 t 1 R d j∈E a(i, j, x) ϕ( j, y)µ(i, j, x)(dy) ρ(t)(i, dx)dt, lim m→∞ T M m ,δt m 7 (ϕ, t 1 , t 2 ) = lim m→∞ T M m ,δt m 8 (ϕ, t 1 , t 2 ) = - i∈E t 2 t 1 R d v(i, x) • ∇ϕ(i, x)ρ(t)(i, dx)dt + R d ϕ(i, x)ρ(t 1 )(i, dx) - R d ϕ(i, x)ρ(t 2 )(i, dx) .
We now pass to the limit in (3.30). Since it is clear that ρ(0)(i, dx) = ρ 0 (i, dx) thanks to (2.4), we deduce that (ρ(t)(i, dx)) t∈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.

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 µ(i, j, x)(dy) and on the initial distribution ρ 0 (•, dx) which allow one to control the expansion of the support of ρ(t)(•, dx) (or alternatively of P M n ,δt m t (i, x)). Indeed, assume for example that the support of ρ 0 (•, dx) is included in some B(0, M 0 )(the closure of B(0, M 0 )) and that for all i, j ∈ E, the support of µ(i, j, x)(dy) is included in B(0, |x|), it is then easy to show that, for all T > 0, all 0 < t < T , the support of ρ(t)(•, dx) is bounded and is included in some B(0, M) where M depends on M 0 , V 1 , V 2 and T . As transition rates are bounded on 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 M m ,δt m t (•, x)dx) m∈N to ρ(t)(•, dx) for all 0 < t < T .

Numerical example

We present the results obtained using the finite-volume scheme (2.4)-(2.10) on a benchmark taken from [START_REF] Labeau | A Monte Carlo estimation of the marginal distributions in a problem of probabilistic dynamics[END_REF]. The states of the system are defined by E = {1, 2, 3}, and the environment variables are described by x ∈ R 2 . The only possible transitions between states are 1 → 2 → 3, which means that a(i, j, x) = 0 only for (i, j) = (1, 2) or (i, j) = (2, 3). We assume that a(1, 2, x) = λ 1 and a(2, 3, x) = λ 2 for all x = (x (1) , x (2) ) ∈ R 2 , where λ 1 = 0.5 and λ 2 = 0.3 and that µ(i, j, x)(dy) = δ 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 = 1.5 0 1 0.5 , A 2 = 1 0 0.75 0.25 and A 3 = 0.75 0 0 0.5 . The initial data ρ 0 is given by the Dirac measure δ x 0 with x 0 = (1, 1) in the state i = 1, which means that i∈E R 2 ϕ(i, x)ρ 0 (i, dx) = ϕ(1, x 0 ), for all ϕ ∈ C 1 b (R 2 , R) E . An analytical solution can be obtained in this case [START_REF] Labeau | A Monte Carlo estimation of the marginal distributions in a problem of probabilistic dynamics[END_REF], which satisfies that ρ(t)(1, •) is a weighted Dirac measure located at the point g(1, x 0 , t) of R 2 , that ρ(t)(2, •) is supported by a line of R 2 and that ρ(t)(3, •) admits a density with respect to the Lebesgue measure, denoted p(t), supported by a triangle.

We consider meshes on R 2 consisting of squares with side |M|. In order to confine the computations as much as possible to the geometrical domain which supports ρ(t), we performed consecutive simulations with the values |M| = 5×10 -3 , |M| = 3.13×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 ρ(0.1)(3, •). We then compare, in Fig. 2, the numerical solution (top) with the analytical one (bottom) for the density p(0.1) of ρ(0.1)(3, •), 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) ) = R p(t)(x (1) , x (2) )dx (2) , for a.e. x (1) ∈ R, and p y (t)(x (2) ) = R p(t)(x (1) , x (2) )dx (1) , for a.e. x (2) ∈ R. We compare in Figs 3 and4, respectively, the analytical and numerical results for p x (t) and p y (t): these figures show an excellent agreement. 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 [START_REF] Labeau | A Monte Carlo estimation of the marginal distributions in a problem of probabilistic dynamics[END_REF]), the computing time is here about three quarters of an hour on a PC, which is longer than those given in [START_REF] Labeau | A Monte Carlo estimation of the marginal distributions in a problem of probabilistic dynamics[END_REF]. The precision in [START_REF] Labeau | A Monte Carlo estimation of the marginal distributions in a problem of probabilistic dynamics[END_REF] 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 [START_REF] Cocozza-Thivent | Marginal distributions of a semi-Markov process and their computations[END_REF]), 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. • 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.

  x)dx k∈N is weakly convergent for all p ∈ N. Such a sub-sequence is again denoted by (M m , δt m ) m∈N in the following and (P M m ,δt m t p (•, x)dx) m∈N is now weakly convergent to a probability distribution ρ(t p )(•, dx) for all p ∈ N.

FIG. 2 .

 2 FIG. 2. Density p(0.1). Top: numerical approximation with |M| = 10 -5 , bottom: analytical solution.

FIG. 3 .

 3 FIG. 3. Analytical and numerical function p x (0.1)(x).

FIG. 4 .

 4 FIG. 4. Analytical and numerical function p y (0.1)(y).

  We could complete the proof of the above lemma assuming only |M| θδt instead of |M|/δt → 0, which would only lead to a different choice for C S .

	.38)
	which means precisely that lim sup δt→0,|M|/δt→0 |S ϕ (M, δt, s, t)| C S |t -s|.
	REMARK 3.3

  (•, x)dx) m∈N is weakly convergent, there exists m 0 such that, for all m, k m 0

	t p	δt m	(i, x)dx	(3.39)
	Let T = t + 1 and let ε ∈ (0, 1) be given. We first choose p ∈ N such that |t -t p |	ε. Since the
	sequence (P M m ,δt m t p


Robert Eymard
email: eymard@univ-mlv.fr

Christiane Cocozza-Thivent

APPROXIMATION OF THE MARGINAL DISTRIBUTIONS OF A SEMI-MARKOV PROCESS USING A FINITE VOLUME SCHEME

Keywords: Mathematics Subject Classification. 60K15, 60K20, 65C20, 65M60 Renewal equation, semi-Markov process, convergence of a finite volume scheme

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, •) (resp. b(0, •)). 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 ≥ 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 ūt (1, •) dx (resp. ūt (0, •) dx). The functions ūt (i, •) ∈ L 1 (R + ), for (i, t) ∈ E × R + , are then (non strictly) positive, and must verify, for t ∈ R + , that i∈E R+ ūt (i, x) dx = 1 (thanks to the law of 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) → ūt (i, x) is the weak solution of the following equations

∂ ūt ∂t (i, x) + ∂ ūt ∂x (i, x) = -b(i, x) u t (i, x), for (i, t, x) ∈ E × R + × R + , (1)
with an initial condition ū0 (i, x) = ūini (i, x), for (i, x)

∈ E × R + , (2)
and a boundary condition

ūt (i, 0) = R+ b(1 -i, x)ū t (1 -i, x) dx, for (i, t) ∈ E × R + . (3)
Equations ( 1)-(3) constitute a system of linear hyperbolic equations on R + × 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 [START_REF] Cinlar | Introduction to Stochastic processes[END_REF] or [START_REF] Cocozza-Thivent | Processus stochastiques et fiabilité des systèmes[END_REF]). 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 [START_REF] Cox | Renewal Theory[END_REF] or [START_REF] Feller | An Introduction to Probability Theory and its Applications[END_REF] 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 [START_REF] Mischler | Stability in a nonlinear population maturation model[END_REF] 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 (η t , X t ), where η t ∈ E denotes the state at time t ∈ 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 ∈ R + is the elapsed time since this component is in the state η t . A given initial probability measure for (η 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 ρ t (i, dx) for (i, t) ∈ E × R + , are measures, solution of a system of linear hyperbolic equations coupled by the boundary conditions, generalizing (1)- [START_REF] Cocozza-Thivent | Processus stochastiques et fiabilité des systèmes[END_REF]. 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 ρ t (i, dx), in order to evaluate the probability A i (t) that a component be in state i ∈ E at time t ∈ R + , given by A i (t) = R+ ρ t (i, dx). Matrix formula can be used for the computation of the Laplace transform of these marginal distributions [START_REF] Cocozza-Thivent | Semi-Markov process for reliability studies[END_REF], but the inversion of the Laplace transform is known to induce some stability problems. Otherwise, the phase method [START_REF] Asmussen | Fitting phase-type pistributions via the EM algorithm[END_REF] 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 [START_REF] Fritz | Notes on the Analytic Description and Numerical Calculation of the Time Dependent Availability[END_REF]) 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.

denotes the state of the component for any time t such that T n ≤ t < T n+1 . We assume that T 0 ≤ 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 ∈ N, for all i 0 , i 1 , • • • , i n-1 , i, j ∈ E, for all real values such that s 0 ≤ 0 < s 1 < • • • < s n and for all bounded measurable function f defined on R + , we have

E 1 {Yn+1=j} f (T n+1 -T n ) / T 0 = s 0 , Y 0 = i 0 , T 1 = s 1 , Y 1 = i 1 , . . . , T n-1 = s n-1 , Y n-1 = i n-1 , T n = s n , Y n = i = E(1 {Yn+1=j} f (T n+1 -T n ) / Y n = i) = R+ f (t) q(i, j, t) dt.
We consider the semi-Markov process (η t ) t∈R+ , defined by

η t = Y n for all n ∈ N and for all t ∈ R + such that T n ≤ t < T n+1 . ( 4
)
Intuitively speaking, such a process (η t ) t∈R+ forgets its past at times T n and satisfies

P(η Tn+1 = j, T n+1 -T n ∈ [t, t + dt] / η Tn = i) = P(η T1 = j, T 1 -T 0 ∈ [t, t + dt] / η 0 = i) = q(i, j, t) dt.
We then define the transition rate a by

a(i, j, t) =        0 i f P(T 1 -T 0 > t / η 0 = i) = +∞ t k∈E q(i, k, u) du = 0, q(i, j, t) P(T 1 -T 0 > t / η 0 = i)
otherwise.

(

) 5
The function a intuitively satisfies the following formula

P(η T1 = j, T 1 -T 0 ∈ [t, t + dt] / T 1 -T 0 > t, η 0 = i) = a(i, j, t) dt.
Defining the function b by

b(i, x) = j∈E a(i, j, x), ∀(i, x) ∈ E × R + , (6)
it can be shown that

P(T 1 -T 0 > t / η 0 = i) = exp - t 0 b(i, u) du . ( 7
)
The above equation means that b(i, •) is the hazard rate of T 1 -T 0 knowing {η 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, •) and that the failure periods have the same probability density function q(0, 1, •). Let T 0 ≤ 0 < T 1 < . . . < T n < . . . be the successive times at which the component changes of state. Then (T n ) n≥0 is an alternative renewal process and the process η 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) = q(i, 1 -i, t) +∞ t q(i, 1 -i, s) ds , ∀i = 0, 1, ∀t ∈ R + , (8)
or, equivalently,

q(i, 1 -i, t) = b(i, t) exp - t 0 b(i, s) ds , ∀i = 0, 1, ∀t ∈ R + , ( 9
)
where a(1, 0, •) and a(0, 1, •) are respectively the failure rate and the repair rate of the component.

For k ∈ N, let us denote by C k (E × R + ) the space C k (R + ) E . Let X t be the elapsed time in the current state at time t, given by

X t = t -T n , for all t ∈ [T n , T n+1 ), for all n ∈ N. ( 10
)
Then the process (η t , X t ) t∈R+ is a Markov process taking its values in E × R + . Let P t be its semi-group, defined by

(P t ξ)(i, x) = E (ξ(η t , X t ) / (η 0 , X 0 ) = (i, x)) , ∀ξ ∈ C 0 (E × R + ). (11)
We then have the following proposition.

Proposition 1.2. Let (Y n , T n ) n≥0 , with T 0 ≤ 0 < T 1 < . . . < T n . . ., 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 (η t ) t≥0 , a, (X t ) t≥0 and P t be respectively defined by ( 4), ( 5), [START_REF] Cox | Renewal Theory[END_REF] and [START_REF] Eymard | Finite Volume Methods, Handbook of Numerical Analysis[END_REF]. We assume that a ∈ C 0 (E × E × R + ) and that there exists R 0 > 0 such that P(X 0 ≤ R 0 ) = 1. Then the following properties hold.

(

) For all ξ ∈ C 0 (E × R + ), P t ξ ∈ C 0 (E × R + ). ( 1
) Let L : C 1 (E × R + ) → C 0 (E × R + ) be the operator defined, for all ξ ∈ C 1 (E × R + ) by (i, x) → Lξ(i, x) = j∈E a(i, j, x)(ξ(j, 0) -ξ(i, x)) + ∂ξ ∂x (i, x). 2
Then we get

P t ξ = ξ + t 0 P s Lξ ds, ∀ξ ∈ C 1 (E × R + ), which leads, for all ξ ∈ C 1 (E × R + ), to E(ξ(η t , X t )) = E(ξ(η 0 , X 0 )) + t 0 E(Lξ(η s , X s )) ds.
(3) Therefore, the probability distribution ρ t of (η t , X t ) satisfies

ξ dρ t = i∈E R+ ξ(i, x) ρ t (i, dx) = i∈E R+ ξ(i, x) ρ ini (i, dx) + i∈E ξ(i, 0) t 0 j∈E R+ a(j, i, x) ρ s (j, dx) ds - i∈E t 0 R+ b(i, x) ξ(i, x) ρ s (i, dx) ds + i∈E t 0 R+ ∂ξ ∂x (i, x) ρ s (i, dx) ds, ∀ξ ∈ C 1 (E × R + ), ( 12
)
with ρ ini (•, dx) = ρ 0 (•, dx).
The proof of Proposition 1.2 is straightforward, following the ideas of [START_REF] Cocozza-Thivent | A general framework for some asymptotic reliability formulas[END_REF][START_REF] Cocozza-Thivent | On the marginal distributions of Markov processes used in dynamic reliability[END_REF].

For k ∈ N, let us denote, in the same manner as above, by

C k (E × R + × R + ) the space (C k (R + × R + )) E . Let C k ct (E × R + × R + ) be the set of functions ϕ of C k (E × R + × R + ) with a compact support in time (i.e. ∃T > 0 : ∀ t > T, ∀ i ∈ E, ∀x ∈ R + , ϕ(i, x, t) = 0).
The following result holds. Proposition 1.3. Let E be a finite set, let a ∈ C 0 (E × E × R + ) be given and b be defined by [START_REF] Cocozza-Thivent | A general framework for some asymptotic reliability formulas[END_REF], and let ρ ini be a probability measure on E × R + , such that there exists R 0 > 0 with i∈E [0,R0] ρ ini (i, dx) = 1. Let us assume that there exists, for all t ∈ R + , a probability measure ρ t on E × R + such that [START_REF] Feller | An Introduction to Probability Theory and its Applications[END_REF] is satisfied. Then

0 = i∈E R+×R+ ∂ϕ ∂t (i, x, t) + ∂ϕ ∂x (i, x, t) -b(i, x) ϕ(i, x, t) ρ t (i, dx) dt + i∈E R+ ϕ(i, x, 0) ρ ini (i, dx) + i,j∈E R+×R+ a(j, i, x) ϕ(i, 0, t) ρ t (j, dx) dt, ∀ϕ ∈ C 1 ct (E × R + × R + ). ( 13
)
Moreover, for all ξ ∈ C 0 (E × R), the function t → i∈E R+ ξ(i, x) ρ t (i, dx) is continuous, and ρ 0 (•, dx) = ρ ini (•, dx) holds.

Proof. We proceed in the same way as in [START_REF] Cocozza-Thivent | On the marginal distributions of Markov processes used in dynamic reliability[END_REF]. We first remark that the continuity of t → i∈E R+ ξ(i, x) ρ t (i, dx) and ρ 0 (•, dx) = ρ ini (•, dx) are consequence of [START_REF] Feller | An Introduction to Probability Theory and its Applications[END_REF], thanks to regularizations in C

1 (E × R) of ξ ∈ C 0 (E × R).
We then apply [START_REF] Feller | An Introduction to Probability Theory and its Applications[END_REF] to the function ξ : (i, x) → ∂ϕ ∂t (i, x, t), and we integrate the result on t ∈ R + . We then get [START_REF] Fritz | Notes on the Analytic Description and Numerical Calculation of the Time Dependent Availability[END_REF] thanks to Fubini's theorem.

Then the following uniqueness result holds.

Theorem 1.4. Let E be a finite set, let a ∈ C 0 (E × E × R + ) be given and b be defined by [START_REF] Cocozza-Thivent | A general framework for some asymptotic reliability formulas[END_REF], and let ρ ini be a probability measure on

E × R + , such that there exists R 0 > 0 with i∈E [0,R0] ρ ini (i, dx) = 1. Then there exists at most one measure µ on E × R + × R + such that, for all T > 0, there exists R T > 0 with i∈E T 0 [0,RT ] µ(i, dx, dt) = i∈E T 0 R+ µ(i, dx, dt) < +∞ and which satisfies 0 = i∈E R+×R+ ∂ϕ ∂t (i, x, t) + ∂ϕ ∂x (i, x, t) -b(i, x) ϕ(i, x, t) µ(i, dx, dt) + i∈E R+ ϕ(i, x, 0) ρ ini (i, dx) + i,j∈E R+×R+ a(j, i, x) ϕ(i, 0, t) µ(j, dx, dt), ∀ϕ ∈ C 1 ct (E × R + × R + ). ( 14
)
Moreover, if, for all t ∈ R + , there exists a measure ρ t on E × R + such that for all ξ ∈ C 0 (E × R + ), the function t → i∈E R+ ξ(i, x)ρ t (i, dx) is continuous and µ(i, dx, dt) = ρ t (i, dx) dt, then the measure ρ t is unique too and satisfies [START_REF] Fritz | Notes on the Analytic Description and Numerical Calculation of the Time Dependent Availability[END_REF].

Proof. We adapt the methods of [START_REF] Cocozza-Thivent | On the marginal distributions of Markov processes used in dynamic reliability[END_REF] to the particular case considered here. Let us suppose that there exist two Radon measures satisfying the hypotheses of the above theorem. The measure µ, now defined as the difference between these two solutions, thus verifies

0 = i∈E R+×R+   ∂ϕ ∂t (i, x, t) + ∂ϕ ∂x (i, x, t) -b(i, x) ϕ(i, x, t) + j∈E a(i, j, x) ϕ(j, 0, t)   µ(i, dx, dt), ∀ϕ ∈ C 1 ct (E × R + × R + ). ( 15
)
We have to prove that µ is equal to 0, which is equivalent to the relation

i∈E R+×R+ ψ(i, x, t) µ(i, dx, dt) = 0, ∀ψ ∈ C 0 ct (E × R + × R + ). ( 16
) For a given ψ ∈ C 0 ct (E × R + × R + ), let T > 0 be such that ψ(•, •, t) = 0 for all t ≥ T and let R T > 0 be such that i∈E T 0 [0,RT ] µ(i, dx, dt) = i∈E T 0 R µ(i, dx, dt) < ∞. We prolong by continuity the function a on E × E × R, we define b as in (6) by b(i, x) = j∈E a(i, j, x), for all (i, x) ∈ E × R and we prolong ψ by continuity on E × R × R + . Let us define the function A : [0, T ] × C 0 (E × [-T, R T ]) → C 0 (E × [-T, R T ]), (s, ξ) → A(s, ξ) by (i, y) → (A(s, ξ))(i, y) = b(i, y + s) ξ(i, y) -j∈E a(i, j, y + s) ξ(j, -s) + ψ(i, y + s, s).
Since the function A is continuous with respect to s and Lipschitz-continuous with respect to ξ, we deduce from the Cauchy-Lipschitz theorem the existence and the uniqueness of a function φ

∈ C 0 (E × [-T, R T ] × [0, T ]), such that ∂ φ ∂s ∈ C 0 (E × [-T, R T ] × [0, T ]) and φ(i, y, T ) = 0, ∀(i, y) ∈ E × [-T, R T ], ∂ φ ∂s (i, y, s) = b(i, y + s) φ(i, y, s) - j∈E a(i, j, y + s) φ(j, -s, s) + ψ(i, y + s, s), ( 17
)
∀(i, y, s) ∈ E × [-T, R T ] × [0, T ].
Prolonging φ in order to get that φ

∈ C 0 ct (E × R × R + ) and ∂ φ ∂s ∈ C 0 ct (E × R × R + ) (by continuity, we can take φ(•, •, s) = 0 for all s ≥ T ), let us prove that 0 = i∈E R+×[0,T ]   ∂ φ ∂s (i, x -t, t) -b(i, x) φ(i, x -t, t) + j∈E a(i, j, x) φ(j, -t, t)   µ(i, dx, dt), (18)
which, thanks to the definition (17) of φ, gives (16). Let β 0 be a real nonnegative indefinitely differentiable function with a compact support in (-1, 1), such that R β 0 (x) dx = 1. We consider, for all n ∈ N , the function φn

∈ C 1 ct (E × R × R + ) defined by (i, y, s) → φn (i, y, s) = R φ(i, x, s)β n (y -x) dx, where β n (x) = nβ 0 (nx).
We then get that φn and ∂ φn ∂s respectively converge to φ and

∂ φ ∂s in C 0 (E × [-T, R T ] × [0, T ]
), using the uniform continuity in a compact set of any continuous function. We then apply (15) to the element

ϕ n ∈ C 1 ct (E × R + × R + ) such that ϕ n (i, x, t) = φn (i, x -t, t) for all (i, x, t) ∈ E × R + × R + . Since we have ∂ϕn ∂t (i, x, t) = -∂ φn ∂y (i, x -t, t) + ∂ φn ∂s (i, x -t, t) and ∂ϕn ∂x (i, x, t) = ∂ φn ∂y (i, x -t, t), this produces 0 = i∈E R+×R+   ∂ φn ∂s (i, x -t, t) -b(i, x) φn (i, x -t, t) + j∈E a(i, j, x) φn (j, -t, t)   µ(i, dx, dt). ( 19
)
Passing to the limit n → ∞ in (19), we thus get that (18) holds, which concludes the first part of the proof of the theorem. Let now µ be a measure satisfying the hypotheses of the above theorem. The uniqueness, for all t ∈ R + , of a measure ρ t on E × R + such that for all ξ ∈ C 0 (E × R), the function t → i∈E R ξ(i, x)ρ t (i, dx) is continuous and µ(i, dx, dt) = ρ t (i, dx) dt, immediately results from this continuity property since we can write

i∈E R ξ(i, x)ρ t (i, dx) = lim τ →0 i∈E 1 τ R×(t,t+τ ) ξ(i, x)µ(i, dx, ds).
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 ūini and ū such that ρ ini (i, dx) = ūini (i, x)dx and ρ t (i, dx) = ū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

∂ ūt ∂t (i, x) + ∂ ūt ∂x (i, x) = -b(i, x) ūt (i, x), ∀(i, t, x) ∈ E × R + × R + , (20)
with the initial conditions ū0 (i,

x) = ūini (i, x), ∀(i, x) ∈ E × R + , (21)
and the boundary conditions

ūt (i, 0) = j∈E R+ a(j, i, x)ū t (j, x)dx, ∀(i, t) ∈ E × R + . ( 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 ξ ∈ C 1 (E × 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

The numerical scheme

In order to compute an approximation of the probability distribution ρ t of (η t , X t ) for t ∈ R + , we use a finite volume scheme (see for example [START_REF] Eymard | Finite Volume Methods, Handbook of Numerical Analysis[END_REF] for a general introduction on such numerical methods). Let h > 0 be given and let us divide R + into intervals of length h. The principle of the method is to approximate the measure ρ t (i, dx) by the measure ūh t (i, x) dx, where the function ū is constant on each square [mh, (m + 1)h[×[nh, (n + 1)h[ for all integers n and m, i.e.

ūh t (i, x) = u h n (i, m), ∀(x, t) ∈ [mh, (m + 1)h[×[nh, (n + 1)h[, ∀(m, n) ∈ N × 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 u h n+1 (i, m) -u h n (i, m) h + u h n (i, m) -u h n (i, m -1) = -h b(i, mh) u h n+1 (i, m) for all m ≥ 1,
or equivalently

u h n+1 (i, m) = u h n (i, m -1) 1 + h b(i, mh) for all m ≥ 1, (23)
and

u h n+1 (i, 0) = j∈E m≥1 h a(j, i, mh) u h n+1 (j, m), ( 24
)
where we recall that b(i, mh) = j∈E a(i, j, mh).

The initial condition is given by the natural discretization of measure ρ ini :

u h 0 (i, m) = 1 h [mh,(m+1)h[ ρ ini (i, dx) for all m ≥ 0. ( 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 ∈ E, the space step is given by N i h with N i ∈ N : consider the scheme given by the initialization

u h 0 (i, m) = 1 N i h [mNih,(m+1)Nih[ ρ ini (i, dx) for all m ≥ 0, (26)
and, for all n ∈ N and i ∈ E, if there exists k ∈ N such that n + 1 = kN i , then ũh n+1 (i, 0) = 0 and ũh

n+1 (i, m) = u h n (i, m -1) 1 + h b(i, mN i h) for all m ≥ 1, else ũh n+1 (i, m) = u h n (i, m) 1 + h b(i, mN i h)
for all m ≥ 0.

                   (27)
Then we set, for all n ∈ N and i ∈ E,

u h n+1 (i, m) = ũh n+1 (i, m) for all m ≥ 1, u h n+1 (i, 0) = ũh n+1 (i, 0) + 1 N i j∈E m≥0 N j h a(j, i, mN j h) ũh n+1 (j, m). ( 28
)
Then Scheme (26)-(28) satisfies the mass conservation, i.e.

i∈E m∈N

N i hu h n+1 (i, m) = i∈E m∈N N i hu h 0 (i, m) = 1.
Moreover, it keeps respecting the non-diffusive effect resulting from the condition CFL= 1, even for all the states i ∈ 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 [START_REF] Cocozza-Thivent | Méthodologie et algorithmes pour la quantification de petits systèmes redondants[END_REF] to different examples, where keeping the same space step for all states would either produce inaccurate results, either lead to inadmissible computing times.

In [START_REF] Cocozza-Thivent | A numerical scheme to solve integro-differential equations in the dynamic reliability field[END_REF], 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 ∈ Z instead of m ∈ N. The numerical scheme ( 23)-(25) can then be written in the following way

a h (j, i, m) = a(j, i, mh), ∀i, j ∈ E, ∀m ∈ Z, m ≥ 1, a h (j, i, m) = 0, ∀i, j ∈ E, ∀m ∈ Z, m < 1, b h (i, m) = j∈E a h (i, j, m), ∀i ∈ E, ∀m ∈ Z, u h 0 (i, m) = 1 h [mh,(m+1)h[ ρ ini (i, dx), ∀i ∈ E, ∀m ∈ N, u h 0 (i, m) = 0, ∀i ∈ E, ∀m ∈ Z \ N, (1 + h b h (i, m)) u h n+1 (i, m) = u h n (i, m -1) + U h n+1 (i, m), ∀i ∈ E, ∀(m, n) ∈ Z × N, U h n+1 (i, m) = 0, ∀i ∈ E, ∀m ∈ Z \ {0}, ∀n ∈ N, U h n+1 (i, 0) = j∈E m∈Z h a h (j, i, m) u h n+1 (j, m), ∀i ∈ E, ∀n ∈ N, ūh t (i, x) = u h n (i, m), ∀(x, t) ∈ [mh, (m + 1)h) × [nh, (n + 1)h), ∀i ∈ E, ∀(m, n) ∈ Z × N. (29)
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 ∈ C 0 b (E × E × R + ). Let ρ ini be a probability measure on the set E × R + , such that there exists R 0 > 0 with i∈E [0,R0] ρ ini (i, dx) = 1. Let h > 0 be given, and let (u h n (i, m)) (i,m,n)∈E×Z×N and ūh : E × R × R + → R be given by Scheme (29). Then

u h n (i, m) ≥ 0, ∀i ∈ E, ∀n ∈ N, ∀m ∈ Z, ( 30
)
u h n (i, m) = 0, ∀i ∈ E, ∀n ∈ N, ∀m ∈ Z, m < 0, ( 31
)
h i∈E m∈Z u h n (i, m) = 1, ∀n ∈ N, ( 32
)
and i∈E R ūh t (i, x) dx = 1, ∀t ∈ R + . ( 33
)
Moreover, denoting, for any positive real x, by [x] the greatest integer lower or equal to x, then

u h n (i, m) = 0, ∀i ∈ E, ∀n ∈ N, ∀m ∈ N, m [R 0 /h] + n. ( 34
)
Proof. Properties (30)-(32) hold for n = 0 since ρ ini is a probability measure with support on E × [0, R 0 ]. Assume that these properties hold for a given n ∈ N. We get from Scheme (29) that Properties (30) and (31) are verified for m + 1. Moreover

h i∈E m∈Z,m≥1 (1 + h b h (i, m))u h n+1 (i, m) = h i∈E m∈Z,m≥0 u h n (i, m) = 1.
It delivers, since we can write

u h n+1 (i, 0) = U h n+1 (i, 0) = j∈E m∈Z h a h (j, i, m) u h n+1 (j, m), ∀i ∈ E, that h i∈E u h n+1 (i, 0) = h i∈E m∈Z,m≥1 h b h (i, m)u h n+1 (i, m),
which gives

h i∈E m∈Z u h n+1 (i, m) = h i∈E m∈Z,m≥1 (1 + h b h (i, m))u h n+1 (i, m) = 1.
This concludes the proof of (32), which immediately produces (33). The proof of (34) then comes from a proof by induction.

Convergence of the numerical scheme

Let us begin with proving a property of continuity with respect to the time of the measures ūh t dx.

Lemma 2.3. Let E be a finite set, and let a ∈ C 0 (E × E × R + ). Let ρ ini be a probability measure on the set E × R + , such that there exists R 0 > 0 with i∈E [0,R0] ρ ini (i, dx) = 1. Let h ∈ (0, 1) be given, and let (u h n (i, m)) (i,m,n)∈E×Z×N and ūh : E × R × R + → R be given by Scheme (29). Let ψ ∈ C 1 (E × R). Then, for all T ∈ R + , there exists a real C 1 (ψ, T ), which only depends on ψ, T , a and ρ ini , such that

i∈E R ψ(i, x) ūh t (i, x) dx - i∈E R ψ(i, x) ūh s (i, x) dx ≤ C 1 (ψ, T ) (|t -s| + h), ∀s, t ∈ [0, T ]. ( 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 ≤ s ≤ t ≤ T . Let us define

∆ h (s, t) = i∈E R ψ(i, x) ūh t (i, x) dx - i∈E R ψ(i, x) ūh s (i, x) dx.
We then have

∆ h (s, t) = i∈E m∈Z h ψ(i, m) u h [t/h] (i, m) - i∈E m∈Z h ψ(i, m) u h [s/h] (i, m),
where we set, for all m ∈ Z,

ψ(i, m) = 1 h (m+1)h mh ψ(i, x) dx. We get ∆ h (s, t) = [t/h]-1 n=[s/h] i∈E m∈Z h ψ(i, m) [u h n+1 (i, m) -u h n (i, m)].
From Scheme (29), we obtain

i∈E m∈Z h ψ(i, m) u h n+1 (i, m) -u h n (i, m) = i∈E m∈Z h ψ(i, m) -u h n (i, m) + u h n (i, m -1) -h b h (i, m) u h n+1 (i, m) + U h n+1 (i, m) = -h i∈E m∈Z ψ(i, m + 1) -ψ(i, m) u h n (i, m) -h 2 i∈E m∈Z ψ(i, m) b h (i, m) u h n+1 (i, m) + h 2 i∈E ψ(i, 0) j∈E m∈Z a h (j, i, m) u h n+1 (j, m).
Hence, we deduce

∆ h (s, t) = -h i∈E m∈Z ψ(i, m + 1) -ψ(i, m) [t/h]-1 n=[s/h] u h n (i, m) -h 2 i∈E m∈Z ψ(i, m) b h (i, m) [t/h]-1 n=[s/h] u h n+1 (i, m) + h 2 i∈E ψ(i, 0) j∈E m∈Z a h (j, i, m) [t/h]-1 n=[s/h] u h n+1 (j, m). Since ψ(i, m + 1)-ψ(i, m) = 1 h (m+1)h mh h 0 ∂ψ ∂x (i, x + y) dy dx and since [t/h]-1 n=[s/h] u h n (i, m) = 0 implies 0 ≤ mh ≤ (m + 2)h ≤ R T , we have |ψ(i, m + 1) -ψ(i, m)| [t/h]-1 n=[s/h] u h n (i, m) ≤ max i∈E sup x∈[0,TT ] ∂ψ ∂x (i, x) [t/h]-1 n=[s/h] u h n (i, m).
Using similar arguments for the other terms, we deduce that

|∆ h (s, t)| ≤ h 2 max i∈E sup x∈[0,RT ] ∂ψ ∂x (i, x) [t/h]-1 n=[s/h] i∈E m∈Z u h n (i, m) + h 2 max i∈E sup x∈[0,RT ] |ψ(i, x)| sup x∈[0,RT ] |b(i, x)| [t/h] n=[s/h]+1 i∈E m∈Z u h n (i, m) + h 2 max i∈E sup x∈[0,1] |ψ(i, x)| max i∈E sup x∈[0,RT ] |b(i, x)| [t/h] n=[s/h]+1 i∈E m∈Z u h n (i, m).
Using (32), we obtain

|∆ h (s, t)| ≤ C 1 (ψ, T ) h ([t/h] -[s/h]) ≤ C 1 (ψ, T ) (t -s + h),
where C 1 (ψ, T ) is a real which only depends on ψ, T , a and ρ ini .

We can now prove the following theorem.

Theorem 2.4. Let E be a finite set, and let a ∈ C 0 (E × E × R + ). Let ρ ini be a probability measure on the set E × R + , such that there exists R 0 > 0 with i∈E [0,R0] ρ ini (i, dx) = 1. Let (h p ) p∈N be a sequence of strictly positive real numbers, which converges to 0 as p → ∞. Then there exists a subsequence of (h p ) p∈N , again denoted (h p ) p∈N , and, for all t ≥ 0, there exists a probability measure ν t on E × R + with a finite support, such that the sequence of measures ((ū hp t (i, x) dx, i ∈ E)) p∈N converges to ν t as p → ∞, where ūhp : E × R × R + → R is given by Scheme (29) for h = h p . The measure convergence means that for all ξ ∈ C 0 (E × R + ), we have

lim p→+∞ i∈E E×R+ ξ(i, x)ū hp t (i, x) dx = E×R+ ξ dν t .
Moreover, for all ξ ∈ C 0 (E × R + ), the function t → E×R+ ξ dν t is continuous.

Proof. Without lack of generality, we assume that h p ∈ (0, 1) for all p ∈ N. Let (t m ) m∈N be a sequence of real numbers, which is dense into R + . For any m ∈ N and p ∈ N, we deduce from ( 32) and (34) that the measure ūhp tm dx is a probability measure on E × 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∈N (depending on m) and a probability measure ν tm such that the sequence of measures ((ū hp tm (i, x) dx, i ∈ E)) p∈N converges to ν tm as p → ∞ in the following sense: for all bounded functions ψ ∈ C 0 (E × R + ):

lim p→+∞ i∈E E×R+ ψ(i, x)ū hp tm (i, x) dx = E×R+ ψ dν tm .
Using the diagonal method, we then can build a subsequence, again denoted (h p ) p∈N such that for all m ∈ N, the sequence (ū hp tm (i, x) dx) p∈N converges to the measure ν tm in the sense given above. Let t ∈ R + and let ξ ∈ C 0 (E × R + ) be given. Let us prove that ( i∈E R ξ(i, x)ū hp t (i, x) dx) p∈N converges as p → ∞. To that purpose, we prove that it is a Cauchy sequence. We set T = t + 1 and R T = R 0 + T + 1. Let ε > 0 be given. We approximate function ξ by a bounded function

ψ ∈ C 1 (E × R + ) such that sup x∈[0,RT ] |ξ(x) -ψ(x)| ≤ ε. Let m ∈ N be such that |t m -t| ≤ min(ε/C 1 (ψ, T ), 1
). We then obtain for all p ∈ N, using Lemma 2.3,

i∈E R ψ(i, x) (ū hp t (i, x) - ūhp tm (i, x)) dx ≤ C 1 (ψ, T ) (|t -t m | + h p ).
We now choose n ∈ N such that, for all p, q ≥ n, i∈E R ψ(i, x) (ū hp tm (i, x)ūhq tm (i, x)) dx ≤ ε (this is possible since the sequence ( i∈E R ψ(i, x)ū hp tm (i, x) dx) p∈N converges) and h p and h q less or equal to ε/C 1 (ψ, T ). We then get, for p, q ≥ n,

i∈E R ξ(i, x) (ū hp t (i, x) - ūhq t (i, x)) dx ≤ i∈E R ψ(i, x) (ū hp t (i, x) - ūhq t (i, x)) dx + 2ε ≤ i∈E R ψ(i, x) (ū hp t (i, x) - ūhp tm (i, x)) dx + i∈E R ψ(i, x) (ū hp tm (i, x) - ūhq tm (i, x)) dx + i∈E R ψ(i, x) (ū hq t (i, x) - ūhq tm (i, x)) dx + 2ε ≤ 7ε.
This proves that ( i∈E R ξ(i, x)ū hp t (i, x) dx) p∈N is a Cauchy sequence. Let us define

S t (ξ) = lim p→∞ i∈E R ξ(i, x)ū hp t (i, x) dx.
The operator S t is positive linear form on C 0 (E × R + ) and S t (1) = 1. Consequently it is a probability denoted by ν t , i.e. S t (ξ) = ξ dν t = i∈E R+ ξ(i, x) ν t (i, dx).

Let us prove the continuity of t → ξ dν t , for any ξ ∈ C 0 (E × R × R + ). We again set T = t + 1 and 35), and we let p → ∞. We thus get, for s

R T = R 0 +T +1. Let ε > 0. We again approximate ξ by ψ ∈ C 1 (E ×R + ) such that sup x∈[0,RT ] |ξ(x)-ψ(x)| ≤ ε. We let h = h p in (
∈ [0, T ], i∈E R ψ(i, x) ρ t (i, dx) - i∈E R ψ(i, x) ρ s (i, dx) ≤ C 1 (ψ, T ) |t -s|. We then have, for all s ∈ [0, T ] such that |t -s| ≤ ε/C 1 (ψ, T ), i∈E R ξ(i, x) ν t (i, dx) - i∈E R ξ(i, x) ν s (i, dx) ≤ i∈E R ψ(i, x) ν t (i, dx) - i∈E R ψ(i, x) ν s (i, dx) + 2ε ≤ 3ε,
which completes this proof of continuity.

We now get the following corollary.

Corollary 2.5. Let be a finite set, and let a ∈ C 0 (E × E × R + ). Let ρ ini be a probability measure on the set E × R + , such that there exists R 0 > 0 with i∈E [0,R0] ρ ini (i, dx) = 1. Let (h p ) p∈N be a sequence of strictly positive real numbers, which converges to 0 as p → ∞. Let us again denote by (h p ) p∈N the subsequence of (h p ) p∈N and, for all t ≥ 0, by ν t the probability measure with a finite support given by Theorem 2.4. Then the sequence of measures (ū hp t (i, x) dx dt) p∈N , where ūhp : E × R × R + → R is given by Scheme (29) for h = h p , converges to the measure µ defined by µ(i, dx, dt) = ν t (i, dx) dt as p → ∞ in the following sense: for all

ξ ∈ C0 ct (E × R × R + ) we have lim p→∞ i∈E R+ R+ ξ(i, x, t) ūhp t (i, x) dx dt = i∈E R+ R+ ξ(i, x, t) ν t (i, dx) dt = E×R+×R+ ξ dµ. Moreover, let ξ ∈ C0 ct (E × R × R + ) and A be such that ξ(i, x, t) = 0, ∀ t > A, ∀ i ∈ E, ∀ x ∈ R, and let R A = R 0 + A + 1. We then have lim p→+∞ h 2 p i∈E m∈Z [A/hp]-1 n=0 ξ(i, mh, nh) u h n (i, m) = E×R+×R+ ξ dµ. ( 36
) Proof. Let ξ ∈ C 0 (E × R × R + ). We get that lim p→∞ i∈E R+ R+ ξ(i, x, t) ūhp t (i, x) dx dt = i∈E R+ R+ ξ(i, x, t)ρ t (i, dx) dt,
thanks to the Lebesgue's dominated convergence theorem.

Let us now prove (36). We can write

i∈E R+×R+ ξ(i, x, t) ūh t (i, x) dx dt = i∈E m≥0 n≥0 [mh,(m+1)h[×[nh,(n+1)h[ ξ(i, x, t) u h n (i, m) dx dt = i∈E m≥0 [A/h]-1 n=0 u h n (i, m) [mh,(m+1)h[×[nh,(n+1)h[ ξ(i, x, t) dx dt + i∈E m≥0 u h n (i, m) [mh,(m+1)h[×[h[A/h],A] ξ(i, x, t) dx dt.
On the one hand:

i∈E m≥0 u h n (i, m) [mh,(m+1)h[×[h[A/h],A] ξ(i, x, t) dx dt ≤ i∈E m≥0 u h n (i, m) ξ L ∞ (E×R+×R+) h (A -h[A/h]) = ξ L ∞ (E×R+×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

δ ξ (h) = h 2 i∈E m∈Z [A/h]-1 n=0 ξ(i, mh, nh) u h n (i, m) - i∈E m≥0 [A/h]-1 n=0 u h n (i, m) [mh,(m+1)h[×[nh,(n+1)h[ ξ(i, x, t) dx dt ≤ i∈E m≥0 [A/h]-1 n=0 u h n (i, m) [mh,(m+1)h[×[nh,(n+1)h[ |ξ(i, mh, nh) -ξ(i, x, t)| dx dt.
Using (34) we obtain

δ ξ h) ≤ i∈E [(R0+A)/h]-1 m=0 [A/h]-1 n=0 u h n (i, m) [mh,(m+1)h[×[nh,(n+1)h[ |ξ(i, mh, nh) -ξ(i, x, t)| dx dt. Let us define A,ξ (h) = sup{|ξ(i, x 1 , t 1 ) -ξ(i, x 2 , t 2 )| : i ∈ E, 0 ≤ x 1 , x 2 ≤ R 0 + A, 0 ≤ t 1 , t 2 ≤ T, |x 1 -x 2 | ≤ h, |t 1 -t 2 | ≤ h}.
Then we have

0 ≤ δ ξ (h) ≤ i∈E [(R0+A)/h]-1 m=0 [A/h]-1 n=0 A,ξ (h) h2 u h n (i, m) ≤ A,ξ (h) [A/h] ≤ A A,ξ (h)
and the right hand side of the above inequalities tends to 0 when h tends to O, which ends the proof.

Let us now prove that the measure µ(i, dx, dt) = ν t (i, dx) dt satisfies [START_REF] Mischler | Stability in a nonlinear population maturation model[END_REF].

Proposition 2.6. Let E be a finite set, and let a ∈ C 0 (E × E × R + ). Let ρ ini be a probability measure on the set E × R + , such that there exists R 0 > 0 with i∈E [0,R0] ρ ini (i, dx) = 1. Let (h p ) p∈N be a sequence of strictly positive real numbers, which converges to 0 as p → ∞. Let us again denote by (h p ) p∈N the subsequence of (h p ) p∈N and, for all t ≥ 0, by ν t the probability measure with a finite support given by Theorem 2.4. Then the measure µ(i, dx, dt) = ρ t (i, dx) dt is such that for all T > 0, there exists

R T > 0 with i∈E T 0 [0,RT ] µ(i, dx, dt) = i∈E T 0
R+ µ(i, dx, dt) < ∞ and that ( 14) is satisfied.

Proof. Let us assume the assumptions of the above proposition, and let

ϕ ∈ C 1 ct (E × R × R + )
. Without lack of generality, we can assume that h p ∈ (0, 1). We consider a time T > 1 such that

ϕ(i, x, t) = 0 ∀ (i, x, t) ∈ E × R × [T -1, +∞), (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 ∈ N is given, h p by h. We also set

µ h (i, dx, dt) = ūh t (i, x) dx dt, ∀(i, x, t) ∈ E × R × R + . ( 38
)
Let us notice that (33) gives

E×R ×[0, T ] dµ = E×R+ ×[0, T ] dµ = T. ( 39
)
Let ∈ (0, 1) and h ∈ (0, 1 -). Let us first prove that

E×R×[0,T ] dµ h = E×R+×[0,T ] dµ h = E×[0,RT -]×[0,T ] dµ h , ( 40
) E×R+×[0,T ] dµ = E×[0,RT ]×[0,T ]
dµ.

We have

0 ≤ E×[RT -,+∞[×[0,T ] dµ h ≤ m≥ R T - h [T /h] n=0 h2 u h n (i, m).
Then Formula (34) gives

u h n (i, m) = 0 for 0 ≤ n ≤ [T /h], and m > [R 0 /h] + [T /h]. It is easy to check that [(R T -)/h] < [R 0 /h] + [T /h]. Therefore we get E×[RT -,+∞[×[0,T ] dµ h = 0, which proves the result for measure µ h . Now, let > 0 and let ϕ ∈ C 0 (E × R + × R + ) a function which satisfies: 0 ≤ ϕ ≤ 1, ϕ = 1 on E × [R T , +∞[×[0, T ] and ϕ = 0 on E × [R T -, +∞[×[0, T + ]. We obtain 0 ≤ E×[RT ,+∞[×[O,T ] dµ ≤ ϕ dµ = lim h→0 ϕ dµ h . For h < 1 -, we have 0 ≤ ϕ dµ h ≤ E×[RT -,+∞[×[T,T + ] dµ h ≤ T + T i∈E R+ ūh t (i, x) dx dt ≤ .
Let us now prove that measure µ satisfies [START_REF] Mischler | Stability in a nonlinear population maturation model[END_REF]. We get, from Scheme (29),

0 = u h n+1 (i, m) -u h n (i, m) + u h n (i, m) -u h n (i, m -1) + h b h (i, m) u h n+1 (i, m) -U h n+1 (i, m).
Multiplying the above equation by -h ϕ(i, mh, nh) and summing on i ∈ E, m ∈ Z and 0 ≤ n ≤ [T /h] -1, we obtain

T h 1 + T h 2 + T h 3 + T h 4 = 0,
with

T h 1 = - i∈E m∈Z [T /h]-1 n=0 (u h n+1 (i, m) -u h n (i, m)) h ϕ(i, mh, nh), T h 2 = - i∈E m∈Z [T /h]-1 n=0 (u h n (i, m) -u h n (i, m -1)) h ϕ(i, mh, nh), T h 3 = - i∈E m∈Z [T /h]-1 n=0 h 2 b h (i, m) u h n+1 (i, m) ϕ(i, mh, nh), T h 4 = i∈E [T /h]-1 n=0 h U h n+1 (i, 0) ϕ(i, 0, nh).
Step 1. Study of T h 1

We have

T h 1 = -h i∈E m∈Z [T /h] n=1 u h n (i, m)ϕ(i, mh, (n -1)h) + h i∈E m∈Z [T /h]-1 n=0 u h n (i, m)ϕ(i, mh, nh) = h i∈E m∈Z [T /h]-1 n=1 u h n (i, m) (ϕ(i, mh, nh) -ϕ(i, mh, (n -1)h)) + h i∈E m∈Z u h 0 (i, m) ϕ(i, mh, 0), since, if h is small enough, i.e. if p is large enough, we have ϕ(i, mh, ([T /h] -1)h) = 0.
Let us define T p 1 and T p 1 by

T p 1 = h i∈E m∈Z [T /h]-1 n=1 u h n (i, m) (ϕ(i, mh, nh) -ϕ(i, mh, (n -1)h)), T p 1 = h i∈E m∈Z u h 0 (i, m) ϕ(i, mh, 0).
We also define ψ h by

ψ h (i, m, n) = ϕ(i, mh, nh) -ϕ(i, mh, (n -1)h) h if n = 0, ψ h (i, m, 0) = ∂ϕ ∂t (i, mh, 0). and ψh by ψh (i, x, t) = ψ h (i, m, n) if (x, t) ∈ [mh, (m + 1)h[×[nh, (n + 1)h[, ∀m ∈ Z, 0 ≤ n ≤ [T /h] -1.
For n ≥ 1, we notice that there exists

t h n ∈ [(n -1)h, nh[ such that ψ h (i, m, n) = ∂ϕ ∂t (i, mh, t h n ).
Since the function ∂ϕ ∂t is continuous, it is uniformly continuous on [0, R T ] × [0, T ], and therefore

lim p→∞ ∂ϕ ∂t -ψhp L ∞ (E×[0,RT ]×[0,T ]) = 0. (41)
We have

E×R×R+ ψh dµ h = i∈E m∈Z [T /h]-1 n=0 h 2 ψ h (i, m, n) u h n (i, m) + E×R×[h[T /h],T ] ψh dµ h = T p 1 + h 2 i∈E m∈Z u h 0 (i, m) ψ h (i, m, 0) + E×R×[h[T /h],T ] ψh dµ h . ( 42
)
On the one hand, using (37), ( 40) and (39), we get On the other hand, using (32) and (39) we obtain

E×R×R+ ψh dµ h - E×R×R+ ∂ϕ ∂t dµ ≤ E×R×[0,T ) ψh - ∂ϕ ∂t dµ h + E×R×R+ ∂ϕ ∂t dµ h - E×R×R+ ∂ϕ ∂t dµ ≤ T ∂ϕ ∂t -ψh L ∞ (E×[0,RT ]×[0,T ]) + E×R×R+ ∂ϕ ∂t dµ h - E×R×R+ ∂ϕ ∂t dµ . Since ∂ϕ ∂t ∈ C0(E × R × R + ), Corollary 2.
h 2 i∈E m∈Z u h 0 (i, m) ψ h (i, m, 0) + E×R×[h[T /h],T ] ψh dµ h ≤ h ψh L ∞ (E×[0,RT ]×[0,T ]) h i∈E m∈Z u h 0 (i, m) + h ψh L ∞ (E×[0,RT ]×[0,T ]) = 2 h ψh L ∞ (E×[0,RT ]×[0,T ]) .
But (41) implies that ψh

(i, m) = R f (x) ρ ini (i, dx).
Indeed, we can write

R f (x) ρ ini (i, dx) -h m∈Z f (mh) u h 0 (i, m) = R f (x) ρ ini (i, dx) - m∈Z [mh,(m+1)h[ f (mh) ρ ini (i, dx) ≤ m∈Z [mh,(m+1)h[ |f (x) -f (mh)| ρ ini (i, dx) ≤ sup x 1 , x 2 ∈ [0, R 0 ] |x 1 -x 2 | ≤ h |f (x 1 ) -f (x 2 )| m∈Z [mh,(m+1)h[ ρ ini (i, dx) = su p x 1 , x 2 ∈ [0, R 0 ] |x 1 -x 2 | ≤ h |f (x 1 ) -f (x 2 )|,
and the uniform continuity of f on [0, R 0 ] implies the convergence to 0 when p tends to ∞. Therefore we get

lim p→∞ T p 1 = i∈E R ϕ(i, x, 0) ρ ini (i, dx).
Gathering the previous results, we have proven that

lim p→∞ T h 1 = E×R+×R+ ∂ϕ ∂t dµ + i∈E R+ ϕ(i, x, 0) ρ ini (i, dx). ( 46
)
Step 2. Study of T h 2

We have

T h 2 = -h i∈E m∈Z [T /h]-1 n=0 u h n (i, m)ϕ(i, mh, nh) + h i∈E m∈Z [T /h]-1 n=0 u h n (i, m)ϕ(i, (m + 1)h, nh) = h i∈E m∈Z [T /h]-1 n=0 u h n (i, m) (ϕ(i, (m + 1)h, nh) -ϕ(i, mh, nh)).
Using the same methods as in step 1, we show that

lim p→∞ T h 2 = E×R+×[0,T ) ∂ϕ ∂x dµ.
Step 3. Study of T h 3

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

T h 3 = -h 2 i∈E m∈Z [T /h]-1 n=1 b h (i, m) u h n (i, m) ϕ(i, mh, (n -1)h).
We notice that lim h→0

T h 3 = lim h→0 T h 3 with T h 3 = -h 2 i∈E m∈Z [T /h]-1 n=1 b h (i, m) u h n (i, m) ϕ(i, mh, nh),
since, using (32) and (34) we obtain

h 2 i∈E m∈Z [T /h]-1 n=1 b h (i, m) u h n (i, m) | ϕ(i, mh, nh) -ϕ(i, mh, (n -1)h) | = h 2 i∈E m∈Z, mh≤R T [T /h]-1 n=1 b h (i, m) u h n (i, m) | ϕ(i, mh, nh) -ϕ(i, mh, (n -1)h) | ≤ sup x∈[0,RT ] sup |t1-t2|≤h |ϕ(i, x, t 1 ) -ϕ(i, x, t 2 | b L ∞ (E×[0,RT ]) h 2 i∈E [T /h]-1 n=0 m∈Z u h n (i, m) ≤ sup x∈[0,RT ] sup |t1-t2|≤h |ϕ(i, x, t 1 ) -ϕ(i, x, t 2 )| b L ∞ (E×[0,RT ]) T,
and the uniform continuity of ϕ on [0, R T ] × [0, T ] gives the result.

In fact, the following relation holds

T h 3 = -h 2 i∈E m≥0 [T /h]-1 n=0 b h (i, m) u h n (i, m) ϕ(i, mh, nh) + h 2 i∈E m∈Z b h (i, m) u h 0 (i, m) ϕ(i, mh, 0).
On the one hand, we have, using (36)

lim p→∞ h 2 p i∈E m≥0 [T /hp]-1 n=0 b hp (i, m) ϕ(i, mh p , nh p ) u hp n (i, m) = E×R+×R+ b ϕ dµ.
On the other hand, thanks to (32) and (34), we have We have

|h 2 i∈E m∈Z b h (i, m) u h 0 (i, m) ϕ(i, mh, 0)| ≤ b L ∞ ([0,RT ]) ϕ L ∞ (E×[0,RT ]×[0,T ]) h 2 i∈E m∈Z u h 0 (i, m) = b L ∞ ([0,RT ]) ϕ L ∞ (E×[0,RT ]×[0,T ]) h,
T h 4 = i∈E [T /h]-1 n=0 h j∈E m∈Z h a h (j, i, m) u h n+1 (j, m) ϕ(i, 0, nh).
and therefore as previously

lim h→0 T h 4 = i∈E j∈E R×[0,T )
a(i, j, x) ϕ(j, 0, t) µ(i, dx, dt)

Step 5. Synthesis Since

0 = T h 1 + T h 2 + T h 3 + T h 4 , we have 0 = lim h→0 (T h 1 + T h 2 + T h 3 + T h 4 ) = i∈E R+×R+ ∂ϕ ∂t (i, x, t) µ(i, dx, dt) + i∈E R+ ϕ(i, x, 0) ρ ini (i, dx) + i∈E R+×R+ ∂ϕ ∂x (i, x, t) µ(i, dx, dt) - i∈E R+×R+ b(i, x) ϕ(i, x, t) µ(i, dx, dt) + i∈E j∈E R+×R+ a(i, j, x) ϕ(j, 0, t) µ(i, dx, dt),
and this is the requested result.

Theorem 2.7. Let E be a finite set, and let a ∈ C 0 (E × E × R + ). Let ρ ini be a probability measure on the set E × R + , such that there exists R 0 > 0 with i∈E [0,R0] ρ ini (i, dx) = 1. Let h > 0 be a strictly positive real number. Then, for all t ∈ R + , the probability measure ūh t (i, x) dx converges to ρ t (i, dx), the unique solution of [START_REF] Feller | An Introduction to Probability Theory and its Applications[END_REF] as h tends to 0.

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

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 [START_REF] Cocozza-Thivent | A numerical scheme to solve integro-differential equations in the dynamic reliability field[END_REF] and [START_REF] Cocozza-Thivent | Méthodologie et algorithmes pour la quantification de petits systèmes redondants[END_REF]. The mean value of the duration in the state i, for i = 0, 1, is then given by m i = R+ t q(i, 1i, t) dt, the standard deviation is given by σ 2 i = R+ t 2 q(i, 1i, t) dtm 2 i and the coefficient of variation C v i is given by C v i = σ i /m i (see for example [START_REF] Cocozza-Thivent | Processus stochastiques et fiabilité des systèmes[END_REF]). We thus approximate the marginal distributions ρ 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) = P(η t = 1) = R+ ρ t (1, dx) and the unavailability is defined by A 0 (t) = R+ ρ 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→∞ A 0 (t) = m 0 m 1 + m 0 • (47)

We examine three numerical examples of the calculation of the unavailability A 0 (t) = P(η 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, 1i, •) kept for both states of the three examples. We also give the mean values m i and the coefficients of variation C v i .

q(1, 0, •) The Weibull distribution, with parameters α and β, is given by q(t) = b(t) exp - the log-normal distribution with parameters α and β, is given by

m 1 C v 1 = σ 1 /m 1 q(0, 1, •) m 0 C v 0 = σ 0 /m 0 Ex
q(t) = 1 √ 2π β t exp -(log t -α) 2 2β 2
, b(t) = q(t)

+∞ t q(s) ds , m = e α+β 2 /2 , σ 2 = e 2α e 2β 2e β 2 , and the gamma distribution, with parameters α and β, is given by q(t) = 1 Γ(α) β α t α-1 e -t/β , b(t) = q(t) +∞ t q(s) ds , m = αβ, σ 2 = αβ 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 ū80/3 (i, x) for i = 0 and i = 1, as a function of x, where ū80/3 (i, x) dx is an approximation of ρ 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 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 [START_REF] Cocozza-Thivent | Marginal distributions of a semi-Markov process and their computations[END_REF][START_REF] Cocozza-Thivent | Méthodologie et algorithmes pour la quantification de petits systèmes redondants[END_REF].

  5 and (41) imply that lim p→∞ E×R×[0,T ) ψhp dµ hp = E×R×[0,T )

Step 4 .

 4 therefore this term tends to 0 when p tends to ∞.We have therefore shown that lim p→∞ Study of T h 4

t 0 b 1 β , σ 2 = α 2 Figure 1 .

 01221 Figure 1. Example 1. Left: numerical value of the unavailability A 0 (t) as a function of the time, right: approximation of ū80/3 (i, x).

Figure 2 .

 2 Figure 2. Numerical values of the unavailability A 0 (t) as a function of the time, Example 2.

Figure 3 .

 3 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.)