A finite-volume scheme for dynamic reliability models

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


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 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 Davis (1984Davis ( , 1993), the family of marginal distributions (ρ(t)(i, dx)) t∈R + is solution of the Chapman-Kolmogorov equations 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)v(i, x) ( j, i, y)u( j, y, t) M( j, i, y, x)dy − u(i, x, t) 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 Cocozza-Thivent & Eymard, 2004). Hence, the functions u and u 0 satisfy in a weak sense the following system of linear scalar hyperbolic equations, only coupled by the boundary condition on {0} × R + : 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 Sinko & Streifer (1967) and Bell & Anderson (1967). In this model, u(1, x, t) represents the density of population of age x at the time t, v(1, x) is the specific growth rate, b(1, x) is the specific mortality at the age x, a(1, 1, x) is the specific fertility at the age x and u 0 (1, x) is the initial density of the population. This model allows for the computation of all the informations concerning the time evolution of the population and has been the object of numerous studies and adaptations to more complex situations (see Mischler et al., 2002, for non-linear improvements, and references therein, and see Feller, 1966, for the probabilistic notions). For such models, whether dynamic reliability or population dynamics model, a challenging problem is to calculate the marginal distribution ρ(t)(·, dx) of the process (I t , X t ). Actually, as can be seen in the literature (see, e.g. Labeau, 1996), such a distribution is analytically calculable in only very simple cases and only numerical approximations are usually possible. One of the most up-to-date methods for such numerical computations appeals to Monte Carlo simulation (see, e.g. Labeau, 1996, and other papers by P. E. Labeau).
We construct here an approximation of ρ(t)(i, dx) using a finite-volume method and we prove that such an approximation weakly converges towards the unique solution of (1.2). Such a method has already been studied and implemented in some restrictive cases of dynamic reliability (see Cocozza-Thivent & Eymard and Cocozza-Thivent & Eymard (2004) for the 'semi-Markov' case). Note that similar numerical approximations for (1.4) have been studied in the 1D case (see Ackleh et al. (2002), Ackleh & Ito (1997) and Abia et al. (2004), and references therein for the population dynamics problem). The study of the scheme on unstructured meshes in 2D (or more) is much more difficult, and the proofs of convergence for the approximation of a linear or non-linear scalar hyperbolic problem, using a finite-volume method, are more recent (they rely on the notion of weak bounded variation inequality, see, e.g. Eymard et al., 2000). Nevertheless, it appears that the suitable functional framework for the study of (1.2) is the set of continuous functions of the time variable valued in P(E × 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: 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: 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: which leads to We then introduce the notion of admissible mesh M of R d , i.e. a family of measurable subsets of R d such that 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 We define the following coefficients for the convection part of the scheme, using the notations given at the beginning of this section: (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 (2.8) we consider the weak solution u of the partial differential equation This function u is then such that u( , 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, (2.9) We now take into account the right-hand side of the partial differential equations (1.3), by setting 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 (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: . 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.

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

LEMMA 3.1 Under Hypothesis (H-3), we have
Proof. Let t 0 and x ∈ R d be given. We can write LEMMA 3.2 (MASS 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 (i,K ) n ) n∈N,i∈E,K ∈M which is the solution of (2.4)-(2.10). Moreover, this family is such that Proof. The existence and uniqueness of the solution ( p 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

This leads to
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., Billingsley (1968)), which allows us to introduce the function f i, j : R + × R + → R + , given by We get the following result.
Then the following properties hold: (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(V 1 δt) − 1 V 1 δt and therefore k) . Then, using that f i, j (R, ε) is nondecreasing withrespect to R, the relations (3.10) and (3.7) yield 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. Then, for all ε 0 > 0, there exists R > 0 only depending on ε 0 , T , V 1 , V 2 , A, CardE, ρ 0 , µ and θ, such that 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 (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 which leads, using (3.14), to q 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 We thus get ( 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 ) we can then write, using (3.11), ( 3.18) Thus, we get from (3.17) and (3.18) (recall that B = A CardE) Since we get a bound of the second term at the right-hand side of (3.16), using we thus get from (3.16) and (3.19) (3.20) Let us now remark that, thanks to the relation (3.12), we obtain This implies using (2.9) and (2.7) Hence, we get from (3.20) 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 Then an easy inductive argument shows that, for all n = 0, . . . , N , we haveq We thus obtain 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.

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

LEMMA 3.5 Under assumptions (H), let
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 which means precisely that lim δt→0,|M| θδt (T M,δt 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 We then have and therefore there exists δt 1 > 0 only depending on ε, f , B and T such that if δt δt 1 , then We then get that . Therefore, there exists δt 2 , only depending on ε and f , such that |T M,δt We then get that, for a given R > 0, T M,δt and Since 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 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, This proves that |T M,δt,R
Proof. We consider the function f defined by which yields that 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 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 We thus get that 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).
We now observe that, thanks to (3.36), 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, 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 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 (3.37) The statement is then the following.
REMARK 3.3 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 .
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). If we compare the use of the finite-volume method for the approximation of the marginal distributions, and Monte Carlo simulations, we can make the following observations: • For both methods, it is possible to adjust the precision, by increasing the number of histories for Monte Carlo simulations and by refining the mesh for the finite-volume method. In the above example (see Labeau (1996)), the computing time is here about three quarters of an hour on a PC, which is longer than those given in Labeau (1996). The precision in Labeau (1996) is, however, lower (this is shown by a comparison between the marginal distributions).
• The finite-volume method has been shown to be efficient in some cases where some other methods have been unsuccessful (see Cocozza-Thivent & Eymard and Cocozza-Thivent & Eymard (2003)), with relatively short computing times on 1D problems. However, one can wonder if these computing times will remain acceptable if the dimension of the problem increases. On the contrary, it seems that Monte Carlo simulation is the method of choice for higher dimensional problems. • 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.