Ignacio Araya
email: iaraya@inf.utfsm.cl

Bertrand Neveu
email: bertrand.neveu@enpc.fr

Gilles Trombettoni
email: gilles.trombettoni@inria.fr

An Interval Extension Based on Occurrence Grouping

Keywords: Intervals, Interval extension, Monotonicity, Occurrence grouping I. Araya, UTFSM, Valparaiso, Chile

In interval arithmetics, special care has been brought to the definition of interval extension functions that compute narrow interval images. In particular, when a function f is monotonic w.r.t. a variable in a given domain, it is well-known that the monotonicity-based interval extension of f computes a sharper image than the natural interval extension does.

This paper presents a so-called "occurrence grouping" interval extension [ f ] og of a function f . When f is not monotonic w.r.t. a variable x in a given domain, we try to transform f into a new function f og that is monotonic w.r.t. two subsets x a and x b of the occurrences of x: f og is increasing w.r.t. x a and decreasing w.r.t. x b .

[ f ] og is the interval extension by monotonicity of f og and produces a sharper interval image than the natural extension does. For finding a good occurrence grouping, we propose a linear program and an algorithm that minimize a Taylor-based over-estimate of the image diameter of [ f ] og . Experiments show the benefits of this new interval extension for solving systems of nonlinear equations.

Introduction

The computation of sharp interval image enclosures is in the heart of interval arithmetics [START_REF] Moore | Interval Analysis[END_REF]. It allows a computer to evaluate a mathematical formula while taking into account in a reliable way round-off errors due to floating point arithmetics. Sharp enclosures also allow combinatorial interval methods to quickly converge towards the solutions of a system of constraints over the reals. At every node of the search tree, a test of existence checks that, for every equation f (X) = 0, the interval extension of f returns an interval including 0 (otherwise the branch is cut). Also, constraint propagation algorithms, used at every node of the search tree to reduce the search space, can be improved when they use better interval extensions. For instance, the Box constraint propagation algorithm [START_REF] Benhamou | Revising Hull and Box Consistency[END_REF] uses a test of existence inside its iterative splitting process.

This paper proposes a new interval extension and we first recall basic material about interval arithmetics [START_REF] Moore | Interval Analysis[END_REF][START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF][START_REF] Jaulin | Applied Interval Analysis[END_REF] to introduce the interval extensions useful in our work.

An (Vectorial variables appear in upper case in this article.)

An interval function [ f ] is a function from IR to IR n , IR being the set of all the intervals over R. When a function f is a composition of elementary functions, an extension of f to intervals must be defined to ensure a conservative image computation.

Definition 1 (Extension of a function to IR; also called inclusion function) Consider a function f : R n → R.

[ f ] : IR n → IR is an extension of f to intervals iff:

∀[V ] ∈ IR n [ f ]([V ]) ⊇ ℑ f ([V ]) ≡ { f (V ), V ∈ [V ]} ∀V ∈ R n f (V ) = [ f ](V )
The first idea is to use interval arithmetics. Interval arithmetics extends to intervals arithmetic operators +, -, ×, / and elementary functions (power, exp, log, sin, cos, ... The optimal image [ f ] opt ([V ]) is the sharpest interval containing ℑ f ([V ]). If f is continuous inside a box [V ], the natural evaluation of f (i.e., the computation of [ f ] n ([V ])) yields the optimal image when each variable occurs only once in f . When a variable appears several times in f , the evaluation by interval arithmetics generally produces an over-estimate of [ f ] opt ([V ]), because the correlation between the occurrences of a same variable is lost. Two occurrences of a variable are handled as independent variables. For example

[x] -[x], with [x] = [0, 1] gives the result [-1, 1], instead of [0, 0], as does [x] -[y], with [x] = [0, 1] and [y] = [0, 1]
. Thus, multiple occurrences of variables render NP-hard the computation of the optimal image of a polynomial [START_REF] Kreinovich | Computational Complexity and Feasibility of Data Processing and Interval Computations[END_REF]. This main drawback of interval arithmetics causes a real difficulty for implementing efficient interval-based solvers. There exist many possible interval extensions for a function, the difficulty being to define an extension that computes a sharp approximation of the optimal image at a low cost.

The first-order Taylor extension [ f ] t of f , also called centered form [START_REF] Moore | Interval Analysis[END_REF], uses the Taylor form of f :

[ f ] t ([V ]) = f (V m ) + n ∑ i=1 ∂ f ∂ x i ([V ]) . ([x i ] -x m i )
where n is the number of variables in f , x m i is the value in the middle of the interval

[x i ], V m is the n-dimensional point in the middle of [V ] (i.e., V m = (x m 1 , ..., x m k )) and ∂ f ∂ x i is an interval extension of ∂ f ∂ x i .
The Taylor extension generally calculates sharp evaluations when the diameters of the partial derivatives are close to 0. In other cases it can be even worse than the natural extension.

A well-known variant of the Taylor extension, called here Hansen extension [ f ] h , computes a sharper image at a higher cost [START_REF] Hansen | Global Optimization using Interval Analysis[END_REF]:

[ f ] h ([V ]) ⊆ [ f ] t ([V ]).
Another extension to intervals uses the monotonicity of a function in a given domain. When f is monotonic w.r.t. a subset of variables, one can replace, in the natural evaluations, the intervals of these monotonic variables1 by degenerated intervals reduced to their maximal or minimal values [START_REF] Moore | Interval Analysis[END_REF][START_REF] Hansen | Global Optimization using Interval Analysis[END_REF].

Definition 2 ( f min , f max , monotonicity-based extension) Let f be a function defined on variables V of domains [V ]. Let X ⊆ V be a subset of monotonic variables. Consider the values x + i and x - i such that: if x i ∈ X is an increasing (resp. decreasing) variable, then x - i = x i and x + i = x i (resp. x - i = x i and x + i = x i ). Consider W = V \ X the set of variables not detected monotonic. Then, f min and f max are functions defined by:

f min (W ) = f (x - 1 , ..., x - n ,W ) f max (W ) = f (x + 1 , ..., x + n ,W ) Finally, the monotonicity-based extension [ f ] m of f in the box [V ] produces the following interval image: [ f ] m ([V ]) = [ f min ] n ([W ]), [ f max ] n ([W ])
The image

[ f ] m ([V ]
) obtained by the monotonicity-based extension is sharper than, or equal to, the image [ f ] n ([V ]) obtained by natural evaluation. That is:

[ f ] opt ([V ]) ⊆ [ f ] m ([V ]) ⊆ [ f ] n ([V ])
In addition, when a function is monotonic w.r.t. each of its variables, i.e., when W is empty in Def. 2, the problem of multiple occurrences disappears and the evaluation (using a monotonicitybased extension) becomes optimal:

[ f ] m ([V ]) = [ f ] opt ([V ]).
Note that the bounds of the evaluation by monotonicity can be computed using any interval extension, not necessarily the natural extension. When the bounds are computed with the Taylor (resp. Hansen) extension, we denote this variant by

[ f ] m+t (resp. [ f ] m+h ).
When the evaluation by monotonicity uses, recursively, the same evaluation by monotonicity for computing the bounds of the image, we call it recursive evaluation by monotonicity and denote it by [ f ] mr . This computes images sharper than or equal to the evaluation by monotonicity does (i.e.,

[ f ] mr ([V ]) ⊆ [ f ] m ([V ]
)) at a cost between 2 and 2n higher. The extended paper illustrates this interval extension [START_REF] Araya | An Interval Extension Based on Occurrence Grouping: Method and Properties[END_REF].

Contribution

This paper explains how to use monotonicity when a function is not monotonic w.r.t. a variable x, but is monotonic w.r.t. a subgroup of occurrences of x. We present in the next section the idea of grouping the occurrences into three sets, increasing, decreasing and non monotonic auxiliary variables. Linear programs for obtaining "interesting" occurrence groupings are described in Sections 3 and 4. In Section 5, we propose an algorithm to solve the linear programming problem presented in Section 4. Finally, in Section 7, experiments show that this new occurrence grouping interval extension function compares favorably with existing ones and show its benefits for solving systems of equations, in particular when we use a filtering algorithm like Mohc [START_REF] Araya | Exploiting Monotonicity in Interval Constraint Propagation[END_REF][START_REF] Chabert | Hull Consistency Under Monotonicity[END_REF] exploiting monotonicity.

Evaluation by monotonicity with occurrence grouping

In this section, we study the case of a function which is not monotonic w.r.t. a variable with multiple occurrences. We can, without loss of generality, limit the study to a function of one variable: the generalization to a function of several variables is straightforward, the evaluations by monotonicity being independent.

Example 1 Consider f 1 (x) = -x 3 + 2x 2 + 14x. We want to calculate a sharp evaluation of this function when x falls in [-2, 1]. The derivative of f 1 is f ′ 1 (x) = -3x 2 + 4x + 14 and contains a positive term [START_REF] Moore | Interval Analysis[END_REF], a negative term (-3x 2 ) and the term 4x that is negative when x ∈ [-2, 0] and positive when

x ∈ [0, 1]. [ f 1 ] opt ([V ]
) is [-13.18, 15], but we cannot obtain it directly by a simple interval function evaluation (one needs to solve f ′ 1 (x) = 0, which is difficult in the general case). In the interval [-2, 1], f 1 is not monotonic. The natural interval evaluation yields [-29, 30], the Horner evaluation [-34, 17] (see [START_REF] Horner | A new Method of Solving Numerical Equations of all Orders, by Continuous Approximation[END_REF]).

When a function is not monotonic w.r.t. a variable x, it sometimes appears that it is monotonic w.r.t. some occurrences. A first naive idea leads to replace the function f by a function f nog , grouping all increasing occurrences into one variable x a , all decreasing occurrences into one variable x b , and the non monotonic occurrences into x c . The domain of the new auxiliary variables is the same:

[x a ] = [x b ] = [x c ] = [x]
. However, the evaluation by monotonicity of the new function f nog always provides the same result as the natural evaluation. The main idea is then to change this grouping in order to reduce the dependency problem and obtain sharper evaluations. We can indeed group some occurrences (increasing, decreasing, or non monotonic) into an increasing variable x a (resp. a decreasing variable x b ) as long as the function remains increasing (resp. decreasing) w.r.t. this variable x a (resp. x b ). If one can move a non monotonic occurrence into a monotonic group, the evaluation will be better (or remain the same). Also, if it is possible to transfer all decreasing occurrences into the increasing part, the dependency problem will now occur only on the occurrences in the increasing and non monotonic parts.

For f 1 , if we group together the positive derivative term with the derivative term containing zero, we obtain the new function:

f og 1 (x a , x b ) = -x 3 b + 2x 2 a + 14x a .
As the interval derivative of the grouping of the first two occurrences (the variable x a ) is positive:

4[x] + 14 = [6, 18], f og
1 is increasing w.r.t. x a . We can then achieve the evaluation by monotonicity and obtain the interval [-21, 24]. We can in the same manner obtain f og 1 (x a , x c ) = -x 3 a + 2x 2 c + 14x a , the evaluation by monotonicity yields then [-20, 21]. We remark that we find sharper images than the natural evaluation of f 1 does.

Interval extension by occurrence grouping

Consider the function f (x) with multiple occurrences of x. We obtain a function f og (x a , x b , x c ) by replacing in f every occurrence of x by one of the three variables x a , x b , x c , such that f og is increasing w.r.t. x a in [x], and f og is decreasing w.r.t. x b in [x]. Then, we define the interval extension by occurrence grouping of f by:

[ f ] og ([V ]) := [ f og ] m ([V ])
Unlike the natural interval extension and the interval extension by monotonicity, the interval extension by occurrence grouping is not unique for a function f since it depends on the occurrence grouping (og) that transforms f into f og .

A 0,1 linear program to perform occurrence grouping

In this section, we propose a method for automatizing occurrence grouping. Using the Taylor extension, we first compute an over-estimate of the diameter of the image computed by [ f ] og . Then, we propose a linear program performing a grouping that minimizes this over-estimate.

Taylor-based over-estimate

First, as f og could be not monotonic w.r.t. x c , the evaluation by monotonicity considers the occurrences of x c as different variables such as the natural evaluation would do.

Proposition 1 ([14]) Let f (x) be a continuous function in a box

[V ] with a set of occur- rences of x: {x 1 , x 2 , ..., x k }. f • (x 1 , .., x k ) is a function obtained from f by considering all the occurrences of x as different variables. Then, [ f ] n ([V ]) = [ f • ] opt ([V ]).
Second, as f og is monotonic w.r.t. x a and x b , the evaluation by monotonicity of these variables is optimal.

Proposition 2 ([14]) Let f (x 1 , ..., x n ) be a monotonic function w.r.t. each of its variables in a box [V ] = {[x 1 ], ..., [x n ]}. Then, the evaluation by monotonicity computes [ f ] opt ([V ]).
Using these propositions, we observe that

[ f og ] m ([x a ], [x b ], [x c ]) is equivalent to [ f • ] opt ([x a ], [x b ], [x c 1 ], ..., [x c ck ]
), considering each occurrence of x c in f og as an independent variable x c j in f • , ck being the number of occurrences of x c in f og . Using the Taylor evaluation, an upper bound of diam([ f ] opt ([V ])) is given by the right side of (1) in Proposition 3.

Proposition 3 Let f (x 1 , ..., x n ) be a function with domains [V ] = {[x 1 ], ..., [x n ]}. Then, diam([ f ] opt ([V ])) ≤ n ∑ i=1 diam([x i ]) . |[g i ]([V ])| (1)
where

[g i ] is an interval extension of g i = ∂ f ∂ x i .
Using Proposition 3, we can calculate an upper bound of the diameter of

[ f ] og ([V ]) = [ f og ] m ([V ]) = [ f • ] opt ([V ]): diam([ f ] og ([V ])) ≤ diam([x]) |[g a ]([V ])| + |[g b ]([V ])| + ck ∑ i=1 |[g c i ]([V ])| where [g a ], [g b ] and [g c i ] are the interval extensions of g a = ∂ f og ∂ x a , g b = ∂ f og ∂ x b and g c i = ∂ f og ∂ x c i respectively. diam([x]) is factorized because [x] = [x a ] = [x b ] = [x c 1 ] = ... = [x c ck ].
In order to respect the monotonicity conditions required by

f og : ∂ f og ∂ x a ([V ]) ≥ 0, ∂ f og ∂ x b ([V ]) ≤ 0, we have the sufficient conditions [g a ]([V ]) ≥ 0 and [g b ]([V ]) ≤ 0, implying |[g a ]([V ])| = [g a ]([V ]) and |[g b ]([V ])| = -[g b ]([V ]). Finally: diam([ f ] og ([V ])) ≤ diam([x]) [g a ]([V ]) -[g b ]([V ]) + ck ∑ i=1 |[g c i ]([V ])| (2)

A linear program

We want to transform f into a new function f og that minimizes the right side of the relation [START_REF] Araya | Exploiting Monotonicity in Interval Constraint Propagation[END_REF]. The problem can be easily transformed into the following integer linear program: Find the values r a i , r b i and r c i for each occurrence x i that minimize

G = [g a ]([V ]) -[g b ]([V ]) + k ∑ i=1 |[g i ]([V ])|r c i (3) subject to: [g a ]([V ]) ≥ 0 (4) [g b ]([V ]) ≤ 0 (5)
r a i , r b i , r c i ∈ {0, 1}; r a i + r b i + r c i = 1 for i = 1, ..., k, (6)
where a value r a i , r b i or r c i equal to 1 indicates that the occurrence x i in f will be replaced, respectively, by x a ,

x b or x c in f og . k is the number of occurrences of x, [g a ]([V ]) = k ∑ i=1 [g i ]([V ])r a i , [g b ]([V ]) = k ∑ i=1 [g i ]([V ])r b i , and [g 1 ]([V ]), ..., [g k ]([V ])

Linear program corresponding to Example 1

We have [START_REF] Moore | Interval Analysis[END_REF][START_REF] Moore | Interval Analysis[END_REF]. Then, the linear program is: Find the values r a i , r b i and r c i that minimize

f 1 (x) = -x 3 + 2x 2 + 14x, f ′ 1 (x) = -3x 2 + 4x + 14 with x ∈ [-2, 1]. The gradient is: [g 1 ]([-2, 1]) = [-12, 0], [g 2 ]([-2, 1]) = [-8, 4] and [g 3 ]([-2, 1]) =
G = 3 ∑ i=1 [g i ]([V ])r a i - 3 ∑ i=1 [g i ]([V ])r b i + 3 ∑ i=1 |[g i ]([V ])|r c i = (4r a 2 + 14r a 3 ) + (12r b 1 + 8r b 2 -14r b 3 ) + (12r c 1 + 8r c 2 + 14r c 3 ) subject to: 3 ∑ i=1 [g i ]([V ])r a i = -12r a 1 -8r a 2 + 14r a 3 ≥ 0; 3 ∑ i=1 [g i ]([V ])r b i = 4r b 2 + 14r b 3 ≤ 0 r a i , r b i , r c i ∈ {0, 1}; r a i + r b i + r c i = 1 for i = 1, ..., 3
We obtain the minimum 22, and the solution

r a 1 = 1, r b 1 = 0, r c 1 = 0, r a 2 = 0, r b 2 = 0, r c 2 = 1, r a 3 = 1, r b 3 = 0, r c 3 = 0,
which is the last solution presented in Section 2. Note that the value of the over-estimate of diam(

[ f ] og ([V ])) is equal to 66 (22 * diam[-2, 1]) whereas diam([ f ] og ([V ]
)) = 41. Although the over-estimate is rough, the heuristic works rather well on this example. Indeed, diam(

[ f ] n ([V ])) = 59 and diam([ f ] opt ([V ])) = 28.18.

A linear programming problem achieving a better occurrence grouping

The linear program above is a 0,1 linear program and is known to be NP-hard in general. We can render it tractable while, more important in practice, improving the minimum G by allowing r a i , r b i and r c i to get real values. In other words, we allow each occurrence of x in f to be replaced by a convex linear combination of auxiliary variables, x a , x b and x c such that f og is increasing w.r.t. x a and decreasing w.r.t.

x b in [x].
Definition 3 (Interval extension by occurrence grouping) Let f (x) be a function with multiple occurrences of the variable x. f og (x a , x b , x c ) is the function obtained by replacing in f every occurrence of x by r a i x a + r b i x b + r c i x c , such that:

-r a i , r b i , r c i ∈ [0, 1] 3 and r a i + r b i + r c i = 1, -∂ f og ∂ x a ([x], [x], [x]) ≥ 0 and ∂ f og ∂ x b ([x], [x], [x]) ≤ 0.

The interval extension by occurrence grouping of f is defined by

[ f ] og ([x]) := [ f og ] m ([x], [x], [x])
Note that f and f og have the same natural evaluation.

In Example 1, we can replace f 1 by f og 1 or f og 2 in a way respecting the monotonicity constraints of x a and x b : 1.

1. f og 1 1 (x a , x b ) = -( 1 2 x a + 1 2 x b )3 + 2x a 2 + 14x a → [ f og 1 1 ] m ([-2, 1]) = [-19.875, 16.125] 2. f og 2 1 (x a , x b , x c ) = -x 3 a +2(0.25x a +0.75x c ) 2 +14x a → [ f og 2 1 ] m ([-2, 1]) = [-
f og 1 2 (x a , x b ) = x 3 a -( 1 4 x a + 3 4 x b ) → [ f og 1 2 ] m ([x]) = [-0.75, 6.375] 2. f og 2 2 (x a , x b ) = ( 11 12 x a + 1 12 x b ) 3 -x b → [ f og 2 2 ] m ([x]) = [-1.756, 6.09] Algorithm 1 Occurrence Grouping (in: f ,[g * ] out: f og ) 1: [G 0 ] ← k ∑ i=1 [g i ] 2: [G m ] ← ∑ 0 ∈[g i ] [g i ] 3: if 0 ∈ [G 0 ] then 4: OG case1([g * ], [r a * ], [r b * ], [r c * ]) 5: else if 0 ∈ [G m ] then 6: OG case2([g * ], [r a * ], [r b * ], [r c * ]) 7: else 8: /* 0 ∈ [G m ] and 0 ∈ [G 0 ] */ 9: if G m ≥ 0 then 10: OG case3 + ([g * ], [r a * ], [r b * ], [r c * ]) 11: else 12: OG case3 -([g * ], [r a * ], [r b * ], [r c * ]) 13: end if 14: end if 15: f og ← Generate New Function( f , [r a * ], [r b * ], [r c * ]) Algorithm 2 OG case2 (in: [g * ] out: [r a * ], [r b * ], [r c * ]) 1: [G + ] ← ∑ [g i ]≥0 [g i ] 2: [G -] ← ∑ [g i ]≤0 [g i ] 3: [α 1 ] ← G + G -+ G -G - G + G --G -G + 4: [α 2 ] ← G + G + + G -G + G + G --G -G + 5: 6: for all [g i ] ∈ [g * ] do 7: if g i ≥ 0 then 8: ([r a i ], [r b i ], [r c i ]) ← (1 -[α 1 ], [α 1 ], 0) 9: else if g i ≤ 0 then 10: ([r a i ], [r b i ], [r c i ]) ← ([α 2 ], 1 -[α 2 ], 0) 11: else 12: ([r a i ], [r b i ], [r c i ]) ← (0, 0, 1) 13:
end if 14: end for Thus, the new linear program that computes convex linear combinations for achieving occurrence grouping becomes: Find the values r a i , r b i and r c i for each occurrence x i that minimize (3) subject to (4), ( 5), [START_REF] Chabert | Hull Consistency Under Monotonicity[END_REF] and

r a i , r b i , r c i ∈ [0, 1] for i = 1, ..., k. (7)
Note that this continuous linear program improves the minimum of the objective function because the integrity conditions are relaxed.

Examples

In Example 1, we obtain the minimum 21 and the new function The image [-0.75, 6.01] obtained by occurrence grouping is sharper than the interval computed by natural and Horner evaluations. In this case, the 0,1 linear program of Section 3 yields the naive grouping.

f og 1 (x a , x b , x c ) = -x 3 a + 2(0.25x a + 0.75x c ) 2 + 14x a : [ f og 1 ] m ([x]) = [-

An efficient Occurrence Grouping algorithm

Algorithm 1 finds r a i , r b i , r c i that minimize G subject to the constraints. At line 15, the algorithm generates symbolically the new function f og that replaces each occurrence

x i in f by [r a i ]x a + [r b i ]x b + [r c i ]x c .
Note that the values are represented by thin intervals, of a few u.l.p. large, for taking into account the floating point rounding errors appearing in the computations. Algorithm 1 uses a vector

[g * ] of size k containing interval derivatives of f w.r.t. each occurrence x i of x. Each component of [g * ] is denoted by [g i ] and corresponds to the interval ∂ f ∂ x i ([V ]).
A symbol indexed by an asterisk refers to a vector (e.g., [g * ], [r a * ]). [START_REF] Araya | An Interval Extension Based on Occurrence Grouping: Method and Properties[END_REF] 9: while g a + g i ≥ 0 do 10:

Algorithm 3 OG case3 + (in:[g * ] out:[r a * ], [r b * ], [r c * ]) 1: [g a ] ← [0, 0] 2: for all [g i ] ∈ [g * ], g i ≥ 0 or g i ≤ 0 do 3: [g a ] ← [g a ] + [g i ] /*
: index ← descending sort ( {[g i ] ∈ [g * ], g i < 0}, criterion → g i -|[g i ]| g i ) 8: j ← 1 ; i ← index
([r a i ] ], [r b i ], [r c i ]) ← (1, 0, 0) 11: [g a ] ← [g a ] + [g i ] 12: j ← j + 1 ; i ← index[ j] 13: end while 14: 15: [α] ← - ga g i 16: ([r a i ], [r b i ], [r c i ]) ← ([α], 0, 1 -[α]) /* [g a ] ← [g a ] + [α][g i ] */ 17: 18: j ← j + 1 ; i ← index[ j] 19: while j ≤ length(index) do 20: ([r a i ], [r b i ], [r c i ]) ← (0, 0, 1) 21: j ← j + 1 ; i ← index[ j] 22: end while
At line 1, the partial derivative [G 0 ] of f w.r.t. x is calculated using the sum of the partial derivatives of f w.r.t. each occurrence of x. At line 2, [G m ] gets the value of the partial derivative of f w.r.t. the monotonic occurrences of x.

According to the values of [G 0 ] and [G m ], we can distinguish 3 cases. The first case is well-known (0 ∈ [G 0 ] in line 3) and occurs when x is a monotonic variable. In the procedure OG case1, all the occurrences of x are replaced by x a (if

[G 0 ] ≥ 0) or by x b (if [G 0 ] ≤ 0)
. The evaluation by monotonicity of f og is equivalent to the evaluation by monotonicity of f . In the second case, when 0 ∈ [G m ] (line 5), the procedure OG case2 (Algorithm 2) achieves a grouping of the occurrences of x. Increasing occurrences are replaced by (1α 1 )x a + α 1 x b , decreasing occurrences by α 2 x a + (1 -α 2 )x b and non monotonic occurrences by x c (lines 7 to 13 of Algorithm 2). The third case occurs when 0 ∈ [G m ] and 0 ∈ [G 0 ]. W.l.o.g., assume that G m ≥ 0. The procedure OG case3 + (Algorithm 3) first groups all the positive and negative occurrences in the increasing group x a (lines 2-5). The non monotonic occurrences are then replaced by x a in an order determined by an array index 2 (line 7) as long as the constraint k ∑ i=1 r a i g i ≥ 0 is satisfied (lines 9-13). The first occurrence x i ′ that cannot be (entirely) replaced by x a because it would make the constraint (4) unsatisfiable is replaced by αx a + (1 -α)x c , with α such that the constraint is satisfied and equal to 0, i.e., [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF][START_REF] Trombettoni | Constructive Interval Disjunction[END_REF]. The rest of the occurrences are replaced by x c (lines 18-22). The extended paper [START_REF] Araya | An Interval Extension Based on Occurrence Grouping: Method and Properties[END_REF] illustrates how the procedure OG case3 + (resp. OG case2) handles the example f 1 (resp. f 2 ). 2 An occurrence x i 1 is handled before

( k ∑ i=1,i =i ′ r a i g i ) + αg i ′ = 0 (lines
x i 2 if g i1 -|[g i1 ]| g i1 ≥ g i2 -|[g i2 ]| g i2
. index[ j] yields the index of the j th occurrence in this order.

Properties

Proposition 4 Algorithm 1 (Occurrence grouping) is correct and solves the linear program that minimizes (3), modulo floating-point roundings, subject to the constraints ( 4), ( 5), ( 6) and [START_REF] Flannery | Numerical Recipes in C[END_REF].

We can check that Algorithm 1 respects the four constraints ( 4)- [START_REF] Flannery | Numerical Recipes in C[END_REF]. We have also proven that the minimum of the objective function ( 3) is reached. The proof concerning OG case3 is sophisticated, due to the sort of indices, and uses known results about the continuous knapsack problem. Special care has been brought to ensure the correctness modulo floatingpoint roundings. A full proof can be found in [START_REF] Araya | An Interval Extension Based on Occurrence Grouping: Method and Properties[END_REF].

Proposition 5 The time complexity of Occurrence Grouping for one variable with k occurrences is O(k log 2 (k)). It is time O(n k log 2 (k)) when a multi-variate function is iteratively transformed by Occurrence Grouping for each of its n variables having at most k occurrences each. A preliminary gradient calculation by automatic differentiation is time O(e)

, where e is the number of unary and binary operators in the expression.

Computing the gradient of a function amounts in two traversals of the tree representing the mathematical expression [START_REF] Benhamou | Revising Hull and Box Consistency[END_REF]. The time complexity of Algorithm 1 is dominated by that of descending sort in the OG case3 procedure.

Instead of Algorithm 1, we may use a standard Simplex algorithm providing that the used Simplex implementation takes into account floating-point rounding errors. A performance comparison between Algorithm 1 and Simplex is shown in Section 7.3. Also, as shown in Section 7.2, the time required in practice by Occurrence Grouping is negligible when it is used for solving systems of equations.

Although Occurrence Grouping can be viewed as a heuristic since it minimizes a Taylor-based over-estimate of the function image diameter, it is important to stress that our new interval extension improves the well-known monotonicity-based interval extension.

Proposition 6 Consider a function f : R n → R, and the previously defined interval natural ([ f ] n ), monotonicity-based ([ f ] m ) and occurrence grouping ([ f ] og ) extensions of f . Let V be the n variables involved in f with domains

[V ]. Then, [ f ] og ([V ]) ⊆ [ f ] m ([V ]) ⊆ [ f ] n ([V ]).

Experiments

Occurrence Grouping has been implemented in the Ibex [START_REF] Chabert | Contractor Programming[END_REF][START_REF] Chabert | Ibex -An Interval Based EXplorer[END_REF] open source interval-based solver in C++. The main goal of these experiments is to show the improvements in CPU time brought by Occurrence Grouping when solving systems of equations.

We briefly recall the combinatorial process followed by an interval-based solver to find all the solutions of a system of equations. The solving process starts from an initial box representing the search space and builds a search tree, following a Branch & Contract scheme:

-Branch: the current box is bisected on one dimension, generating two sub-boxes.

-Contract: filtering (also called contraction) algorithms reduce the bounds of the box with no loss of solution.

The process terminates with boxes of size smaller than a given positive ω.

Contraction algorithms comprise multidimensional interval Newton algorithms (or variants) issued from the numerical interval analysis community [START_REF] Moore | Interval Analysis[END_REF][START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF] along with algorithms from constraint programming. In all our experiments, at each node of the search tree, i.e., between two bisections, the solving strategy uses two types of contractors (all available in Ibex) in sequence:

1. (from constraint programming) the shaving/slicing contractor 3BCID [START_REF] Lhomme | Consistency Techniques for Numeric CSPs[END_REF][START_REF] Trombettoni | Constructive Interval Disjunction[END_REF] and a recent constraint propagation algorithm, Mohc [START_REF] Araya | Exploiting Monotonicity in Interval Constraint Propagation[END_REF][START_REF] Chabert | Hull Consistency Under Monotonicity[END_REF], exploiting monotonicity of functions; 2. (from interval analysis) an interval Newton using a Hansen-Sengupta matrix, a leftpreconditioning of the matrix, and a Gauss-Seidel method to solve the interval linear system [START_REF] Hansen | Global Optimization using Interval Analysis[END_REF].

Sixteen systems of equations have been used in our experiments. They are issued from the COPRIN team website [START_REF] Merlet | The COPRIN examples page[END_REF], most of them being also known in the COCONUT benchmark suite devoted to interval analysis and global optimization. 3 They correspond to square systems with a finite number of zero-dimensional (isolated) solutions of at least two constraints involving multiple occurrences of variables and requiring more than 1 second to be solved (considering the times appearing in the COPRIN website). All experiments have been performed on a same computer (Intel 6600 2.4 GHz).

Comparison between interval extensions

We first report a comparison between the evaluation by occurrence grouping (i.e., the diameter of [ f ] og ([V ])) and several existing interval evaluations ([ f ] ext ), including Taylor, Hansen and monotonicity-based extensions (see Section 1). Since the Taylor and Hansen extensions are not comparable with the natural extension, to obtain more reasonable comparisons, we have redefined

[ f ] t ([V ]) = [ f ] ′ t ([V ]) ∩ [ f ]([V ]) and [ f ] h ([V ]) = [ f ] ′ h ([V ]) ∩ [ f ] n ([V ]), where [ f ] ′
t and [ f ] ′ h are the actual Taylor and Hansen extensions respectively. The first column of Table 1 indicates the name of each instance. The other columns refer to existing extensions [ f ] ext and report an average of ratios

System [ f ] n [ f ] t [ f ] h [ f ] m [ f ] mr [ f ] mr+h [ f ] mr+og Brent 0.
ρ ext = Diam([ f ] og ([V ])) Diam([ f ] ext ([V ]))
. These ratios are measured while the Branch & Contract solving strategy mentioned above4 is run to solve the tested system of equations. They are calculated every time a constraint is handled inside the constraint propagation, at every node of the search tree, thus avoiding biases.

The table highlights that [ f ] og computes, in general, sharper interval images than all the competitors. (Only [ f ] mr+og achieves better evaluations, but it also uses occurrence grouping.) The improvements w.r.t. the two evaluation methods by monotonicity (i.e., [ f ] m and [ f ] mr ) corroborate the benefits of our approach. For example, in Fourbar, [ f ] og obtains an evaluation diameter which is 42.7% of the evaluation diameter provided by [ f ] mr .

[ f ] mr+h obtains the sharpest evaluations in three benchmarks (Caprasse, Fourbar and Virasoro). However, [ f ] mr+h is more expensive than [ f ] og . [ f ] mr+h requires computing 2n interval partial derivatives, thus traversing 4n times the expression tree if an automatic differentiation method is used [START_REF] Hansen | Global Optimization using Interval Analysis[END_REF].

The extension using occurrence grouping [ f ] mr+og provides necessarily a better evaluation than, or equal to, [ f ] og . However, the experiments on the tested systems show that the gain in evaluation diameter is only 1.6% on average (between 0% and 6.2%), so that we do not believe it constitutes a promising extension due to its additional cost.

Occurrence Grouping inside a monotonicity-based contractor

These experiments are significant in that they underline the benefits of occurrence grouping for improving the solving of systems of constraints. Mohc [START_REF] Araya | Exploiting Monotonicity in Interval Constraint Propagation[END_REF][START_REF] Chabert | Hull Consistency Under Monotonicity[END_REF] is a new constraint propagation contractor (like HC4 or Box [START_REF] Benhamou | Revising Hull and Box Consistency[END_REF]) that exploits the monotonicity of a function to improve the contraction of the related variable intervals.

Table 2 shows the results of Mohc without the OG algorithm (¬OG), and with Occurrence-Grouping (OG), i.e., when the function f is transformed into f og before applying the main MohcRevise( f og ) procedure. We observe that, for 7 of the 16 benchmarks, Occurrence-Grouping is able to improve significantly the results of Mohc; in Butcher, Fourbar, Virasoro and Yamamura1 the gains in CPU time ( ¬OG OG ) obtained are greater than 30, 11, 7.5 and 5.4 respectively.

Practical time complexity

We have first compared the performance of two Occurrence Grouping implementations: using our ad-hoc algorithm (Occurrence Grouping) and using a Simplex method. The basic Simplex implementation [START_REF] Flannery | Numerical Recipes in C[END_REF] we have used is not rigorous, i.e., it does not take into account rounding errors due to floating point arithmetic. Adding this feature should make the algorithm run more slowly. Two important results have been obtained. First, we have checked experimentally that our algorithm is correct, i.e., it obtains the minimum value for the objective function G. Second, just as we expected, the performance of the generalpurpose Simplex method is worse than the performance of our algorithm. It runs between 2.32 (Brent) and 10 (Virasoro) times more slowly.

Two results highlight that the CPU time required by Occurrence Grouping is very interesting in practice. In Table 7.2 indeed, for instances like Brent, Eco9 or Trigexp2, Occurrence Grouping is called a large number of times with roughly no effect on solving: similar number of choice points is reported with and without occurrence grouping. However, the overall CPU is also very similar. Another experiment is reported in [START_REF] Araya | An Interval Extension Based on Occurrence Grouping: Method and Properties[END_REF].

Table 2 Experimental results using the monotonicity-based contractor Mohc inside a solving strategy. The first column indicates the name of each instance, along with its number of variables/equations (left) and solutions (right). The second column reports the CPU time (first row of a multi-row) and the number of nodes (second row) obtained by a strategy using 3BCID(Mohc) without OG. The third column report the results of our strategy using 3BCID(OG+Mohc). The fourth column indicates the (large) number of calls to Occurrence Grouping during the solving process.

System

Conclusion

We have proposed a new method to improve the monotonicity-based evaluation of a function f . This Occurrence Grouping method creates for each variable of f three auxiliary, respectively increasing, decreasing and non monotonic variables. It then transforms f into a function f og that replaces the occurrences of every variable by a convex linear combination of these auxiliary variables. The evaluation by monotonicity of f og defines the evaluation by occurrence grouping of f and is better than the evaluation by monotonicity of f .

The extension by occurrence grouping computes at a low cost sharper interval images than existing interval extensions do. The main benefits of occurrence grouping lie in the improvement of an efficient contraction algorithm, called Mohc, that exploits monotonicity of functions. Occurrence grouping transforms, with nearly no overhead in practice, the constraints on the fly, during the constraint propagation and the solving process.

  interval [x] = [a, b] is the set of real numbers between a and b. x = a denotes the minimum of [x] and x = b denotes the maximum of [x]. The diameter/width of an interval is: diam([x]) = xx, and the absolute value of an interval is: |[x]| = max(|x|, |x|). A Cartesian product of intervals is named a box, and is denoted by a vector. [V ] = {[x 1 ], [x 2 ], ..., [x n ]}.

  ). For instance, [a, b] + [c, d] = [a + c, b + d]. The natural interval extension [ f ] n of a real function f simply replaces arithmetic over the reals by interval arithmetic.

Example 2

 2 20, 16.125] Consider the function f 2 (x) = x 3x and the interval [x] = [0.5, 2]. f 2 is not monotonic and the optimal image [ f 2 ] opt ([x]) is [-0.385, 6]. The natural evaluation yields [-1.975, 7.5], the Horner evaluation [-1.5, 6]. We can replace f 2 by one of the following functions (among others).

  20, 16.25]. The minimum 21 is inferior to 22 (obtained by the 0,1 linear program). The evaluation by occurrence grouping of f 1 yields [-20, 16.25], which is sharper than the image [-20, 21] obtained by the 0,1 linear program presented in Section 3. In Example 2, we obtain the minimum 11.25 and the new function f og 2 (x a , x b ) = ( 44 45 x a + 1 45 x b ) 3 -( 11 15 x a + 4 15 x b ).

  All positive and negative derivatives are absorbed by [g a ] */ 4: ([r a i ], [r b i ], [r c i ]) ← (1, 0, 0)

	5: end for
	6:
	7

Table 1

 1 Different interval extensions compared to [ f ] og

For the sake of conciseness, we sometimes write that a "variable x is monotonic" instead of writing that f is monotonic w.r.t. x.

See www.mat.univie.ac.at/ ˜neum/glopt/coconut/Benchmark/Benchmark.html

Similar results are obtained by other strategies.


Ignacio Araya
email: iaraya@inf.utfsm.cl

Gilles Trombettoni
email: gilles.trombettoni@inria.fr

Bertrand Neveu
email: bertrand.neveu@enpc.fr

Gilles Chabert
email: gilles.chabert@emn.fr

Upper Bounding in Inner Regions for Global Optimization under Inequality Constraints

Keywords: global optimization, upper bounding, intervals, branch and bound, inner regions

In deterministic constrained global optimization, upper bounding the objective function generally resorts to local minimization at the nodes of the branch and bound. The local minimization process is sometimes costly when constraints must be respected.

We propose in this paper an alternative approach when the constraints are inequalities or relaxed equalities so that the feasible space has a non-null volume. First, we extract an inner region, i.e., an (entirely feasible) convex polyhedron or box in which all points satisfy the constraints. Second, we select a point inside the extracted inner region and update the upper bound with its cost.

. This upper bounding shows good performance in medium-sized systems proposed in the COCONUT suite.

Upper bounding in inner regions

In deterministic constrained global optimization, upper bounding the objective function consists in nding a feasible point that improves the best cost already found in the branch and bound. Most global optimizers resort to local minimization1 using a Lagrangian relaxation. The considered function is sometimes big, which may render the local minimization slow. This paper describes an alternative approach for global optimization under inequality constraints dened by: min x∈[x] f (x) subject to g(x) ≤ 0, where f : R n → R is the real-valued objective function and g : R n → R m is a vector-valued function. x = {x 1 , ..., x i , ...x n } is a vector of variables varying in a box [x]. 2 x is said to be feasible if it satises the constraints.

The main idea is to exploit so called inner regions, i.e., subsets of the search space in which all points are feasible. Denition 1. Consider a system (f, g, x, [x]) comprising only inequality constraints. An inner region r in is a feasible subset of [x], i.e., r in ⊂ [x] and all points x ∈ r in satisfy g(x) ≤ 0.

At every node (iteration) of our interval B&B named IbexOpt [START_REF] Trombettoni | Inner Regions and Interval Linearizations for Global Optimization[END_REF], the cost is bounded above by using two inner region extraction algorithms, called InHC4 and InnerPolytope. InHC4 is described in Section 2. It tries to extract an inner box from the current outer box. If it fails, one simply picks a point randomly inside the outer box and checks its feasibility. If it succeeds, a simple monotonicity analysis of f replaces the intervals of the monotonic variables by the adequate bounds in the found inner box and the other values are randomly chosen. Then, InnerPolytope [START_REF] Trombettoni | Inner Regions and Interval Linearizations for Global Optimization[END_REF] builds a hyperplane for every inequality constraint. The hyperplane is produced by a special convex form of interval Taylor where the expansion point is chosen at a corner of the studied outer box. If it succeeds in building an inner polytope, the point minimizing the linearized form of the objective function is used to update the upper bound.

Contribution and limits

Contrarily to existing approaches, the proposed inner region extraction algorithms separate the feasibility part (handled rst, by inner region extraction) and the computation of the cost (handled next, inside the found inner region).

It is important to highlight that, like the other inner region extraction algorithms, ours can fail in nding an inner region even if one such region exists. However, they are rather inexpensive. In particular, InHC4 is faster when it fails in nding an inner region for a given constraint because the loop on all the constraints can be prematurely stopped (see below).

This upper bounding based on inner region extraction could also apply to thick and relaxed equations that dene a feasible space with a non-null volume. A thick equation has at least one coecient that can be modeled by an interval constant. This parameter corresponds to a bounded uncertainty, e.g., an imprecision on a measurement, or an irrational constant, like π. A pure equality f k (x) = 0 can also be handled with a relaxation as a thick equation f k (x) ∈ [-ǫ eq , +ǫ eq ], i.e, two inequalities -ǫ eq ≤ f k (x) ≤ ǫ eq . In this case of course, we can only guarantee the global optimum of the relaxed system, but ǫ eq can often be chosen almost arbitrarily small.

The InHC4 inner box extraction algorithm

InHC4 follows the simple and general scheme proposed in [START_REF] Benhamou | Universally Quantied Interval Constraints[END_REF][START_REF] Collavizza | Comparing Partial Consistencies[END_REF]. A main loop handles every constraint once in sequence and intersects incrementally the dierent boxes built.

3

The handling of the j th constraint uses as input the inner box returned by the handling of constraint g j-1 (x) ≤ 0. This box is inner w.r.t. the rst j -1 contraints. Handling the rst constraint g 1 (x) ≤ 0 is achieved with the outer box.

Handling the constraint g j (x) ≤ 0 consists in nding, inside the input box under construction, a box which is inner w.r.t. this single constraint. Bottom-up evaluation (see Fig. 1left): The tree is traversed from the leaves to the root and intervals associated to an operator are computed with interval arithmetics. For example, the node pointed by the arrow is initialized with the interval [0, 10]-[0, 15] = [-15, 10]. Thus, every node contains an interval corresponding to the natural interval evaluation of the subexpression. Top-down inner projection (see Fig. 1right): In each node related to a binary operator op and to an interval [z], the 2-dimensional box corresponding to its children x 1 and x 2 is reduced to an inner box

[x 1 ] in × [x 2 ] in such that: ∀(x 1 , x 2 ) ∈ [x 1 ] in × [x 2 ] in : x 1 op x 2 ∈ [z]
(1) If op corresponds to a unary operator, its unique child is reduced to [x] in , such that:

∀x ∈ [x] in : op(x) ∈ [z]
(2) If the inner projection returns an empty box (i.e., no box satises relation ( 1) or ( 2)), then the top-down process is interrupted. It means that InHC4R has failed in nding a box that is inner w.r.t. g j (x) ≤ 0. Since the approach is not complete because not all the feasible space is extracted during the top-down traversal, an inner box could be indeed missed by the process.

Consider the product operator of Fig. 1right and its two children. The reduced intervals appear in bold in the left side of each node. After the reduction of the product operator, its interval becomes [0, 5]. Before reduction, its children are associated to the intervals [10, 10] and [0, 1]. They are then reduced to [10, 10] and [0, 0.5] respectively. The reduction agrees with relation (1), i.e, ∀y ∈ [0, 0.5] : 10 * y ∈ [0, 5].

Inner projection for unary and binary basic operators

For unary, monotonic and continous operators, like log and exp, the inner projection is trivial and computes the (maximum) inner interval (i.e, no feasible point is lost, modulo roundos). It is very close to a standard projection in HC4R. However, for managing oating-point roundo errors, the outward rounding of HC4R is replaced by inward rounding.

For non monotonic unary operator like x 2 or sinus, a union of intervals is computed by HC4R, before returning the hull of these intervals. For an inner projection in InHC4R instead, only a single interval is kept since holes between these intervals contain inconsistent points.

For binary operators, Chabert and Beldiceanu in [START_REF] Chabert | Sweeping with Continuous Domains[END_REF] proposed inner projections, but with a case-by-case approach. We have extended their approach and built a more generic projection based on monotonicity properties. There usually exists an innite number of maximal boxes (as depicted in Fig. 2), and we have succeeded in designing inner projection operators that select randomly one maximal inner box. Note that these inner projections lead to heuristic choices since a single box cannot include the whole inner/feasible space. Also note that the two inequalities z ≤ (x 1 op x 2 ) ≤ z are handled in sequence, the inner box computed for one inequality being used as input of the second one.

For binary (or n-ary) operators that are monotonic w.r.t. each of their variables, a generic procedure, called MonoMaxInnerBox, can compute randomly one maximal inner box, if one such box exists, as shown in Fig. 2. This procedure is of course used for the addition and subtraction operators. It is also used for handling several (monotonic) subcases of the non monotonic binary operators: the multiplication and the division. Fig. 3 illustrates the two main cases for the multiplication x 1 * x 2 ∈ [z], depending whether 0 belongs or not to [z].

Handling the division operator amounts in rewriting

x 1 /x 2 = x 1 * 1 x 2 ∈ [z]
, although a direct implementation would also be possible.

Properties

We have proven that every implemented unary and binary operator computes a maximal inner box, modulo the loss involved by inward roundos. In case a constraint contains only a single occurrence of each variable, InHC4R thus computes a maximal inner box, when one such box is found. The result nally holds for a system of inequality constraints handled by InHC4.

Experiments

We have tested our original upper bounding procedure on a sample of about thirty constrained global optimization found in the COCONUT benchmark suite. 24 of them correspond to the most dicult systems selected by Ninin et al. [START_REF] Ninin | A Reliable Ane Relaxation Method for Global Optimization[END_REF]. Equations f k (x) = 0 are relaxed by inequalities -ǫ eq ≤ f k (x) ≤ ǫ eq , with ǫ eq = 1.e-8. The main results are the following. A rst experiment highlights the benets of this upper bounding, compared to a simple probing in every explored outer box. A second experiment underlines that, in a large majority of the tested systems, the upper bounding is satisfactory since the upper bound converges faster than the lower bound towards the nal value. Third, a qualitative study determines which of the two inner region extraction heuristics is the most useful in every system. A last study analyzes the size of the outer boxes in which the algorithms succeed in extracting an inner region.

  Thus, if a box is returned by the handling of the last constraint, this box is inner w.r.t. all the constraints. The handling of an individual constraint in InHC4 is radically dierent from[START_REF] Collavizza | Comparing Partial Consistencies[END_REF][START_REF] Benhamou | Universally Quantied Interval Constraints[END_REF]. Contrarily to their refutation process, our InHC4-Revise procedure tries to extract an inner region at each operator of the constraint. Like the main procedure of the state-of-the-art constraint propagation algorithm HC4[START_REF] Benhamou | Revising Hull and Box Consistency[END_REF][START_REF] Messine | Méthodes d'Optimisation Globale basées sur l'Analyse d'Intervalle pour la Résolution des Problèmes avec Contraintes[END_REF], our InHC4-Revise procedure (InHC4R) works with a tree representation of the constraint, as illustrated in Fig.1.

Figure 1 .

 1 Figure 1. Binary tree representation of the constraint 10yxy 2 ≤ 0. Left: First bottom-up evaluation phase. Right: top-down inner projection phase. Let us denote by [x] the input box and g j (x) ≤ 0 the constraint. Each node of the tree is associated to an interval, the intervals related to the leaves are initialized with the corresponding values in [x]. Then, the following two phases are performed:

Figure 2 .

 2 Figure 2. The dotted box corresponds to a maximal inner box of [z] w.r.t. the constraint g(x1, x2) ≤ 0. A point ẋ1 is randomly picked inside the range of allowed values illustrated by the horizontal segment. Only one remaining value ẋ2 can then make the computed inner box maximal.

Figure 3 .

 3 Figure 3. Inner projection for the binary multiplication. Left: Two maximal boxes that can indierently be computed by MonoMaxInnerBox in the two disjoint inner regions (quadrants) dened by the operator x1 * x2 ∈ [z] ≥ 0. Middle and right: Maximal box computed for x1 * x2 ∈ [z] ∋ 0 (z ≥ -z) with four calls to MonoMaxInnerBox (boxes in grey).

We consider minimization in this paper without loss of generality.

An interval [x i ] = [x i , x i ] denes the set of reals x i s.t. x i ≤ x i ≤ x i . A box [x] is the Cartesian product of intervals [x 1 ] × ... × [x i ] × ... × [xn].


Anne-Céline Boulanger

Clément Cancès

Hélène Mathis

Khaled Saleh

Nicolas Seguin

OSAMOAL: optimized simulations by adapted models using asymptotic limits

come

Introduction 1. Models and asymptotic limits 1.1. The Goldstein-Taylor model 1.2. The p-system 2. Asymptotic preserving schemes 2.1. The Goldstein-Taylor model 2.2. The p-system 3. Interface coupling between the fine model and the coarse model

Introduction

In simulation of complex flows for industrial purposes, it is natural to replace an accurate model by a simpler one in the regions of the computational domain where the flow is supposed to be almost at rest for the sake of reducing the computational cost. Let us give an example. Assume that the flow under consideration is a twophase flow, composed by air and water. In the highly heterogeneous regions, a complex model which accounts for phenomena with very small time and space scales has to be used such as the Baer-Nunziato model [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] where each phase is described by its own density, velocity and pressure and global relaxation terms represent the effects of the differences of velocity and pressure. In some smoother regimes, drift-flux models can be used [START_REF] Hibiki | One-dimensional drift-flux model and constitutive equations for relative motion between phases in various two-phase flow regimes[END_REF], where the pressures are equal and the relative velocity between the phases is a given function of the characteristic variable describing the flow (densities of each phases, mean velocity and mean pressure).

These regions where the simpler (or coarser) model can be used can be determined empirically or by the use of preliminary computations. In this work, we continue the program initiated in [START_REF] Mathis | Model adaptation for hyperbolic systems with relaxation[END_REF] and [START_REF] Cancès | Model adaptation for hyperbolic systems with relaxation[END_REF] about model adaptation. This techniques rely on automatically selecting the "good" model among a hierarchy of models, compatible by asymptotic arguments. The notion of "good" model means that we want to maximize the local use of coarse models without deteriorating the accuracy of the numerical results. Here, we focus on parabolic limits of hyperbolic models, while hyperbolic limits were studied in [START_REF] Mathis | Model adaptation for hyperbolic systems with relaxation[END_REF] and [START_REF] Cancès | Model adaptation for hyperbolic systems with relaxation[END_REF]. We restrict ourselves to standard systems and we do not investigate the dynamic part of the adaptation: we focus on the construction of numerical schemes which preserve the parabolic asymptotic limit and on the interfacial coupling between the hyperbolic system and its parabolic limit.

In the first section, we present the two models we are dealing with: the Goldstein-Taylor model and the p-system with friction. We also provide their asymptotic limit, which turns out to be a parabolic equation. In the second section, we derive and analyze numerical schemes which are able to reproduce the parabolic limits at the discrete level, namely asymptotic preserving schemes. Up to our knowledge, the results on the p-system are new. Note that we try to provide accurate stability estimates in order to fully understand the robustness of the CFL condition when ε → 0. The next part is devoted to the coupling between the hyperbolic model and its parabolic limit through a fixed thin interface. It has to be done carefully because classical methods of coupling fails to give accurate results. In order to illustrate our developments, preliminary numerical illustrations are also provided.

Models and asymptotic limits

The Goldstein-Taylor model

Let us consider the following fine model, called the Telegraph equation, the Goldstein-Taylor system or the p-system with friction for a Chaplygin gas:

   ε∂ t v + ∂ x u = 0, ε∂ t u + a 2 ∂ x v = -σ ε u, (M GT f )
where σ is a positive friction coefficient and a the sound speed.

An asymptotic analysis when ε → 0, using the Chapman-Enskog expansion for instance, leads to a coarse model, which is nothing but the heat equation:

   ∂ t v - a 2 σ ∂ xx v = 0, u = 0. (M GT c )

The p-system

Let us now consider the p-system with a friction term

   ε∂ t τ -∂ x u = 0, ε∂ t u + ∂ x P (τ ) = -σ ε u, (M p-sys f )
where τ is the specific volume, u the velocity and σ is a positive friction coefficient. The function P denotes the pressure law and satisfies classical assumptions: the function P is decreasing, convex, and tends to 0 as τ → ∞.

Denoting P a primitive form of P , the homogeneous p-system

ε∂ t τ -∂ x u = 0, ε∂ t u + ∂ x P (τ ) = 0 (1)
admits an entropy η defined by η(τ, u) = u 2 2 -P(τ ). Remark 1.1. For τ ⋆ fixed in (0, +∞), the strictly convex function

η τ⋆ (u, τ ) = u 2 2 -(P(τ ) -P(τ ⋆ )) + P (τ ⋆ )(τ -τ ⋆ ) ( 2
)
is also an entropy for the system (1) (one could also use the term of relative entropy) and satisfies η τ⋆ (u, τ ) ≥ 0, and η τ⋆ (0, τ ⋆ ) = 0.

The associated limit model, once again obtained by a Chapman-Enskog expansion, is

   ∂ t τ + 1 σ ∂ xx P (τ ) = 0, u = 0. (M p-sys c )
Since the pressure law P is a decreasing function, the first equation is a nonlinear parabolic equation.

Asymptotic preserving schemes

The aim of this part is to present the design of numerical schemes which are compatible with the asymptotic limits of the previous section, which commonly are called asymptotic preserving schemes [START_REF] Franck | Asymptotic preserving finite volumes discretization for non-linear moment model on unstructured meshes[END_REF][START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF][START_REF] Gosse | Space localization and well-balanced schemes for discrete kinetic models in diffusive regimes[END_REF][START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF][START_REF] Jin | Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations[END_REF].

Definition 2.1. A numerical scheme for system (M GT f ) (respectively (M p-sys f
)) is said to be asymptotic preserving if it is stable and consistent with the solutions of the hyperbolic model

(M GT f ) (resp. (M p-sys f
)) for all ε > 0 and if, at the limit ε → 0, it converges to a stable and consistent numerical scheme with the solutions of the limit parabolic model (M GT c ) (resp. (M p-sys c

)).

One of the important points is that the CFL condition for the approximation of the fine model becomes a CFL condition for the approximation of the coarse model. In our case, since we only are dealing with explicit schemes, this means that this condition for ε > 0 is of the form ∆t ≤ C∆x while at the limit ε → 0, it is of the form ∆t ≤ C∆x 2 .

In both cases, we construct first a well-balanced scheme, following [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF] and [START_REF] Gosse | Un schéma-équilibre adapté aux lois de conservation scalaires non-homogènes[END_REF]. We see that the CFL condition induced by a classical stability analysis is too restrictive to be valid for all ε > 0. Therefore, we propose an implicit version of this scheme in order to make it asymptotic preserving. Note that this new scheme is explicit in practice (the matrix to invert actually is diagonal). We then recover the pioneer work of Gosse and Toscani [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF][START_REF] Gosse | Space localization and well-balanced schemes for discrete kinetic models in diffusive regimes[END_REF] in the case of the Goldstein-Taylor and combining their method with a HLL (or equivalently relaxation) scheme, we easily extend it to the nonlinear case of the p-system. In both cases, we provide energy or entropy estimates together with CFL conditions and we carefully examine the behavior of these CFL conditions as ε → 0.

The Goldstein-Taylor model

We first present the well-balanced scheme for the Goldstein-Taylor model and show that it is not asymptotic preserving. We then propose an implicit modification to correct it, recovering the Gosse-Toscani scheme [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF][START_REF] Gosse | Space localization and well-balanced schemes for discrete kinetic models in diffusive regimes[END_REF].

The well-balanced scheme

The scheme of [START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter-dependent source. The case of Euler system with friction[END_REF] and the scheme of [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF] are both based on the construction of well-balanced scheme, following [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF] and [START_REF] Gosse | Un schéma-équilibre adapté aux lois de conservation scalaires non-homogènes[END_REF]. As a consequence, applied to system (M GT f ), these schemes are identical and the authors of these works claim that this scheme is asymptotic preserving. Their analysis is based on arguments of consistency and on numerical tests, because the models they investigate are nonlinear systems and a complete analysis of stability seems to be out of reach.

Actually, this scheme has been initially introduced in [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF] in the context of system (M GT f ). Since this system is linear, a complete analysis of the scheme can be carried out and the conclusion is that this scheme is not asymptotic preserving in the sense of definition 2.1.

In the sequel, ∆x, ∆t denote some positive parameters corresponding respectively to the space and time steps. Let u 0 , v 0 ∈ L 2 (R) be initial data, then we discretize them as

u 0 i = 1 ∆x (i+1/2)∆x (i-1/2)∆x u 0 (x)dx, v 0 i = 1 ∆x (i+1/2)∆x (i-1/2)∆x v 0 (x)dx, ∀i ∈ Z. (3)
Thanks to Jensen's inequality, u 0 i i and v 0 i i belong to ℓ 2 (Z) and

u 0 i i ℓ 2 (Z) ≤ u 0 L 2 (R) , v 0 i i ℓ 2 (Z) ≤ v 0 L 2 (R) .
For i ∈ Z and n ∈ N * , we denote by

W n i = (u n i , v n i )
T the discrete unknowns of the well-balanced finite volume scheme. The derivation of such a scheme, which is described in the works cited above, is based on the introduction of the function χ(t, x) = x. System (M GT f ) is replaced by the fully first-order system

       ε∂ t v + ∂ x u = 0, ε∂ t u + a 2 ∂ x v + σ ε u∂ x χ = 0, ∂ t χ = 0,
and the associated Riemann problem can be solved, with χ r -χ l = ∆x in order to mimic the jump of χ at each interface of the mesh. The eigenvalues of this system are 0 and ±a/ε and the two intermediate states are (v -, ū/K ε , χ l ) for -a/ε < x/t < 0 and (v + , ū/K ε , χ r ) for 0 < x/t < a/ε, where

K ε = 1 + σ∆x 2aε and ū(W l , W r ) = u l + u r 2 - a 2 (v r -v l ), v-(W l , W r ) = v l - 1 a (ū(W l , W r )/K ε -u l ), v+ (W l , W r ) = v r + 1 a (ū(W l , W r )/K ε -u r ).
As a result, the well-balanced finite volume scheme writes

v n+1 i = v n i - ∆t εK ε ∆x ū(W n i , W n i+1 ) -ū(W n i-1 , W n i ) , (4)
u n+1 i = u n i - a 2 ∆t ε∆x v-(W n i , W n i+1 ) -v+ (W n i-1 , W n i ) . (5)
Equation ( 5) can also be written as

u n+1 i = u n i - a 2 ∆t εK ε ∆x v(W n i , W n i+1 ) -v(W n i-1 , W n i ) - σ∆t ε 2 K ε u n i , (6)
where

v(W l , W r ) = v l + v r 2 - 1 2a (u r -u l ). (7)
Besides, the scheme (4) for the discretization of v writes

v n+1 i = v n i - ∆t εK ε ∆x u n i+1 -u n i-1 2 + a 2 (2v n i -v n i+1 -v n i-1 ) . (8)
Proposition 2.1. Under the CFL condition

1 ε + σ∆x 2a a∆t ∆x + σ∆t 2ε ε + σ∆x 2a ≤ 1, ( 9
)
the scheme (6)-( 8) is L 2 -diminishing, i.e. for all n ∈ N,

i∈Z ∆x (u n+1 i ) 2 + (v n+1 i ) 2 ≤ i∈Z ∆x (u n i ) 2 + (v n i ) 2 . ( 10
)
Proof. In order to prove stability of our scheme under a given CFL condition, we adopt the Von Neumann analysis. Assuming the solution is L 2 , we can introduce its Fourier modes in the following manner:

u n (x) = k∈Z ûn (k)e i2πkx .
Hence,

u n (x + ∆x) = u n (x)e i2πk∆x
, and the application of the Fourier transform on ( 6)-( 8) leads to:

  ûn+1 (k) vn+1 (k)   = A exp (k)   ûn (k) vn (k)   ,
where, denoting by

λ := a∆t εK ε ∆x , µ := σ∆t ε 2 K ε , α k := kπ∆x,
the amplification matrix A exp (k) can be written

A exp (k) :=   1 -2λ sin 2 (α k ) -µ -iλa sin(2α k ) -i λ a sin(2α k ) 1 -2λ sin 2 (α k )   .
Define

δ ′ := µ 2 4 -λ 2 sin 2 (2α k ),
then the eigenvalues ν ± ∈ C of A exp (k) are given by

ν ± := 1 -2λ sin 2 (α k ) - µ 2 ± √ δ ′ .
Assume first that δ ′ ≥ 0 (this is the case for ε large enough). In this case, one has

|δ ′ | ≤ µ 2 then ν + ≤ 1, while ν -≥ 1 -2λ -µ.
One obtains thus that under the CFL condition ( 9), one has ν -≥ -1, ensuring the L 2 stability of the scheme. Assume now that δ ′ < 0, so that

|ν + | 2 = |ν -| 2 = 1 -2λ sin 2 (α k ) - µ 2 2 + λ 2 sin 2 (2α k ) - µ 2 4 ≤ 1 + 4λ sin 2 (α k ) λ + µ 2 -1 .
Hence, the scheme is L 2 -diminishing under the CFL condition [START_REF] Franck | Asymptotic preserving finite volumes discretization for non-linear moment model on unstructured meshes[END_REF].

It particularly follows from ( 9) that ∆t has to fulfill the ε-dependent estimate

∆t ≤ 2ε ε + σ∆x 2a σ .
Hence, for a fixed space step ∆x, letting ε tend to 0 requires the time step ∆t to tend also to 0. This can also be remarked by taking constant initial data. Using the form ( 8)-( 6), one can directly see that this numerical scheme becomes a modified explicit Euler method for the equation ε∂ t u = -(σ/ε)u. It is well-known that such a method can only be stable under a restrictive condition on the time step ∆t, which makes it tend to 0 when ε → 0.

Implicit modification of the well-balanced scheme

As mentioned above, the time step ∆t tends to zero when ε tends to 0. The simple, but crucial, idea is then to take an implicit discretization of the source term, as done in [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF] (see also [START_REF] Gosse | Space localization and well-balanced schemes for discrete kinetic models in diffusive regimes[END_REF]). This leads to the new numerical scheme

u n+1 i = u n i - ∆t εK ε ∆x a 2 v n i+1 -v n i-1 2 + a 2 (2u n i -u n i+1 -u n i-1 ) - σ∆t ε 2 K ε u n+1 i , (11)
v n+1 i = v n i - ∆t εK ε ∆x u n i+1 -u n i-1 2 + a 2 (2v n i -v n i+1 -v n i-1 ) . ( 12
)
The purpose of the following proposition is to show that the new scheme ( 11)-( 12) is L 2 -diminishing for a choice of ∆t depending on ∆x but not on ε.

Proposition 2.2. Under the CFL condition

1 ε + σ∆x 2a a∆t ∆x ≤ 1, ( 13
)
the scheme [START_REF] Godlewski | The numerical interface coupling of nonlinear hyperbolic systems of conservation laws. II. The case of systems[END_REF], ( 12) is L 2 -diminishing in the sense of [START_REF] Gallouët | Compactness of discrete approximate solutions to parabolic PDEs-Application to a turbulence model[END_REF]. The inequality [START_REF] Gosse | Un schéma-équilibre adapté aux lois de conservation scalaires non-homogènes[END_REF] holds for all ε > 0 if ∆t satisfies the usual parabolic stability condition

∆t ≤ σ 2a 2 ∆x 2 , ( 14
)
which does not depend on ε.

Proof. The scheme ( 11)-( 12) can be obtained by splitting the convection and the source term, i.e.

• step 1: compute the convective contribution, i.e.

u n+1/2 i = u n i - ∆t εK ε ∆x a 2 v n i+1 -v n i-1 2 + a 2 (2u n i -u n i+1 -u n i-1 ) , (15)
v n+1/2 i = v n i - ∆t εK ε ∆x u n i+1 -u n i-1 2 + a 2 (2v n i -v n i+1 -v n i-1 ) ; (16)
• step 2: take the reaction term into account, i.e.

u n+1 i 1 + σ∆t ε 2 K ε = u n+1/2 i , v n+1 i = v n+1/2 i . (17)
Clearly, the step 1 is L 2 -diminishing under the CFL condition [START_REF] Gosse | Un schéma-équilibre adapté aux lois de conservation scalaires non-homogènes[END_REF]. Indeed, the scheme ( 15), ( 16) corresponds to the former scheme ( 6)-( 8) in the particular case where σ = 0. Therefore, its stability follows from Proposition 2.1. The step 2 is obviously L 2 -diminishing since for all i, u n+1 i

2 ≤ u n+1/2 i 2 .
We establish now that the scheme ( 11)-( 12) is asymptotic preserving, in the sense that, for fixed ∆x and for a convenient choice of ∆t, then, letting ε tend to 0, the solution

W ε,n+1 i = u ε,n+1 i , v ε,n+1 i T
of ( 11)-( 12) tends towards W 0,n+1 i which is the solution of a finite volume scheme which is consistent with the limit problem (M GT c ).

Proposition 2.3. Let W ε,n+1 i i∈Z,n∈N
be the solution of (11)-( 12) corresponding to data u 0 i i∈Z and v 0 i i∈Z belonging to ℓ 2 (Z). Assume that (14) holds, then, for all i ∈ Z, u ε,n+1 i tends to 0 as ε → 0 and v ε,n+1 i tends towards the solution v 0,n+1 i of the consistent discretization of (M GT c ) given by

v 0,n+1 i = v n i - a 2 ∆t σ∆x 2 2v n i -v n i-1 -v n i+1 , if n ≥ 1, ( 18
) v 0,1 i = v 0 i - a 2 ∆t σ∆x 2 2v 0 i -v 0 i-1 -v 0 i+1 - a∆t σ∆x 2 [u 0 i+1 -u 0 i-1 ]. ( 19
)
As a consequence, the numerical scheme (11)-( 12) under the CFL condition (13) is asymptotic preserving in the sense of definition 2.1.

Proof. Assume that [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF] holds, so that it follows from the proof of Proposition 2.2 that for all ε belonging to an unlabeled sequence of R * + tending to 0,

i∈Z ∆x u ε,n+1/2 i 2 + v ε,n+1 i 2 ≤ i∈Z ∆x u 0 i 2 + v 0 i 2 < +∞,
where u ε,n+1/2 i is given by [START_REF] Gosse | Space localization and well-balanced schemes for discrete kinetic models in diffusive regimes[END_REF]. Since ∆x is fixed, we know that there exists C ≥ 0 not depending on ε such that

∀i ∈ Z, u ε,n+1/2 i , v ε,n+1 i ≤ C. (20)
Therefore, since σ > 0 and since ε 2 K ε → 0 as ε → 0, it follows from (17) that u ε,n+1 i → 0 as ε → 0.

Thanks to [START_REF] Jin | Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations[END_REF], there exists an unlabeled subsequence of v ε,n+1 i ε converging towards some v 0,n+1

i with v 0,n+1 i i ∈ ℓ 2 (Z).
Using the fact that εK ε → σ∆x 2a , we obtain that ( 18) and ( 19) are fulfilled at the limit ε → 0.

Remark 2.1. As it appears in the first iteration of the limit scheme [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF], an additional quantity

r 0 i := a∆t σ∆x 2 u 0 i-1 -u 0 i+1
has to be taken in consideration. It is easy to check that for u 0 ∈ L 2 (R), thanks to [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF], the quantity r 0 i tends to 0 in ℓ 2 (Z) as ∆x, ∆t → 0. More precisely, if u 0 ∈ H 1 (R), there exists C depending neither on ∆x nor on ∆t such that r 0 i ℓ 2 (Z) ≤ C∆x. Denote by

A :    ℓ 2 (Z) → ℓ 2 (Z) w i → w i -a 2 ∆t σ∆x 2 [2w i -w i-1 -w i+1 ] ,
and assume that there exists ζ ∈ (0, 1) such that the following more restrictive CFL condition holds:

4a 2 ∆t σ∆x 2 ≤ ζ.

Then the mapping A is invertible and

A -1 ≤ 1 1-ζ .
Therefore, there exists a unique ṽ0 i ∈ ℓ 2 (Z), such that Aṽ 0 i = r 0 i . Moreover, ṽ0 i tends to 0 in ℓ 2 (Z) as ∆t, ∆x → 0 since r 0 i does. The solution (v n i ) i,n of the scheme (18)-( 19) is then the solution of the usual approximation of the heat equation [START_REF] Hibiki | One-dimensional drift-flux model and constitutive equations for relative motion between phases in various two-phase flow regimes[END_REF] for the slightly modified initial data v 0 i + ṽ0 i , the correction ṽ0 i i∈Z tending to zero with the discretization parameters ∆x, ∆t. Thanks to the continuity of the solution v to the heat equation with respect to the initial data v 0 , this allows to claim that the discrete solution provided by the scheme (18)-( 19) converges in L 2 (R) towards the unique solution v of the heat equation

∂ t v - a 2 σ ∂ xx v = 0, v |t=0 = v 0 .

The p-system

Let us now consider the case of the p-system.

The well-balanced scheme

Once again, the schemes proposed in [START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter-dependent source. The case of Euler system with friction[END_REF] and in [START_REF] Berthon | Asymptotic preserving HLL schemes[END_REF] follow the construction of well-balanced schemes and lead to the same numerical scheme. It is based on the HLL scheme [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] and can also be derived by the mean of a relaxation procedure (see e.g. [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF][START_REF] Chalons | Godunov-type schemes for hyperbolic systems with parameter-dependent source. The case of Euler system with friction[END_REF][START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF]). It writes, for all i ∈ Z and n ∈ N,

τ n+1 i = τ n i - ∆t ε∆x u n i-1/2 -u n i+1/2 , u n+1 i = u n i - ∆t ε∆x Π n,- i+1/2 -Π n,- i-1/2 , (21)
where u n i+1/2 , Π n,± i+1/2 and τ n,± i+1/2 are given by

u n i+1/2 = 1 1 + σ∆x 2εa u n i + u n i+1 2 + P (τ n i ) -P (τ n i+1 ) 2a , , (22)
Π n,- i+1/2 =P (τ n i ) -a(u n i+1/2 -u n i ), (23)
Π n,+ i+1/2 =P (τ n i+1 ) -a(u n i+1 -u n i+1/2 ). ( 24
)
In the above system, a > 0 is a parameter to be fixed later (cf. Proposition 2.4). By taking ( 22)-( 24) in consideration in [START_REF] Mathis | Model adaptation for hyperbolic systems with relaxation[END_REF] and introducing L ε = 1 + σ∆x 2εa , this scheme can be reformulated as

τ n+1 i = τ n i - ∆t εL ε ∆x u n i-1 -u n i+1 2 + 1 2a -2P (τ n i ) + P (τ n i+1 ) + P (τ n i-1 ) , (25)
u n+1 i = u n i - ∆t εL ε ∆x P (τ n i+1 ) -P (τ n i-1 ) 2 + a 2 2u n i -u n i-1 -u n i+1 - σ∆t ε 2 L ε u n i . (26)
We fail to obtain a complete analysis of stability for the scheme (25), (26). Nonetheless, a simple case can be used to check that this numerical scheme also suffer from a too restrictive stability condition of the same type as [START_REF] Franck | Asymptotic preserving finite volumes discretization for non-linear moment model on unstructured meshes[END_REF]. Assume that ∀i ∈ Z τ n i = τ , u n i = ū, where τ and ū are two positive constant. Then, the scheme (25)-(26) writes

∀i ∈ Z, τ n+1 i = τ , u n+1 i = 1 - σ∆t ε 2 L ε ū.
This formula corresponds to an explicit discretization of the ordinary differential equation

τ ′ (t) = 0, u ′ (t) = - σ ε u.
As a consequence, a necessary condition of the stability (in L 2 for instance) of the scheme is

σ∆t 2ε 2 L ε ≤ 1.
Here again, this condition implies that ∆t → 0 when ε → 0, which means that the scheme ( 21) is not asymptotic preserving.

Implicit modification of the well-balanced scheme

Following the idea used for the Goldstein-Taylor system, we replace the explicit discretization of the source term of ( 26) by an implicit one, which yields

τ n+1 i = τ n i - ∆t εL ε ∆x u n i-1 -u n i+1 2 + 1 2a -2P (τ n i ) + P (τ n i+1 ) + P (τ n i-1 ) , (27)
u n+1 i = u n i - ∆t εL ε ∆x P (τ n i+1 ) -P (τ n i-1 ) 2 + a 2 2u n i -u n i-1 -u n i+1 - σ∆t ε 2 L ε u n+1 i . ( 28
)
While we have not been able to completely study the stability of the numerical scheme (21) (but we proved that the associated CFL condition is too restrictive), we can prove that this new scheme is stable under an ε-independent CFL condition and therefore, it is asymptotic preserving. Note that the stability we invoke is the decrease of the entropy.

Let τ ⋆ ∈ (0, ∞), let η τ⋆ be the entropy of (1) defined by ( 2), and let τ 0 , u 0 be measurable functions such that

R η τ⋆ (τ 0 , u 0 )dx < +∞, (29)
and such that (τ 0 -τ ⋆ ) ∈ L ∞ ∩ L 1 (R). We moreover assume that τ 0 is essentially bounded away from 0. The initial data τ 0 , u 0 are discretized by

τ 0 i = 1 ∆x (i+1/2)∆x (i-1/2)∆x τ 0 (x)dx, u 0 i = 1 ∆x (i+1/2)∆x (i-1/2)∆x u 0 (x)dx,
so that, thanks to Jensen inequality,

i∈Z η τ⋆ (τ 0 i , u 0 i )∆x ≤ R η τ⋆ (τ 0 , u 0 )dx < +∞.
In particular, there exists τ , τ > 0 such that 0 < τ ≤ τ ⋆ ≤ τ < ∞, and such that τ ≤ τ 0 ≤ τ a.e. in R.

Proposition 2.4. Under the Whitham condition

a 2 ≥ 2 max n∈N max i∈Z -P ′ (τ n i ), -P ′ (τ n,± i+1/2 ) , (30)
where the quantities τ n,± i+1/2 are defined by

τ n,- i+1/2 = τ n i + 1 a u n i+1/2 -u n i , τ n,+ i+1/2 = τ n i+1 + 1 a u n i+1 -u n i+1/2 ,
and the CFL condition

2a∆t ∆x ε + σ∆x 2a ≤ 1, ( 31
)
then the scheme is entropy dissipative, i.e. the solution τ n+1 i , u n+1 i i∈Z of the scheme (27), (28)

satisfies i∈Z η τ⋆ (τ n+1 i , u n+1 i )∆x ≤ i∈Z η τ⋆ (τ n i , u n i )∆x. ( 32
)
Proof. As previously in the proof of Proposition 2.2, making the source term implicit allows to write the scheme ( 27), (28) as a splitting scheme.

• Step 1. We compute the convective contribution given by:

τ n+1/2 i = τ n i - ∆t εL ε ∆x u n i-1 -u n i+1 2 + 1 2a -2P (τ n i ) + P (τ n i+1 ) + P (τ n i-1 ) , u n+1/2 i = u n i - ∆t εL ε ∆x P (τ n i+1 ) -P (τ n i-1 ) 2 + a 2 2u n i -u n i-1 -u n i+1 .
This scheme is nothing but the usual HLL scheme to solve

∂ t τ -∂ x u = 0, ∂ t u + ∂ x P (τ ) = 0
with a time step ∆t = ∆t εLε . Following [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF], this scheme is entropy dissipative under the Whitham condition (30) and under the CFL condition (31).

• Step 2. The contribution of the source term is now computed with the implicit scheme

τ n+1 i = τ n+1/2 i , u n+1 i = 1 + σ∆t ε 2 L ε -1 u n+1/2 i . ( 33
)
In view of the definition (2) of the entropy η τ⋆ , this step is also entropy dissipative, so that the relation (32) holds.

Remark 2.2. As it was noticed in Proposition 2.2 for the case of the Goldstein-Taylor model, the CFL condition (31) is fulfilled for all ε > 0 if one requires that

∆t ≤ σ∆x 2 4a 2 . ( 34
)
As done previously on the Goldstein-Taylor model, we investigate the limit as ε tends to 0 for fixed discretization parameters of the discrete solution, denoted by (τ ε,n i , u ε,n i ) i,n , of the scheme ( 27)-(28).

Proposition 2.5. Under the Whitham condition (30), the CFL condition (34) and the finite entropy condition (29), then for all i ∈ Z, the limit of the scheme (27)-( 28) when ε → 0 writes

τ 0,n+1 i = τ 0,n i - ∆t σ∆x 2 -2P (τ 0,n i ) + P (τ 0,n i-1 ) + P (τ 0,n i+1 ) ∀n ≥ 1, (35)
τ 0,1 i = τ 0 i - ∆t σ∆x 2 -2P (τ 0 i ) + P (τ 0 i-1 ) + P (τ 0 i+1 ) + a(u 0 i-1 -u 0 i+1 ) . (36)
Proof. It follows from Proposition 2.4 and Remark 2.2 that under the Whitham condition (30), the CFL condition (34) and the finite entropy condition (29) on the initial data, one has, for all n ≥ 1,

i∈Z ∆x [η(τ ε,n i , u ε,n i )] ≤ R η(τ 0 , u 0 )dx < ∞. (37)
This ensures that there exists C ≥ 0 not depending on ε such that for all ε > 0, one has

∀i ∈ Z, ∀n ≥ 1, |u ε,n i | ≤ C, -P(τ ε,n i ) ≤ C.
Hence, there exists u 0,n i i,n and P 0,n i i,n such that up to an unlabeled subsequence,

∀i ∈ Z, u ε,n i → u 0 i , P(τ ε,n i ) → P 0,n i as ε → 0.
Since P is strictly increasing and continuous, then, denoting by τ 0,n

i = P -1 P 0,n i , we obtain that ∀i ∈ Z, τ ε,n i → τ 0,n i as ε → 0.
Thanks to (33), we obtain that for all n ≥ 1 and for all i ∈ Z, u 0,n i = 0. Since 2aεL ε → σ∆x as ε → 0 (recall that L ε = 1 + σ∆x 2εa ), considering the limit as ε → 0 in (27) yields (35), while a contribution coming from u 0 still appears in the computation (36) at the first time step. Obviously, since (35)-(36) define τ 0,n i i,n in a unique way, the whole unlabeled sequence (τ ε,n i ) i,n ε converges towards τ 0,n i i,n when ε → 0.

We provide now the proof of convergence of the scheme (35), (36) towards the unique weak solution of the problem

∂ t τ + 1 σ ∂ xx P (τ ) = 0, τ |t=0 = τ 0 , lim x→±∞ τ = τ ⋆ . ( 38
)
Lemma 2.6. Assume that a 2 > -P ′ (τ ), then there exists ω ∈ (0, τ ) and ∆x ω > 0 such that for all ∆x ≤ ∆x ω , one has

τ 0,n i ∈ [τ -ω, τ + ω], ∀i ∈ Z, ∀n ≥ 1. ( 39
)
Proof. Since a 2 > -P ′ (τ ), there exists ω > 0 such that a 2 = -P ′ (τ -ω). Denote by τ 0,1

i i the quantities defined by τ 0,1 i = τ 0 i - ∆t σ∆x 2 -2P (τ 0 i ) + P (τ 0 i-1 ) + P (τ 0 i+1 ) , ∀i ∈ Z.
Then it follows from classical monotonicity arguments that

τ ≤ τ 0,1 i ≤ τ , ∀i ∈ Z. Now, since τ 0,1 i = τ 0,1 i + aσ∆t ∆x 2 (u 0 i+1 -u 0 i-1
), we obtain, thanks to the CFL assumption (34) and from the regularity of u 0 that for ∆x small enough,

τ 0,1 i ∈ [τ -ω, τ + ω], ∀i ∈ Z.
Hence, the condition a 2 = -P ′ (τ -ω) ensures that the scheme is still monotone on this new range of values of τ 0,n i . Therefore, from a usual discrete maximum principle, we obtain by a straightforward induction that (39) holds.

We now state a L 1 -stability result, whose proof relies on classical monotonicity arguments (see e.g. [START_REF] Eymard | Finite volume methods[END_REF]).

Lemma 2.7. Let τ 0,n i i,n , τ 0,n i i,n
be solutions of (35)-( 36), corresponding to initial data τ 0 i i , τ 0 i i ∈ [τ , τ ], then, for a 2 > -P ′ (τ ), under the CFL condition (34), and for ∆x small enough, one has, for all n ≥ 0,

i∈Z ∆x τ 0,n+1 i -τ 0,n+1 i ≤ i∈Z ∆x τ 0,n i -τ 0,n i .
Denote by τ ∆x the piecewise constant function defined almost everywhere by

τ ∆x (x, t) = τ 0,n i if (x, t) ∈ ((i -1/2)∆x, (i + 1/2)∆x) × [n∆t, (n + 1)∆t).
From previous lemma, we deduce directly the following regularity result by choosing τ 0

i = τ 0 i-1 : Lemma 2.8. Assume that τ 0 ∈ L ∞ (R; [τ , τ ]) with τ 0 -τ ⋆ ∈ L 1 (R)
, assume that a 2 > -P ′ (τ ) and that (34) holds, and assume that ∆x is small enough, then, for all n ≥ 1,

i∈Z |τ 0,n i -τ 0,n i-1 |∆x ≤ i∈Z |τ 0 i -τ 0 i-1 |∆x + i aσ∆t ∆x 2 u 0 i+1 -u 0 i -(u 0 i-1 -u 0 i-2 ) .
In particular, under the CFL assumption (34) and if τ 0 ∈ BV (R) and u 0 ∈ W 1,∞ ∩ L 2 (R), then there exists C not depending on ∆x, ∆t such that

T V (τ ∆x (•, t)) ≤ T V (u 0 ) + C < ∞, ∀t ∈ R + . ( 40
)
Assume that the assumptions of Lemma 2.8 hold, then, thanks to the tools introduced in [10], we can claim that, up to an unlabeled subsequence of (τ ∆x ) ∆x , then there exists τ ∈ L ∞ ∩ BV (R) that satisfies

τ ∆x → τ in L 1 loc (R × R + ) as ∆x, ∆t → 0. (41)
It remains to check that τ satisfies (38), in particular that for all ψ ∈ D(R× R + ), the following weak formulation holds:

R×R+ τ ∂ t ψdxdt + R τ 0 ψ(•, 0)dx - R×R+ P (τ )∂ xx ψdxdt = 0. ( 42
)
The condition at x → ±∞ is fulfilled in the sense that, as a direct consequence of Lemma 2.7, one has

τ -τ ⋆ ∈ L 1 (R). Proposition 2.9. Assume that u 0 ∈ W 1,∞ (R), that τ 0 ∈ BV (R) such that τ 0 ∈ L ∞ (R; [τ , τ ]) and τ 0 -τ ⋆ ∈ L 1 (R)
, then, under the condition (34) with a 2 > -P ′ (τ ), the limit τ of τ ∆x exhibited in (41) is a solution of (38)

in the sense that τ -τ ⋆ ∈ L ∞ (R + ; L 1 (R))
and that the weak formulation (42) is satisfied. As a consequence, the numerical scheme (27)-( 28) is asymptotic preserving.

Proof. The convergence of τ ∆x towards τ in the good functional space has been already justified from compactness results. It remains to check that the weak formulation (42) is satisfied. Let ψ ∈ D(R × R + ), then we denote by

ψ n i = 1 ∆x (i+1/2)∆x (i-1/2)∆x ψ(x, n∆t)dx, , ∀i ∈ Z, ∀n ≥ 0.
Multiplying the scheme (35) by -ψ n+1 i ∆x, (36) by -ψ 1 i ∆x, then summing on i ∈ Z and n ∈ N yields, after reorganization of the sums,

A + B + C + D = 0, (43) where A
= i∈Z n≥0 τ 0,n i ψ n+1 i -ψ n i ∆t ∆x∆t, B = i∈Z τ 0 i ψ 0 i ∆x, C = i∈Z n≥0 P (τ n i ) 2ψ n i -ψ n i+1 -ψ n i-1 ∆x 2 ∆x∆t, D = i∈Z a∆t σ∆x 2 (u i+1 -u i-1 )ψ 1 i ∆x.
It follows from the L 1 convergence of τ ∆x towards τ and Lemma 2.6 that lim ∆t,∆x→0

A = R×R+ τ ∂ t ψdxdt, (44)
lim ∆x→0 B = R τ 0 ψ(•, 0)dx, ( 45
) lim ∆t,∆x→0 C = - R×R+ P (τ )∂ xx ψdxdt. (46)
Now, we deduce from the CFL condition (34) and from the regularity of u 0 that

|D| ≤ C∆x. (47)
Therefore, by considering the limit ∆x, ∆t → 0 in (43), it follows from (44)-(47) that τ satisfies (42).

Interface coupling between the fine model and the coarse model

In this section, we investigate the coupling through a thin interface between the hyperbolic model (M GT f ) or (M p-sys f ) with the parabolic model corresponding to its large time asymptotic limit, namely (M GT c ) or (M p-sys c ). As shown later in figure 2 for instance, a naive coupling by transmission of the information by Dirichlet boundary conditions in the spirit of [START_REF] Godlewski | The numerical interface coupling of nonlinear hyperbolic systems of conservation laws. II. The case of systems[END_REF][START_REF] Godlewski | The numerical interface coupling of nonlinear hyperbolic systems of conservation laws. I. The scalar case[END_REF] leads to inaccurate results. Therefore, we solve at the thin interface a coupled system in order to connect the fluxes of both models.

The Goldstein-Taylor model

Let us first illustrate the construction of the coupling condition on the simplest example, that is the Goldstein-Taylor model. Consider a fixed interface at x I+1/2 ∈ R. We aim at giving a numerical value for the fluxes at the interface x I+1/2 for the variables u, v, where, on the left side of the interface, the fine model is solved:

               ε∂ t v + ∂ x u = 0, ε∂ t u + a 2 ∂ x v = -σ ε u, u(x, 0) = u L , v(x, 0) = v L , for (x, t) in D f = (-∞, x I+1/2 ) × [0, ∞), (48)
while, on the other side of the interface, the coarse model is considered:

         ∂ t v -a 2 σ ∂ xx v = 0, u = 0, v(x, 0) = v R , for (x, t) in D c = (x I+1/2 , ∞) × [0, ∞). (49)
A first naive idea would be to solve the Riemann problem for the fine problem on both side of the interface, with the initial data u R = 0 and use the solution to construct numerical fluxes for both sides. As illustrated on Figure 2-right, this algorithm produces quite bad results. In this test case, the regime is such that the Goldstein-Taylor model (M GT f ) and its parabolic limit (M p-sys f ) provide similar results (actually, they are superposed on Figure 2-right). Nevertheless the naive coupling provides a very different solution which is discontinuous at x I+1/2 .

Therefore, another strategy has to be developed in order to improve the behavior of the coupled numerical solution and recover the continuity of v at x I+1/2 . To do so, we introduce some auxiliary unknowns U * = (u * , v * ) T on the left-hand side of the interface and solve a partial Riemann problem, see Figure 1. The Rankine-Hugoniot condition through the 1-wave of speed -a ε yields

u * -u L = a(v L -v * ). ( 50
)
By requiring the continuity of the flux of v through the interface x I+1/2 , we impose that This condition is very natural since it ensure the conservation of v. Combining (50) and (51), we obtain the intermediate values

u * ε = - a 2 σ v R -v * ∆x 2 . ( 51
) parabolic part t=0 t = ∆t U L = (u L , v L ) T v R U * hyperbolic part
v * (W L , v R ) = 1 1 + 2εa σ∆x v L + 2εa σ∆x v R + u L a . ( 52
) u * (W L , v R ) = 1 K ε (u L + a(v L -v R )) , (53)
It remains to plug these formulas in the numerical scheme ( 4) and ( 5) setting i = I and replacing the numer-

ical flux ū(W n I , W n I+1 )/K ε and v-(W n I , W n I+1 ) by u * (W n I , v n I+1 ) and v * (W n I , v n I+1
). Then, the source term is once again treated in an implicit way in order to obtain an asymptotic preserving scheme which, after some computations, leads to

v n+1 I = v n I - ∆t εK ε ∆x u n I + a(v n I+1 -v n I ) -ū(W n I-1 , W n I ) , (54)
u n+1 I = u n I - a 2 ∆t εK ε ∆x v n I+1 -v(W n I-1 , W n I ) - σ∆t ε 2 K ε u n+1 I , (55)
while the parabolic flux is the cell I + 1 is given by (34) and (53). Note that the numerical flux at the interface x I+1/2 in (54) satisfies

lim ε→0 u * ε ((v n I , 0), v n I+1 ) = 2a 2 σ v L -v R ∆x ,
which is not a consistent approximation for the flux of the expected heat equation that v should satisfy. Nevertheless, the resulting coupling method provides very good results. We illustrate this new coupling method in a simple test case. The 1D space domain is split in two parts: (-1, 0) and (0, 1). The model in the left part is the fine one (M GT f ) and in the right part, the model is the coarse one (M GT c ). We compare the numerical coupling method described above with the basic one, which amounts to set u n I+1 in the hyperbolic scheme ( 4)- [START_REF] Cancès | Model adaptation for hyperbolic systems with relaxation[END_REF] to define the numerical coupling flux. We also plot the results obtain with the full hyperbolic solution and the full parabolic solution. The initial condition is of Riemann type, with v L = 1 and v R = u L = u R = 0. The values of the parameters are set to: σ = 2 and a = 1. Each part is discretized with 30 cells. One can easily check that the basic coupling method introduces a jump at the coupling interface. Moreover, for ε = 0.1 (right figure), though the full hyperbolic solution and the parabolic solution are superposed, the basic coupling still provides a jump at the coupling interface. Actually, this test case was the motivation of the construction of this new coupling method. We can check that in this latter case the results of the new coupling method are very good. If we go back to the case ε = 0.2 (left figure), the hyperbolic solution

The p-system

We perform a similar coupling technique for the p-system. Here again the fine model (M p-sys f ) is solved on the left side of the interface and the coarse model (M p-sys c ) is solved on the right side of x I . Because of the nonlinearity of both models, the design of an exact coupling strategy with a partial Riemann problem would lead to a nonlinear method. In order to avoid this, we make use of a relaxation method [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF][START_REF] Coquel | Some new Godunov and relaxation methods for two-phase flow problems[END_REF] to linearize the p-system and we study the associated coupled partial Riemann problem in order to define the numerical fluxes. It takes the form

         ε∂ t τ -∂ x u = 0, ε∂ t u + ∂ x π(T, τ ) = -σ ε u, ∂ t T = λ(τ -T ), (τ, u, T )(x, 0) = (τ L , u L , τ L ), for (x, t) in D f = (-∞, x I ) × [0, ∞), (56)
where π(τ, T ) = P (T ) + a 2 (T -τ ) denotes the relaxed pressure and we assume that the Whitham condition

a 2 > max s |P ′ (s)| is satisfied. The coarse model (M p-sys c
) to be solved on the right side of the interface reads

     ∂ t τ + 1 σ ∂ xx P (τ ) = 0, u = 0, τ (x, 0) = τ R , for (x, t) in D c = (x I , ∞) × [0, ∞). (57)
Following the same strategy as the one used for the Goldstein-Taylor model, we introduce an interfacial state U * = (τ * , u * , π * ) T on the left hand side of the interface. The left state and the interfacial state U * are related by the Rankine-Hugoniot condition through the 1-wave of speed -a ε , i.e.

a(τ * -τ L ) = (u * -u L ) (58) -a(u * -u L ) = (π * -P (τ L )) (59)
since π L = P (τ L ). First, we assume that the pressure is continuous through the interface, that is to say that there is only one pressure at the interface, namely π * . In order to obtain a conservative method for τ , the associated flux has to be the same through the interface x I , which leads to the equality

- 1 ε u * = 1 σ P (τ R ) -π * ∆x 2 , (60)
which is nothing but the counterpart to (51). Combining equations ( 58), ( 59) and (60), we obtain

τ * = τ L + 1 a (u * -u L ), (61)
u * = 1 1 + σ∆x 2εa u L - 1 a (P (τ R ) -P (τ L )) , (62)
π * = 1 1 + 2εa σ∆x P (τ L ) + 2aε σ∆x P (τ R ) + au L . (63)
As in the case of the Goldstein-Taylor model, the values u * and π * are plugged in the well-balanced scheme [START_REF] Mathis | Model adaptation for hyperbolic systems with relaxation[END_REF] for the cell i = I. Following the same computations which leads to the scheme (25)-( 26) and after the implicit modification of the source term, one obtains

τ n+1 I = τ n I - ∆t εL ε ∆x
u n I + a(P (τ n I+1 ) -P (τ n I )) -u n I-1/2 , (64)

u n+1 I = u n I - ∆t εL ε ∆x P (τ n I+1 ) -Π(W n I-1 , W n I ) - σ∆t ε 2 L ε u n+1 I , (65)
where u n I-1/2 is given by ( 22) and Π(W L , W R ) = 1 2 (P (τ L ) + P (τ R )) -a 2 (u R -u L ).

Remark 3.1. Other coupling conditions could have been chosen, for instance the continuity of the specific volume τ through the interface instead of the continuity of the pressure. However this condition requires to solve a nonlinear problem in order to determine the interfacial state.

Let us now present the numerical results of two test cases. The pressure law is P (τ ) = τ -2 and σ = 2. The initial data corresponds to a Riemann problem with τ L = 2, τ R = 1 and u L = u R = 0. Once again, each part of the computational domain is discretized in 30 cells. Two tests with different values of ε are provided in Figure 3. In the right figure, ε is sufficiently small to obtain very close solutions for the fully hyperbolic case and for the fully parabolic case. The solution of our coupling method is also very close to these two solutions, as expected. In the contrary, the basic coupling method, which amounts to use the hyperbolic flux with u I+1 = 0, introduces a jump at the coupling interface. In the left figure, the solutions of the fully hyperbolic case and of the fully parabolic case are different since ε is greater than in the previous case. A jump of τ appears with the basic coupling method which is larger than in the case of a smaller ε. Indeed, some complementary numerical experiments let us think that the amplitude of this jump to tend to 0 when ε goes to 0. The result provided by our coupling method is very satisfying. It follows the shape of the hyperbolic solution in the left part and the shape of the parabolic solution in the right part. As in the Goldstein-Taylor case, the transmission through the coupling interface with our new coupling method is however not fully perfect. Indeed, a small jump can be detected for both values of ε which can be probably related to the lack of consistency of the coupling flux when ε tends to 0.

Conclusion

During this project, we have proposed asymptotic preserving schemes for the parabolic limit of the Goldstein-Taylor system and of the p-system. The core of these schemes and their analysis was already present in [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF][START_REF] Gosse | Space localization and well-balanced schemes for discrete kinetic models in diffusive regimes[END_REF] as far as the Goldstein-Taylor model is concerned. We also have provided an original coupling method between a hyperbolic system (with source term) and its parabolic limit in order to remove the jumps that the classical coupling method makes appear. The remaining work about the model adaptation concerns the dynamical part: definition of error indicators and development of the fully dynamic adaptive algorithm, following [START_REF] Cancès | Model adaptation for hyperbolic systems with relaxation[END_REF][START_REF] Mathis | Model adaptation for hyperbolic systems with relaxation[END_REF]. This is under investigation, as well as the extension to more complex models.

Figure 1 .

 1 Figure 1. The interface states U * = (u * , v * ) are linked to (u L , v L ) by the 1-wave of speed -a/ε, and to v R by imposing the continuity of the numerical flux of v.

Figure 2 .

 2 Figure 2. Approximate v with ε = 0.2 (left) and with ε = 0.1 (right) at t = 0.05.

Figure 3 .

 3 Figure 3. Approximate τ with ε = 0.4 (left) and with ε = 0.1 (right) at t = 0.2.

This project has been supported by the LRC Manon (Modélisation et Approximation Numérique Orientées pour l'énergie Nucléaire -CEA/DM2S-LJLL).