R Abgrall

A general framework to construct schemes satisfying additional conservation relations. Application to entropy conservative and entropy dissipative schemes

 in the sense that no particular constraints are set on quadrature formula and that a priori maximum accuracy can still be achieved. We study the behaviour of the method on a non linear scalar problem. However, the method is not restricted to scalar problems.

Introduction

In this paper, we are interested in the approximation of non-linear hyperbolic problems. To make things more precise, our target are the Euler equations in the compressible regime, other examples are the MHD equations. The case of parabolic problems in which the elliptic terms play an important role only in some area of the computational domain, such as the Navier-Stokes equations in the compressible regime, or the resistive MHD equations, can be dealt with in a similar way.

Let D an open subset of R p and d C 1 functions f j , j = 1, . . . , d defined on D. We consider the hyperbolic problem, with u :

Ω ⊂ R d → D ⊂ R p , defined by div f (u) = d j=1 ∂f j ∂x j (u) = 0 (1a)
subjected to

∇ u f (u) • n(x), 0 - (u -u b ) = 0 on ∂Ω. (1b)
In (1b), n(x) is the outward unit vector at x ∈ ∂Ω and u b is a regular enough function. We have also used standard notations: if A is a p × p diagonalisable matrix in R, A = LΛR where L and R are two p × p matrices with LR = Id p×p and Λ = diag(λ i ) i=1,...,p diagonal, then

A ± = LΛ ± R where Λ ± = diag(λ ± i ) i=1,...,p
and, for x ∈ R, x + = max(x, 0) and x -= min(x, 0). The weak formulation of (1) is: u ∈ L 2 (Ω) d ∩ L ∞ (Ω) d is a weak solution of (1) if u ∈ D and for any ϕ ∈ C 1 0 (Ω) d ,

-

Ω ∇ϕ • f (u) dx + ∂Ω ϕ F n (u, u b ) -f (u) • n dγ = 0 (2)
where F n is a flux that is almost everywhere the upwind flux:

F n (u, u g ) = ∇ u f (u) • n(x) + u + ∇ u f (u) • n(x) -u g
In general, schemes are constructed starting from the conservation relation. In some cases, such as when the solution is smooth enough, it appears that the solution of (1) also satisfies an additional conservation relation. There is no reason that the initial discretisation will also be consistent with the new conservation relation. A typical example is the entropy. If U , defined on D, is a convex and C 1 function and if there exists g, a C 1 function, g = (g 1 , . . . , g d ) : D → R d , such that: for all j = 1, . . . , d,

∇ u U • ∇ u f j = ∇ u g j , (3a)
then when u is smooth we get the additional conservation relation

div g(u) = d j=1 ∂g j ∂x j (u) = 0. ( 3b
)
An entropy solution of (1) will satisfy div g(u) ≤ 0 in the sense of distributions. Hence a rather natural question is: given a scheme for [START_REF] Chen | Entropy stable high order discontinuous Galerkin methods with suitable quadrature rules for hyperbolic conservation laws[END_REF], how can we modify it so that the modified scheme is consistent with an additional conservation relation, i.e. in the case of an entropy, how to modify it so that it is also compatible with an entropy inequality? This is not a new question: there has been lots of work on that particular question, see [START_REF] Tadmor | Handbook of Numerical Methods for Hyperbolic Problems, Basic and fundemental issues[END_REF][START_REF] Tadmor | The numerical viscosity of entropy stable schemes for systems of conservation laws[END_REF][START_REF] Tadmor | Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems[END_REF][START_REF] Jiang | On a cell entropy inequality for discontinuous Galerkin methods[END_REF][START_REF] Bouchut | A MUSCL method satisfying all the numerical entropies[END_REF][START_REF] Lefloch | Fully discrete, entropy conservative schemes of arbitrary order[END_REF], and more recently see [START_REF] Hou | Solutions of multi-dimensional hyperbolic systems of conservation laws by square entropy condition satisfying discontinuous Galerkin method[END_REF][START_REF] Ismail | Affordable, entropy-consistent Euler flux functions II: entropy production at shocks[END_REF][START_REF] Hicke | Multidimensional summation-by-part operators: general theory and application to simplex elements[END_REF][START_REF] Gassner | A skew-symmetric discontinuous Galerkin spectral element discretisation and its relation to SBP-SAT finite differennce methods[END_REF][START_REF] Hiltebrand | Entropy stable shock capturing spacetime discontinuous Galerkin schemes for systems of conservation laws[END_REF][START_REF] Guermond | Entropy viscosity method for nonlinear conservation laws[END_REF][START_REF] Chen | Entropy stable high order discontinuous Galerkin methods with suitable quadrature rules for hyperbolic conservation laws[END_REF][START_REF] Ranocha | Summation-by-parts operators for correction procedure via reconstruction[END_REF][START_REF] Parsani | Entropy stable staggered grid discontinuous spectral collocation methods of any order for the compressible Navier-Stokes equations[END_REF][START_REF] Svärd | Review of summation-by-parts schemes for initial-boundary-value problems[END_REF][START_REF] Hiltebrand | Entropy stable shock capturing spacetime discontinuous Galerkin schemes for systems of conservation laws[END_REF][START_REF] Svärd | Entropy-stable schemes for the Euler equations with far-field and wall boundary conditions[END_REF][START_REF] Ray | Entropy stable schemes for compressible Euler equations[END_REF][START_REF] Ray | Entropy stable scheme on two-dimensional unstructured grids for Euler equations[END_REF][START_REF] Gassner | A well balanced and entropy conservative discontinuous Galerkin spectral element method for the shallow water equations[END_REF][START_REF] Chandrasekar | Kinetic energy preserving and entropy stable finite volume scheme for compressible Euler and Navier-Stokes equations[END_REF] and references therein for an un-complete list of contributions. One of the salient aspects of research on this problem is a focus on the discrete version of the schemes: for example, when some integral appear in the formulation, as in finite element method, quadrature formula are needed, and the integrals are never computed exactly. The central question is then how to translate, from a discrete point of view, the 'continuous' differential relations (3a) and the algebra associated to them at the discrete level. Up to our knowledge, entropy stable schemes can be rigorously designed only for special quadrature formula. The objective of this paper is to provide a different point of view within the Residual Distribution framework. In particular, we show that we can consider general quadrature formula.

The outline of the paper is as follows. In a first part, we recall the class of schemes (nicknamed as Residual distribution schemes or RD or RDS for short) we are interested in, and show their link with more classical methods such as finite element ones. Then we recall a condition that guaranties that the numerical solution will converge to a weak solution of the problem, and a second one about entropy condition, and recall a systematic way to check the formal accuracy of the scheme. In the third part, we give a general recipe that show how to modify a scheme in order to fullfil an additional conservation law. The fourth part is devoted to a discussion on the entropy conservation: we show explicitly how to modify a given scheme. The next section is devoted to the construction of entropy dissipative schemes, and we provide numerical examples. A conclusion follows.

Notations

From now on, we assume that Ω has a polyhedric boundary in R d , d = 1, 2, 3. This simplification is by no mean essential. We denote by E h the set of internal edges/faces of a triangulation T h , and by F h those contained in ∂Ω. K stands either for an element K or a face/edge e ∈ E h ∪ F h . The boundary faces/edges are denoted by Γ. The mesh is assumed to be shape regular, h K represents the diameter of the element K.

Similarly, if e ∈ E h ∪ F h , h e represents its diameter. Throughout the text, we will restrict ourselves to the case where the elements K are simplicies, mostly for simplicity reasons.

Throughout this paper, we follow Ciarlet's definition [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF][START_REF] Ern | Theory and practice of finite elements[END_REF] of a finite element approximation: for any element K, we have a set of linear forms Σ K acting on the set P k of polynomials of degree k such that the linear mapping

q ∈ P k → σ 1 (q), . . . , σ |Σ K | (q)
is one-to-one. The space P k is spanned by the basis function {ϕ σ } σ∈Σ K defined by ∀σ, σ , σ(ϕ σ ) = δ σ σ .

We have in mind either the Lagrange interpolation where the degrees of freedom are associated to the Lagrange points in K, or other type of polynomials approximation such as using Bézier polynomials where we will also do the same geometrical identification. More specifically, if {A 1 , . . . , A d+1 } are the vertices of K (i.e K is the convex hull of the vertices {A 1 , . . . , A d+1 }), the Lagrange points of degree k are defined by their barycentric coordinates

L α = d+1 l=1 i l d + 1 A l , α = (i 1 , . . . , i d+1 ) (4)
with i l positive integers such that i 1 + . . . + i d+1 = k Any point M of R d can be uniquely represented by its barycentric coordinates {λ 1 (M), . . . , λ d+1 (M)} with [START_REF] Hicke | Multidimensional summation-by-part operators: general theory and application to simplex elements[END_REF], are the polynomials

M = d+1 l=1 λ l (M)A l , d+1 l=1 λ l (M) = 1. Note that λ l (M) ≥ 0 if and only if M ∈ K. The Bézier polynomials of degree k, associated to the multi-index α = (i 1 , • • • , i d+1 ), i l positive integer, i 1 + . . . + i d+1 = k as in
B α = k! i 1 ! • • • i d+1 ! λ i1 1 . . . λ i d+1 d+1 .
They sum up to unity, are positive on K, are interpolatory on the vertices only, span P k as well as the Lagrange polynomials of degree k. For the Lagrange interpolation of degree k, the set of degrees of freedom is the set of Lagrange points of degree k, i.e. the L α where α is as in [START_REF] Hicke | Multidimensional summation-by-part operators: general theory and application to simplex elements[END_REF]. For the Bézier approximation, we make the geometrical identification between the muti-index α that defines the polynomial B α and the Lagrange point L α , though the Bézier polynomials are not interpolatory in general. However, the Bézier polynomials of degree k span P k , as well as the Lagrange polynomials of degree k. Considering all the elements covering Ω, the set of degrees of freedom is denoted by S and a generic degree of freedom by σ. We note that for any K,

∀x ∈ K, σ∈K ϕ σ (x) = 1.
For any element K, #K is the number of degrees of freedom in K. If Γ is a face or a boundary element, #Γ is also the number of degrees of freedom in Γ.

The solution is sought for in the space V k h (or V h for short since the degree k will be assumed to be constant) defined by, setting

V h = K {v ∈ L 2 (K), v |K ∈ P k }.
We will consider two cases V h = V h ∩ C 0 (Ω) in which case the triangulation needs to be conformal, and V h = V h where the continuity requirement is released: the triangulation does not need to be conformal anymore.

Throughout the text, we need to integrate functions. This is done via quadrature formula, and the symbol used in volume integrals

K v(x) dx := |K| q ω K q v(x K q ) (5)
or boundary/face integrals

∂K v(x) dγ := |f | q ω f q v(x f q ) (6)
means that these integrals are done via user defined numerical quadratures. If e ∈ E h , represents any internal edge, i.e. e ⊂ K ∩ K + for two elements K and K + , we define for any function

ψ the jump [∇ψ] = ∇ψ |K -∇ψ |K + . Similarly, {v} = 1 2 v |K + v |K + . If
x and y are two vectors of R q , for q integer, x, y is their scalar product. In some occasions, it can also be denoted as x • y or x T y. Similarly, if x = (x 1 , . . . , x d ) is a vector and = ( 1 , . . . , d ) with j ∈ R m , then we denote

x • f = x T f = d i=1 x i f i .
Last, v h ∈ V h will denote any test function, u will be the conserved quantities defined in (1), V will be the entropy variable and V h ∈ V h will be its interpolant in P k . Depending on the context, u h ∈ V h will denote either an approximation of u or u(V h ). Similarly, π h (u) ∈ V h is the interpolant of u or u(π h (V )), depending on the context.

3 Schemes, conservation, entropy dissipation

Schemes

In order to integrate the steady version of (1) on a domain Ω ⊂ R d with the boundary conditions (1b), on each element K and any degree of freedom σ ∈ S belonging to K, we define residuals Φ K σ (u h ). Following [START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF][START_REF] Abgrall | High-order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF], they are assumed to satisfy the following conservation relations: For any element K, and any

u h ∈ V h , σ∈K Φ K σ (u h ) = ∂K fn (u h |K , u h |K -) dγ. ( 7a
)
where u |K is the restriction of u h in the element K while u h |K -is the restriction of u h on the other side of the local edge/face of K. In addition, fn is a consistant numerical flux, i.e. f (u, u) = f (u) • n. Similarly, we consider residuals on the boundary elements Γ. On any such Γ, for any degree of freedom σ ∈ S ∩ Γ, we consider boundary residuals Φ Γ σ (u h ) that will satisfy the conservation relation

σ∈Γ Φ Γ σ (u h ) = Γ fn (u h |Γ , u b ) -f (u h ) • n dγ (7b)
where, for x ∈ Γ, u h |Γ (x) = lim y→x,y∈Ω u h (y). Once this is done, the discretisation of (1) is achieved via: for any σ ∈ S, K⊂Ω,σ∈K

Φ K σ (u h ) + Γ⊂∂Ω,σ∈Γ Φ Γ σ (u h ) = 0. (8)
In [START_REF] Jiang | On a cell entropy inequality for discontinuous Galerkin methods[END_REF], the first term represents the contribution of the internal elements. The second exists if σ ∈ Γ and it represents the contribution of the boundary conditions. In fact, the formulation ( 8) is very natural. Consider a variational formulation of the steady version of (1):

find u h ∈ V h such that for any v h ∈ V h , a(u h , v h ) = 0.
Let us show on three examples that this variational formulation leads to [START_REF] Jiang | On a cell entropy inequality for discontinuous Galerkin methods[END_REF]. They are

• The SUPG [27] variational formulation, with u h , v h ∈ V h = V h ∩ C 0 (Ω): a(u h , v h ) := - Ω ∇v h • f (u h ) dx + K⊂Ω h K K ∇f (u h ) • ∇v h τ K ∇f (u h ) • ∇u h dx + ∂Ω v h F n (u h , u b ) -f (u h ) • n dγ. (9)
and τ K > 0.

• The Galerkin scheme with jump stabilization, see [START_REF] Burman | Edge stabilization for Galerkin approximation of convection-diffusionreaction problems[END_REF] for details. We have

a(u h , v h ) := - Ω ∇v h • f (u h ) dx + e⊂Ω θ e h 2 K e ∇v h • ∇u h dγ + ∂Ω v h F n (u h , u b ) -f (u h ) • n dγ. (10)
Here,

u h , v h ∈ V h = V h ∩ C 0 (Ω)
, and θ e is a positive parameter.

• The discontinuous Galerkin formulation: we look for

u h , v h ∈ V h = V h such that a(u h , v h ) := K⊂Ω - K ∇v h • f (u h )dx + ∂K v h • fn (u h |K , u h |K -) dγ . (11)
In [START_REF] Hou | Solutions of multi-dimensional hyperbolic systems of conservation laws by square entropy condition satisfying discontinuous Galerkin method[END_REF], the boundary integral is a sum of integrals on the faces of K, and here for any face of K u h

|K -

represents the approximation of u on the other side of that face in the case of internal elements, and u b when that face is on ∂Ω. In [START_REF] Hou | Solutions of multi-dimensional hyperbolic systems of conservation laws by square entropy condition satisfying discontinuous Galerkin method[END_REF], f is a consistent numerical flux. Note that to fully comply with (8), we should have defined for boundary faces u h |K -= u h |K , and then ( 11) is rewritten as

a(u h , v h ) := K⊂Ω - K ∇v h • f (u h )dx + ∂K v h • fn (u h |K , u h |K ) dγ + Γ⊂∂Ω Γ v h • fn (u h Γ , u b ) -f (u h ) • n dγ. ( 12
)
with by abuse of language, f = F on the boundary edges.

In the SUPG, Galerkin scheme with jump stabilisation or the DG scheme, the boundary flux can be chosen different from F. This can lead to boundary layers if these flux are not "enough" upwind, but we are not interested in these issues here. Using the fact that the basis functions that span V h have a compact support, then each scheme can be rewritten in the form [START_REF] Jiang | On a cell entropy inequality for discontinuous Galerkin methods[END_REF] with the following expression for the residuals:

• For the SUPG scheme [START_REF] Bouchut | A MUSCL method satisfying all the numerical entropies[END_REF], the residual are defined by

Φ K σ (u h ) = ∂K ϕ σ f (u h ) • n dγ - K ∇ϕ σ • f (u h ) dx + h K K ∇ u f (u h ) • ∇ϕ σ τ K ∇ u f (u h ) • ∇u h dx (13) with τ K > 0.
• For the Galerkin scheme with jump stabilization [START_REF] Lefloch | Fully discrete, entropy conservative schemes of arbitrary order[END_REF], the residuals are defined by:

Φ K σ (u h ) = ∂K ϕ σ f (u h ) • n dγ - K ∇ϕ σ • f (u h ) dx + e faces of K θ e h 2 e ∂K [∇u] • [∇ϕ σ ] dγ (14)
with θ e > 0. Here, since the mesh is conformal, any internal edge e (or face in 3D) is the intersection of the element K and an other element denoted by K + .

• For the discontinuous Galerkin scheme,

Φ K σ (u h ) = - K ∇ϕ σ • f (u h )dx + ∂K ϕ σ • fn (u h |K , u h |K -) dγ (15)
using the second definition of u h |K -. • The boundary residuals are

Φ Γ σ (u h ) = Γ ϕ σ fn (u h Γ , u b ) -f (u h ) • n dγ (16)
All these residuals satisfy the relevant conservation relations, namely (7a) or (7b), depending if we are dealing with element residuals or boundary residuals.

For now, we are just rephrasing classical finite element schemes into a purely numerical framework. However, considering the pure numerical point of view and forgetting the variational framework, we can go further and define schemes that have no clear variational formulation. These are the limited Residual Distribution Schemes, see [START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF][START_REF] Abgrall | High-order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids[END_REF], namely

Φ K σ (u h ) = β σ ∂K f (u h ) • n dγ (17)
or

Φ K σ (u h ) = β σ ∂K f (u h ) • n dγ + h K θ S K ∇ u f (u h ) • ∇ϕ σ τ K ∇ u f (u h ) • ∇u h dx (18) or Φ K σ (u h ) = β σ ∂K f (u h ) • n dγ + θ J h 2 K ∂K [∇u h ] • [∇ϕ σ ] dγ. (19)
or, combining the two kind of stabilisations

Φ K σ (u h ) = β σ ∂K f (u h )•n dγ+θ J h 2 K ∂K [∇u h ]•[∇ϕ σ ] dγ+h K θ S K ∇ u f (u h )•∇ϕ σ τ ∇ u f (u h )•∇u h dx (20)
where the parameters β σ are defined to guarantee conservation, σ∈K β σ = 1 and such that (18) without the streamline term and ( 19) without the jump terms satisfy a discrete maximum principle. The parameters τ K , θ J and θ S in ( 18), ( 19) and ( 20) are positive real numbers. The streamline term and jump term are introduced because one can easily see that spurious modes may exist, but their role is very different compared to ( 13) and ( 14) where they are introduced to stabilize the Galerkin scheme: if formally the maximum principle is violated, experimentally the violation is extremely small, if not existant. See [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF][START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF] for more details.

A similar construction can be done starting from a discontinuous Galerkin scheme without non-linear stabilisation such as limiting, has been applied.

The non-linear stability is provided by the coefficient β σ which is a non-linear function of u h . Possible values of β σ are described in the appendix A.

Conservation

It is time to state the geometrical assumptions we make on the quadrature formula: Assumption 3.1 (Assumption on the quadrature formula). In (7b) and (7a), numerical integrations are made on the faces of the elements of T h , boundary of Ω included. Hence when we write, for an element K ∂K ψ • n dγ, or for a face on ∂Ω Γ ψ dγ, we are adding contributions 'living' on the squeleton of T h . We request that the quadrature formula depend only on the squeleton elements, so that if

f = K 1 ∩ K 2 then f ∩∂K1 ψ • n dγ + f ∩∂K2 ψ • n dγ = 0. ( 21
)
In the previous relation, the first integral is a numerical formula on f seen from K 1 , and the second term is an integral seen from K 2 . Since the normal are opposite, this means that the geometrical location of the quadrature points are the same on f seen from K 1 and f seen from K 2 . If f is on the boundary of K, this is also the intersection of an element K and the boundary of Ω, and again, the quadrature points have the same geometrical locations whatever be the way we consider f .

In practice, this is not a restriction at all since quadrature points are defined by their barycentric coordinates and symmetric with respect to permutations of the vertices that respect orientation.

From [START_REF] Jiang | On a cell entropy inequality for discontinuous Galerkin methods[END_REF], using the conservation relations (7b) and (7a), we obtain for any

v h ∈ V h , v h = σ∈S v σ ϕ σ , the following relation: 0 = - Ω ∇v h • f (u h ) dx + ∂Ω v h fn (u h , u b ) -f (u h ) • n dγ + e∈E h e [v h ] fn (u h |K , u h |K -) dγ + K⊂Ω 1 #K σ,σ ∈K (v σ -v σ ) Φ K σ (u h ) -Φ K,Gal σ (u h ) + Γ⊂∂Ω 1 #Γ σ,σ ∈Γ (v σ -v σ )(Φ Γ σ u h , u b ) -Φ Gal,Γ σ (u h , u b ) (22a)
where

Φ K,Gal σ (u h ) = - K ∇ϕ σ • f (u h ) dx + ∂K ϕ σ fn (u h |K , u h |K -) dγ, Φ Γ,Gal σ (u h , u b ) = Γ ϕ σ fn (u h , u b ) -f (u h ) • n dγ. (22b) Note that σ∈K Φ K,Gal σ (u h ) = ∂K fn (u h |K , u h |K -) dγ and σ∈Γ Φ Γ,Gal σ (u h , u b ) = fn (u h , u b ) -f (u h ) • n dγ.
The proof of (22a) is given in appendix C. The relation (22a) is instrumental in proving the following results. The first one is proved in [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF], and is a generalisation of the classical Lax-Wendroff theorem.

Theorem 3.2. Assume the family of meshes T = (T h ) is shape regular. We assume that the residuals {Φ K σ } σ∈K , for K an element or a boundary element of T h , satisfy:

• For any M ∈ R + , there exists a constant C which depends only on the family of meshes T h and M such that for any

u h ∈ V h with ||u h || ∞ ≤ M , then Φ K σ (u h |K ) ≤ C σ,σ ∈K |u h σ -u h σ |
• The conservation relations (7a) and (7b) with quadrature formulas satisfying assumption 3.1.

Then if there exists a constant C max such that the solutions of the scheme (8

) satisfy ||u h || ∞ ≤ C max and a function v ∈ L 2 (Ω) such that (u h ) h or at least a sub-sequence converges to v in L 2 (Ω), then v is a weak solution of (1)
Proof. The proof can be found in [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF], it uses (22a) and some adaptation of the ideas of [START_REF] Kröner | A Lax-Wendroff type theorem for upwind finite volume schemes in 2-d. East-West[END_REF]. One of the key arguments comes from the consistency of the flux f as well as ( 21) applied to f

Another consequence of (22a) is the following result on entropy inequalities: Proposition 3.3. Let (U, g) be a entropy-flux pair for (1) and ĝn be a numerical entropy flux consistent with g • n. Assume that the residuals satisfy: for any element K,

σ∈K V σ , Φ K σ ≥ ∂K ĝn (V h |K , V h |K -) dγ (23a)
and for any boundary edge e,

σ∈e V σ , Φ e σ ≥ e ĝn (V h Γ , V b ) -g(V h ) • n dγ. ( 23b
)
Then, under the assumptions of theorem 3.2, the limit weak solution also satisfies the following entropy inequality: for any ϕ ∈ C 1 (Ω), ϕ ≥ 0,

- Ω ∇ϕ • g(V ) dx + ∂Ω - ϕ g(V b ) • n dγ ≤ 0.
Proof. The proof is similar to that of theorem 3.2 once we are writing the conservative variables in term of the entropy variable Another consequence of (22a) is the following condition under which one can guarantee to have a k + 1-th order accurate scheme. We first introduce the (weak) truncation error

E u h , ϕ = σ∈S h ϕ σ K⊂Ω,σ∈K Φ K σ + Γ⊂∂Ω,σ∈Γ Φ Γ σ . ( 24
)
Proposition 3.4. Assuming that:

1. The solution of the steady problem u is smooth enough, 2. The approximation π h (u) of u is accurate with order k + 1, and f Lipschitz continuous, 3. The quadrature formula (5) and (6) are of order k + 1: for any ψ of class at least C k+1 , for any element K and face f ,

K ψ(x) dx = K ψ(x) dx + O(h d+k+1 ), f ψ dγ = f ψ dγ + O(h d+k ),
4. The residuals, computed with the interpolant π h (u) of the solution, are such that for any element K and boundary element Γ

Φ K σ (π h (u)) = O h k+d , Φ Γ σ (π h (u)) = O h k+d-1 (25)
then the truncation error satisfies the following relation

|E π h (u), ϕ | ≤ C(f , u) ||ϕ|| H 1 (Ω) h k+1 ,
with C a constant which depends only on f , and ||u|| ∞ .

The proof can be found in appendix C.

Local conservation. In [START_REF] Abgrall | Some remarks about conservation for residual distribution schemes[END_REF], we show that schemes satisfying the conservation requirements (7a) and (7b) are indeed locally conservative. More precisely, for each degree of freedom, one can define a control volume Cσ of polygonal type, and numerical flux fn σσ (u h ) that depend on the values u σ in a neighborhood of C σ and the local normal of C σ ∪ C σ which is planar. This set of connection defines a dual graph. By construction, we have fn σσ (u n ) = -fn σ σ (u h ), which ensure local conservation. Details can be found in [START_REF] Abgrall | Some remarks about conservation for residual distribution schemes[END_REF], but we give two examples in appendix B. This means that (8) can be rewritten as

σ ∈Nσ fn σσ (u h ) = 0
where N σ is the set of degrees of freedom connected to σ by the dual graph defined from the cells.

Construction of entropy conservative schemes

In the rest of the paper, V h ∈ V h represents the approximation of the entropy variable V = ∇ u U (u) in V h , and since, the mapping u → V (u) is one-to-one, u h := u(V h ) here. Hence in general it is not a polynomial.

Starting from a scheme {Φ σ }, we show in this section how to construct a new scheme, with residuals

{Φ σ } such that σ∈K V σ , Φ σ = ∂K ĝn (V h |K , V h |K -) dγ. ( 26
)
This relation is motivated by Proposition 3.3 where the inequality is replaced by an equality, so that combined with Theorem 3.2, we know that if the scheme converges, the limit solution will satisfy two conservation relations: one on u and one on the entropy. The equality, together with the remark of the appendix B shows that we get indeed local conservation of the entropy, too.

Following [START_REF] Abgrall | Staggered grid residual distribution scheme for lagrangian hydrodynamics[END_REF][START_REF] Abgrall | A high-order nonconservative approach for hyperbolic equations in fluid dynamics[END_REF],the idea is to write:

Φ σ = Φ σ + r σ ( 27
)
such that conservation is still guaranteed

σ∈K Φ K σ + r σ = ∂K fn (V h |K , V h |K -) dγ, i.e. σ∈K r σ = 0, (28a)
and ( 26) holds true, i.e.

σ∈K V σ , r σ = ∂K ĝn (V h |K , V h |K -) dγ - σ∈K V σ , Φ σ := E (28b)
The two relations (28) defines a linear system with always at least two unknowns: the solution which is always valid is

r σ = α V σ -V with V = 1 #K σ∈K V σ . (29a) Since σ∈K V σ , V σ -V = σ∈K V σ -V , V σ -V + V , σ∈K (V σ -V ) = σ∈K V σ -V , V σ -V = σ∈K (V σ -V ) 2 since σ∈K (V σ -V ) = 0 by construction. Thus we get α = E σ∈K (V σ -V ) 2 . ( 29b
)
Remark 3.5.

1. Since (28b) means that the new set of residuals satisfy a conservation relation for the entropy on each element, from [START_REF] Abgrall | Some remarks about conservation for residual distribution schemes[END_REF] the new scheme will be locally conservative for the variable u and the entropy.

2. Relation (29b) becomes singular for a constant state. In practice, we use

α = E σ∈K (V σ -V ) 2 + ε ( 30
)
and we have taken ε = 10 -20 .

3. Interpretation of [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF]. Since α has a priori no sign, depending on the local entropy production of the original scheme, one adds or removes the right amount to be entropy conservative.

In fact adding α(V σ -V ) = α #K σ ∈K (V σ -V σ ) to
Now, let us have a look at the accuracy of such a scheme. Let us remind we are looking at steady problems. According to the accuracy conditions [START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF], we assume Φ σ = O(h k+d ). Similarly, since the numerical flux ĝ is Lipschitz continuous and consistent, we also have point-wise, for a smooth steady exact solution V , div g(V ) = 0, so that

0 = K div g(V ) dx = ∂K g(V ) • n dγ and then ∂K ĝn (V h |K , V h |K -) dγ = ∂K ĝn (V h |K , V h |K -) dγ - ∂K g(V ) • n dγ = ∂K ĝn (V h |K , V h |K -) -g(V ) • n dγ + ∂K g(V ) • n dγ - ∂K g(V ) • n dγ = O(h k+d ) (31)
provided the quadrature formula on the boundary is at least of order k. This is the case by assumption. In general, we cannot have more than O(h k+d ) because even if the integration is exact, we cannot have

∂K ĝn (V h |K , V h |K -) -g(V ) dγ = O(h k+d+1 ).
Hence, E = O(h k+d ). However, for smooth problems, we have

α = O(h k-2+d ),
so that the correction r σ is only O(h k-1+d ) a priori: we loose one order. Thus we get Proposition 3.6. The scheme defined by [START_REF] Hughes | A new finite element formulation for CFD: I. symmetric forms of the compressible Euler and Navier-Stokes equations and the second law of thermodynamics[END_REF], ( 29) is localy entropy conservative and at least of order k at steady state.

Up to now, we have shown, starting form a general additional conservation constraint, how to modify a scheme so that the modified scheme will satisfy one additional conservation relation, and what we see is that a priori we loose one order of accuracy. This is a general statement that assumes there is no particular connection between the 'old' conservation relation and the new one. However, there is a connection for the entropy flux, it is given by the following relation, see [START_REF] Tadmor | Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems[END_REF]: the entropy flux g is related to the flux f by the relation

g(V ) = V , f (V ) -θ(V )
In addition, the potential satisfies

∇ V θ(V ) = f (V ).
Using this, and a proper definition of the entropy numerical flux, it is possible to do better. We explore this question for the discontinuous Galerkin method, and the non linear RD scheme. The discussion on the discontinuous Galerkin method will also embed the continuous Galerkin approximation.

Discontinuous Galerkin methods

In the case of the discontinuous Galerkin methods, we get

E = -∇V h • f (V h ) dx + ∂K V h , f n (V h |K , V h |K -) dγ - ∂K ĝn (V h |K , V h |K -) dγ.
The precise definition of ĝn will be done later in this section. Since point-wise1 , we have

∇V h • f (V h ) = div ( V h , f (V h ) ) -V h • ∇f (V h ) and V h • div f (V h ) = div g(V h
), we can write:

E = - K ∇V h • f (V h ) dx + ∂K V h , fn (V h |K , V h |K -) dγ - ∂K ĝn (V h |K , V h |K -) dγ = - K div V h , f (V h ) dx + K V h • div f (V h ) dx + ∂K V h , fn (V h |K , V h |K -) dγ - ∂K ĝn (V h |K , V h |K -) dγ = - K div V h , f (V h ) dx + K div g(V h ) dx + ∂K V h , fn (V h |K , V h |K -) dγ - ∂K ĝn (V h |K , V h |K -) dγ = - K div ( V h , f (V h ) ) dx + K div ( V h , f (V h ) ) dx (I) + K div g(V h ) dx - K div g(V h ) dx (II) - K div ( V h , f (V h ) -g(V h ) dx + ∂K V h , fn (V h |K , V h |K -) dγ -ĝn (V h |K , V h |K -) dγ = (I) + (II) + ∂K V h f n (V h ) -g n (V h ) dγ - ∂K V h • f n (V h ) -g n (V h ) dγ III + ∂K V h • fn (V h |K , V h |K -) -f n (V h ) -ĝn (V h |K , V h |K -) -g n (V h ) dγ (32)
Since V h is a polynomial, if the quadrature formulas are of order k, we get that I, II and III in [START_REF] Abgrall | Some remarks about conservation for residual distribution schemes[END_REF] are

I = O(h d+k+1 ), II = O(h d+k ), III = O(h d+k ).
Let us have a look at the last term. If we define the entropy flux to be

ĝn (V h , V h,-) = {V h }, fn (V | Kp h , V h |K -) -θ({V h }), (33)
we have

-V h , f (V h ) + V h , fn (V h |K , V h |K -) + g(V h ) • n -ĝn (V h |K , V h |K -) = V h -{V h }, fn (V h |K -, V h |K -) + θ({V h }) -θ(V h ) • n = V h -{V h }, fn (V h |K , V h |K -) -∇ v θ( V ), V h -{V h } = V h -{V h }, fn (V h |K , V h |K -) -∇ V θ( V )
where V = t{V h } + (1t)V by abuse of language 2 . Since ∇ V θ = f , we have indeed

-V h , f (V h ) + V h , fn (V h |K , V h |K -) + g(V h ) • n -ĝn (V h |K , V h |K -) = V h -{V h }, fn (V h |K , V h |K -) -f ( V ) = O(h 2(k+1) )
because the numerical flux is Lipschitz continuous. Hence

∂K -V h , f (V h ) + V h , fn (V h |K , V h |K -) + g(V h ) • n -ĝn (V h |K -, V h |K -) dγ = O(h d+2k+1 )
provided suitable quadrature formulas. In the end we get Proposition 3.7. The scheme defined by ( 27), ( 29) is localy entropy conservative. In addition

E = σ∈K V σ , Φ σ -∂K ĝn (V h |K , V h |K -) dγ, with ĝn (V h |K , V h |K -) = {V }, f -θ({V h }) satisfies E = O(h d+k+1 ) (34)
if the quadrature formula are of order k for the surface/volume integrals and k + 1 for the boundary, so that the scheme (27)-eqrefent is stil of order k + 1 at steady state.

Remark 3.8. 1. We do not use the fact we are solving a steady problem: the same would be true for the semi-discrete unsteady problem.

2. We have only used approximation properties: the same would be true when non linear limiting is applied, provided the formal order is kept for sommth solutions

Non linear RD schemes

Let us recall that for the non linear RD schemes defined in appendix A, the accuracy condition is that the boundary integrals are of order k -1, since then

∂K fn (V h |K , V h |K -) dγ = O(h k+d ).
2 We write θ({V h }) -θ(V h ) component by component. For the component #l, we get

θ l ({V h }) -θ l (V h ) = ∂θ l ∂V (t l {V h } + (1 -t l )V ), {V h } -V , t l ∈ [0, 1] Hence V h -{V h }, fn(V h |K , V h |K -) -∇ V θ( V ) means, by abuse of language p l=1 V h -{V h } l fn(V h |K , V h |K -) l - ∂θ l ∂V (t l {V h } + (1 -t l )V ) .
If ĝn is a Lipschitz continuous flux that is consistent with the entropy flux, we also have

∂K ĝn (V h |K , V h |K -) dγ = O(h k+d ).
Let us look at E in more detail when Φ σ = β σ Φ, σ∈K β σ = Id. We have, defining again ĝ = {V }, f -θ({V h }) and V = σ∈K β σ V σ , and taking into account we assume a smooth steady solution

E = ∂K ĝn (V h |K , V h |K -) dγ - ∂K V, fn (V h |K -, V h |K -) dγ = ∂K ĝn (V h |K , V h |K -) -g(V ) dγ - ∂K V, fn (V h |K , V h |K -) -V, f (V ) dγ.
We look at each term. First, using Taylor formula, and the fact that

∇ V θ(V ) = f (V ), ĝn (V h |K , V h |K -) -g(V ) = {V h }, fn (V h |K , V h |K -) -θ({V h } -V , f (V ) + θ(V ) = {V h } -V , f (V ) + {V h }, fn (V h |K , V h |K -) -f (V ) + ∇ V θ(V ), V -{V h } + O (V -{V h }) 2 = {V h }, fn (V h |K , V h |K -) -f (V ) + O (V -{V h }) 2 so that ĝn (V h |K , V h |K -) -g(V ) -V , fn (V h |K -, V h |K -) -f (V ) = {V h } -V , fn (V h |K , V h |K -) -f (V ) + O (V -{V h }) 2 = O(h) × O(h k+1 )
if the boundary quadrature formula is of order k. We have Proposition 3.9. The scheme defined by [START_REF] Hughes | A new finite element formulation for CFD: I. symmetric forms of the compressible Euler and Navier-Stokes equations and the second law of thermodynamics[END_REF], ( 29) is localy entropy conservative. In addition,

E = σ∈K V σ , Φ σ -∂K ĝn (V h |K , V h |K -) dγ, with ĝn (V h |K , V h |K -) = {V }, f -θ({V h }) satisfies E = O(h k+d+1 ) ( 35
)
if the boundary quadrature formula is of order k.

Remark 3.10. Here we must use the fact we are solving a steady problem.

Construction of entropy stable schemes

Starting form the residuals {Φ K σ } K∈T h ,σ∈K , we want to construct a scheme {Φ K σ } K∈T h ,σ∈K , that is entropy stable, i.e., for any K,

σ∈K V σ , Φ K σ (V h ) ≥ ∂K ĝn (V h |K , V h |K -) dγ.
Starting from the previous construction, and dropping the superscript K because we will consider only quantities in K, the idea is to write the new residuals as:

Φ σ (V h ) = Φ σ + r σ + Ψ σ ( 36
)
where the r σ are defined by ( 29) and the Ψ σ satisfy σ∈K Ψ σ = 0 and σ∈K V σ , Ψ σ (V h ) ≥ 0 without violating the accuracy condition [START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF]. Let us consider these two expressions for Ψ σ (V h ) which obviously satisfy the conservation requirement

σ∈K Ψ σ (V h ) = 0 since σ∈K ϕ σ = 1:
1. With jumps: We define

Ψ σ = θh 2 K ∂K [∇ x ϕ σ ] • [∇ x V h ] dγ (37)
We have

σ∈K V σ , Ψ σ (V h ) = θ h 2 K e edge of K K [∇ x V h ] 2 dγ ≥ 0
for any θ > 0. It is also easy to check that if V h is an approximation of order k of a smooth V , we have, provided the quadrature formula is of order k that

Ψ σ = O(h k+d )
so that the accuracy requirement are met.

2. With streamline: we define

Ψ σ = θh K K ∇ V f (V h ) • ∇ x ϕ σ τ K ∇ V f (V h ) • ∇ x V h dx. (38)
Clearly,

σ∈K V σ , Ψ σ (V h ) = θh K K ∇ V f (V h ) • ∇ x V σ τ K ∇ V f (V h ) • ∇ x V h dx ≥ 0
for any θ > 0 and any τ > 0 In practice, we take τ K as the scalar/matrix defined by

τ -1 = σ ∈K K - σ , K σ = 1 |K| K ∇ u f • ∇ x ϕ σ dx, K + σ = max(K σ , 0), K - σ = min(K σ , 0).
It is also easy to check that if V h is an approximation of order k of a smooth V , we have, provided the quadrature formula is of order k that Ψ σ = O(h k+d ) so that the accuracy requirement are met.

In practice, the quadrature formula involved in the edge integrals and the surface integrals need to be such that:

• For e [∇ x V h ] 2 dγ: if V h is not constant, e [∇V h ] 2 dγ > 0. • For K ∇ V f (V h ) • ∇ x V h τ K ∇ V f (V h ) • ∇ x V h dx: if ∇ V f (V h ) • ∇ x V h is not identicaly 0, then K ∇ V f (V h ) • ∇ x V σ τ K ∇ V f (V h ) • ∇ x V h dx > 0.
In [START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF], we have provided minimal conditions on these quadrature formula, such that this conditions are met. Note that these quadrature formulas may be non-consistent: the only things we need is to keep the accuracy, and have a strictly positive entropy production.

Numerical illustrations

In this section, we illustrate the behaviour of the methods on non linear examples. We first study the conservation issues, for mass and entropy, and then show the behavior of the scheme on a Burger-like problem.

Conservation

In Ω = [-20, 20] × [-20, 20] we look at the problem

∂u ∂t + ∂ √ u ∂x + ∂u ∂y = 0 (39a) with u 0 (x, y) = 1 if r = (x + 5) 2 + (y + 5) 2 > 0.5 1 + 2 sin π 4 (1 -2r) else.
u(x, y, t) = u 0 (x, y) for t ≥ 0, and x ∈ ∂Ω.

(39b)

We integrate for t ∈ [0, 5]: at time t = 5, the non constant part of the solution has not reached the boundary, this is why we have taken such a large domain. We want to avoid any interference with the boundary conditions. The flux is non polynomial, so that any standard quadrature rule will never be exact. Following the notations of (8), we integrate in time with an Euler forward method

u n+1 σ = u n σ - ∆t C σ K⊂Ω,σ∈K Φ K σ (u n,h ) + Γ⊂∂Ω,σ∈Γ Φ Γ σ (u n,h ) , (40)
where C σ is a dual control volume. This method is not accurate in time, but we can check numerically the conservation of mass,

M = σ∈Ω |C σ |u σ .
An entropy for (39a) is

U (u) = u 2 2 , the entropy flux is g(u) = ( u 2 2 , 1 3 u √ u). Since V (u n ), u n+1 -u n ≤ U (u n+1 ) -U (u n )
we will not get conservation of entropy, but we will check that, the entropy residuals defined by ( 27) Φ σ are such that σ∈Ω K⊂Ω,σ∈K

V σ , Φ K σ (u n,h ) + Γ⊂∂Ω,σ∈Γ V σ , Φ Γ σ (u n,h ) = 0. ( 41
)
Up to our knowledge, any of the so-called entropy conservative schemes show rigorous conservation of the spatial discretisation terms only. To get also entropy conservation at the time-discrete level, it seems that the scheme must be implicit in time, see [START_REF] Hiltebrand | Entropy stable shock capturing spacetime discontinuous Galerkin schemes for systems of conservation laws[END_REF]. This is the reason why we are not showing

E = Ω U (u h ) dx
but how well (41) is fullfiled in figure 1. We see that the mass is always conserved up to machine accuracy, while the correction has a clear effect on the spatial conservation of entropy. This was shown for the pure Galerkin scheme.

Accuracy

In this section we want to show that the correction does not spoil the accuracy. This is done on the following test case

: in Ω = [0, 1] × [0, 1], ∂ sinh u ∂x + ∂u ∂y = 0 (42a)
with the boundary condition, for y = 0,

u(x, 0) = x - 1 2 (42b)
The solution is regular and computed by the method of characteristics. We have chosen the exponential flux to make sure that none of the standard quadrature formula is exact. The quadrature formula are:

• Linear approximation. the quadrature point in the triangle K is its centroid, and the weight is 1. On the boundary, we have chosen the Gaussian formula associated to the Legendre polynomials of degree 2: there are two points of weight 1 2 that are the image of ± 1 √ 3 in the mappings form [-1, 1] to any of the edges of K.

• Quadratic approximation. The surface integrals are computed with 7 quadrature points, the formula is exact for degree 5 The integral on the edges are evaluated by the same Gaussian formula as before. We have used the scheme 36 where Φ σ is obtained by the Galerkin formulation using the above quadrature and Ψ σ is defined by [START_REF] Abgrall | High order schemes for hyperbolic problems using globally continuous approximation and avoiding mass matrices[END_REF] with θ = 0.01. The boundary conditions are implemented with a local Lax Friedrich numerical scheme. The errors are given in table 1 and2.

-
As expected the optimal orders are reached. We can also observe the improved accuracy for a given number of freedom, between second and third order.

Behavior on discontinuous solutions

The non linear example is an adaptation of the Burgers equation:

∂u ∂y + ∂ sinh u ∂x = 0 if x ∈ [0, 1] 2 u(x, y) = 1 2 -x on the inflow boundary. ( 43
)
If the y-coordinate corresponds to time, we are back to a more standard formulation. The exact solution consists in a fan that merges into a shock. We have performed the simulations using a quadratic approximation. The mesh (i.e. all the degrees of freedom) is displayed in figure 2. The numerical solutions are obtained by several methods:

• The SUPG scheme with different values of θ in [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: second order case[END_REF].

• The Galerkin scheme with jumps and different values of θ in (37)

• The non-linear RD scheme where the first order scheme is the Rusanov one (see appendix A)

• The non-linear RD scheme with jump filtering, [START_REF] Ray | Entropy stable schemes for compressible Euler equations[END_REF], the jumps being tuned by θ

• The non-linear RD scheme with streamline filtering, [START_REF] Svärd | Entropy-stable schemes for the Euler equations with far-field and wall boundary conditions[END_REF], the jumps being tuned by θ

These schemes are nicknamed respectively as SUPG-θ, GalJ-θ, RD, RD-J-θ and RD-S-θ. We also have made a second set of simulations where the schemes are modified by adding the entropy correction. More precisely, this is done after the evaluation of the Galerkin term for the SUPG-θ and GALJ-θ schemes, and after the nonlinear RD step of the RD-J-θ and RD-S-θ steps, i.e. before the filtering step in all cases, as described in section 3.4. These schemes are nicknamed as SUPG-E-θ, GalJ-E-θ, RD, RD-J-E-θ and RD-S-E-θ respectively. Since the numerical flux is not polynomial, the Galerkin step is only approximate: the entropy correction is a priori acting. We have used the same quadrature points as in section 4.2.

Our goal is to see if we can lower the parameter θ in the filtering term because we have a fine tuning of the local entropy production. All the simulations are done with the Euler forward time stepping, CFL=0.3. The plots use the same isolines between 2 and -2, 100 isolines: we can see the numerical oscillations, if any. From figure 3, we see that the SUPG scheme without correction and θ S = 0.1 provides very good results in the fan but is almost unstable in the shock (the calculation has been stopped since it is meaningless). When we lower the parameter θ to 0.01, the shock does not improve, and some wiggles exist in the fan, so we are not filtering enough. A relatively good compromise is reached for θ = 0.02. This is in contrast with the GalJ-0.01 result that are comparable to the SUPG-0.02 ones (with smaller oscillations at the shock, by the way).

From figure 4, we see that θ = 0 leads to wiggles in the fan. They are due to spurious modes, and not a stability problem, see [START_REF] Abgrall | Essentially non oscillatory residual distribution schemes for hyperbolic problems[END_REF][START_REF] Abgrall | About non linear stabilization for scalar hyperbolic problems[END_REF]. We have tried several values of θ for the streamline and the jump term, it seems that the best results are obtained for θ = 0.1, both for the fan and the discontinuity (remember that we have 100 isolines between -2 and 2, so 25 between -0.5 and 0.5.

From figure 5, we see that Gal-E-0.01 behaves much better that SUPG-E-0.01. This is an other indication (comparing to figure 3-(b) that the streamline term is not suffisant. The solution is much better with θ = 0.02, see figure 5-(c), but it is similar to 3-(c). We have also experienced similar disapointments with the unsteady version of these schemes in [START_REF] Abgrall | High order schemes for hyperbolic problems using globally continuous approximation and avoiding mass matrices[END_REF] where the jump term seems to filter out more efficiently. Comparing figures 4-(c) and 5-(d), we see that the quality of the fan has improved with the entropy correction, while, of course, the shock becomes wiggly. If we decrease the value of θ, then the solution becomes similar to 5-(b).

Conclusion

We have developed a general technique to guaranty that the discretisation of a conservative hyperbolic problem is consistent with an additional conservation relation. This is a generalisation of the work [START_REF] Abgrall | A high-order nonconservative approach for hyperbolic equations in fluid dynamics[END_REF] and a sequel of [START_REF] Abgrall | Some remarks about conservation for residual distribution schemes[END_REF]. This construction uses in depth a reinterpretation of conservative schemes as Residual distribution schemes. Of course, the modifications will depend on the problem. Aiming at illustrating this idea, this is why we have focussed on an entropy equality and have constructed explicit modifications of the original discretisation that guaranty a consistency with the original PDE and the additional entropy inequality. We also have shown that independently of the choice of the entropy flux, one lose one order of accuracy. However, in the case of discontinuous Galerkin schemes, and stabilized residual distribution schemes, we can provide an explicit form of the entropy flux that does not lead to a degradation of the accuracy. To achieve this goal, we do not need any particular structure of the quadrature formula, in particular they do not need to be exact nor to achieve the summation-by-part property.

We have tested the method on scalar problem, the formula we have given also work for systems. We have shown, numerically, at least for finite element like methods, that the optimal accuracy is kept. We have also checked that the method is globally conservative for the original conserved variable and the entropy. Using the results of [START_REF] Abgrall | Some remarks about conservation for residual distribution schemes[END_REF], it is easy to see that, at least for the steady case, that these schemes can be reformulated as finite volume schemes. In the case of entropy stable scheme, we can get an explicit form of the numerical flux for the conserved variable and the entropy.

We have also shown that adding this entropy correction can help to lower the amount of artificial diffusion that is needed to filter spurious modes for the Galerkin method and the non linear residual distribution ones. Our study also confirm that the stabilisation by jump is more efficient than the streamline diffusion one. One reason is probaly that the jump stabilisation allow inter-element exchange of informations, contrarily to the streamline one.

The emphasis of this paper is put on the steady case, but the unsteady state is similar, see [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: second order case[END_REF] and [START_REF] Abgrall | High order schemes for hyperbolic problems using globally continuous approximation and avoiding mass matrices[END_REF][START_REF] Abgrall | High-order residual distribution scheme for the timedependent Euler equations of fluid dynamics[END_REF]. However, if the unsteady case can be considered as well, this is to the price of having implicit formulations, in the spirit of [START_REF] Abgrall | A high-order nonconservative approach for hyperbolic equations in fluid dynamics[END_REF]. This issue will be considered in future work. Future work will also deal with the question of how to have an entropy conservative scheme in smooth regions, and how to guaranty non oscillatory behaviour in discontinuous ones

We conclude this work in noticing that one of the ingredients of the method is that any element has at least 2 vertices: the conservation relation translates by a linear relations on the residuals, and the modification by another one. Having more than 2 degrees of freedom, one solution is a priori possible. As a consequence, it could be possible, using the same kind of technique, to guaranty more that one additional conservation relation, this will also be considered for future work.

A Construction of the non linear stabilisation for RD schemes

Here we consider a globally continuous approximation:

u h ∈ V h ∩ C 0 (Ω).
Consider one element K. Since there is no ambiguity, the drop, for the residuals, any reference to K in the following. The total residual is defined by

Φ(u h ) = ∂K f (u h ) • n dγ,
and we assume to have monotone residuals {Φ L σ } σ∈K . By this we mean

Φ L σ (u h ) = σ ∈K c L σσ (u σ -u σ )
with c σσ ≥ 0 that also satisfies

σ∈K Φ L σ = Φ.
It can easily be shown that the condition c σσ ≥ 0 garanties that the scheme is monotone under a CFL like condition. One example is given by the Rusanov residuals:

Φ Rus σ (u h ) = - K ∇ϕ σ • f (u h ) dx + ∂K ϕ σ f (u h ) • n dγ + α(u σ -ū),
where ū is the arithmetic average of of the u σ s on K and α satisfies:

α ≥ #K max σ,σ ∈K K ϕ σ ∇ u f (u h ) dx .
Here #K is the number of degrees of freedom in K. Indeed, this residual can be rewritten as

Φ Rus σ (u h ) = σ ∈K c σσ (u σ -u σ ) with c σσ = K ϕ σ • ∇ u f (u h ) dx + α #K .
Under the condition above, c σσ ≥ 0 and hence we have a maximum principle.

The coefficients β σ introduced in the relations ( 18) and ( 19) are defined by:

β σ = max(0, Φ L σ Φ ) σ ∈K max(0, Φ L σ Φ )
.

and can be shown to be always defined, to guaranty a local maximum principle for ( 18) and [START_REF] Ray | Entropy stable schemes for compressible Euler equations[END_REF], see [START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF].

Remark A.1 (About the coefficients c L σσ ). All the examples of monotone residual we are aware of are such that for linear problems, te c L σσ are independant of u h . Then one can show that for any σ ∈ K,

σ ∈K c L σσ -c L σ σ = K ∇f (u h ) • ∇ϕ σ dx.
This relation implies the conservation relation (7a).

B Examples of finite volume writing of Residual distribution schemes

This part is taken from [START_REF] Abgrall | Some remarks about conservation for residual distribution schemes[END_REF]. As a motivation, we first show that any finite volume scheme can be written as a residual distribution sccheme. Then we give two example of the converse. The general statement can be found in [START_REF] Abgrall | Some remarks about conservation for residual distribution schemes[END_REF], it applies also to discontinuous representations, see the same reference

B.1 Finite volume as Residual distribution schemes

The notations are defined in Figure 7. We specialize ourselves to the case of triangular elements, but exactly the same arguments can be given for more general elements, provided a conformal approximation space can be constructed. This is the case for triangle elements, and we can take k = 1.

The control volumes in this case are defined as the median cell, see figure 7. We concentrate on the approximation of div f , see equation [START_REF] Chen | Entropy stable high order discontinuous Galerkin methods with suitable quadrature rules for hyperbolic conservation laws[END_REF]. Since the boundary of C σ is a closed polygon, the scaled outward normals n γ to ∂C σ sum up to 0:

γ⊂∂Cσ n γ = 0
where γ is any of the segment included in ∂C σ , such as P G on Figure 7. Hence

γ⊂∂Cσ fnγ (u σ , u -) = γ⊂∂Cσ fnγ (u σ , u -) - γ⊂∂Cσ n γ • f (u σ ) = K,σ∈K γ⊂∂Cσ∩K fnγ (u σ , u -) -f (u σ ) • n γ
To make things explicit, in K, the internal boundaries are P G, QG and RG, and those around σ ≡ 1 are P G and RG. We set

Φ K σ (u h ) = γ⊂∂Cσ∩K fnγ (u σ , u -) -f (u σ ) • n γ = γ⊂∂(Cσ∩K) fnγ (u σ , u -). ( 44
)
The last relation uses the consistency of the flux and the fact that C σ ∩ K is a closed polygon. The quantity Φ K σ (u h ) is the normal flux on C σ ∩ K. If now we sum up these three quantities and get:

σ∈K Φ K σ (u h ) = fn12 (u 1 , u 2 ) -fn13 (u 1 , u 3 ) -f (u 1 ) • n 12 + f (u 1 ) • n 31 + fn23 (u 2 , u 3 ) -fn12 (u 2 , u 1 ) + f (u 2 ) • n 12 -f (u 2 ) • n 23 + -fn23 (u 3 , u 2 ) + fn31 (u 3 , u 1 ) -f (u 3 ) • n 23 + f (u 3 ) • n 31 = f (u 1 ) • n 12 -n 31 + f (u 2 ) • -n 23 + n 31 + f (u 3 ) • n 31 -n 23 = f (u 1 ) • n 1 2 + f (u 2 ) • n 2 2 + f (u 3 ) • n 3 2
where n j is the scaled inward normal of the edge opposite to vertex σ j , i.e. twice the gradient of the P 1 basis function ϕ σj associated to this degree of freedom. Thus, we can reinterpret the sum as the boundary integral of the Lagrange interpolant of the flux. The finite volume scheme is then a residual distribution scheme with residual defined by (44) and a total residual defined by

Φ K := ∂K f h • n, f h = σ∈K f (u σ )ϕ σ . (45)
Let K be a fixed triangle. We are given a set of residues {Φ K σ } σ∈K , our aim here is to define a flux function such that relations similar to (44) hold true. We explicitly give the formula for P 1 and P 2 interpolant.

B.2 General case

One can deal with the general case, i.e when K is a polytope contained in R d with degrees of freedoms on the boundary of K. The set S is the set of degrees of freedom. We consider a triangulation T K of K whose vertices are exactly the elements of S. Choosing an orientation of K, it is propagated on T K : the edges are oriented. The problem is to find quantities fn σ,σ for any edge [σ, σ ] of T K such that:

Φ σ = edges [σ,σ ] fn σ,σ + f b σ (46a) with fn σ,σ = -fn σ ,σ (46b)
and

f b σ = ∂K ϕ σ fn (u h , u h -) dγ. (46c)
The control volumes will be defined by their normals so that we get consistency. The discontinuous nature of u h , if any, is handled through the numerical flux fn . Note that (46b) implies the conservation relation

σ∈K Φ σ = σ∈K f b σ , (46d)
so that we rewrite (46a), introducint Ψ σ = Φ σ -f b σ as:

Ψ σ = edges [σ,σ ] fn σ,σ (46e)
The general P 1 case. A solution is

fn13 = Ψ 1 -Ψ 3 3 fn23 = Ψ 2 -Ψ 3 3 fn32 = Ψ 3 -Ψ 2 3
In order to describe the control volumes, we first have to make precise the normals n σ in that case. It is easy to see that in all the cases described above, we have

N σ = - n σ 2 .
Then a short calculation shows that  

n 12 n 23 n 31   = 1 6   n 1 -n 2 n 2 -n 3 n 3 -n 1   .
Using elementary geometry of the triangle, we see that these are the normals of the elements the dual mesh. For example, the normal n 12 is the normal of P G, see figure 7.

Relying more on the geometrical interpretation (once we know the control volumes), we can recover the same formula by elementary calculations, see [START_REF] Abgrall | On a class of high order schemes for hyperbolic problems[END_REF].

The general example of the P 2 approximation.

Using a similar method, we get (see figure 8 The result on the boundary term is similar since the boundary numerical flux is upwind and the boundary of Ω is not characteristic: only two types of boundary faces exists, the upwind and downwind ones. On the downwind faces, the boundary flux vanishes. On the upwind ones, we get the estimate for the Galerkin boundary residuals thanks to the same approximation argument.

The mesh is assumed to be regular: the number of elements (resp. edges) is O(h -d ) (resp. O(h -d+1 )). Let us assume [START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF]. Let v ∈ C 1 0 (Ω). Using (22a) for π h (u), where π h (u) -represents the interpolant of u on K -. We have, using

- Ω ∇π h (v) • f (u) dx + ∂Ω π h (v) F n (u, u b ) -f (u) • n dγ = 0, - Ω ∇π h (v) • f (π h (u)) dx + ∂Ω π h (v) fn (π h (u), u b ) -f (π h (u)) • n dγ = - Ω ∇π h (v) • f (π h (u)) -f (u) dx + ∂Ω π h (v) fn (π h (u), u b ) -fn (u, u b ) dγ - ∂Ω π h (v) f (π h (u)) -f (u) • n dγ = O(h k+1 )
since the flux on the boundary is the upwind flux F n , and using the approximation properties of π h (u). Then thanks to the regularity of the mesh, that π h (v) is the interpolant of a C 1 function and the previous estimates.

Figure 1 :

 1 Figure 1: Plot of (41) for the pure galerkin scheme, with and without the correction

--

  The centroid with weight 0.225, The points defined by the barycentric coordinates (a, b, b), a = 0.797426985353087, a + 2b = 1 and cyclic permutations of the coordinates. The weights are ω = 0.125939180544827. The points defined by the barycentric coordinates (a, b, b), a = 0.059715871789770, a + 2b = 1 and cyclic permutations of the coordinates. The weights are ω = 0.132394152788506.

Figure 2 :

 2 Figure 2: Mesh for the numerical experiments. All the Lagrange degrees of freedom are represented. From the representation, we evaluate the solution at the P 2 degrees of freedom and then plot the solution using a local P 1 representation in all the triangles represented on the figure.

  (a) SUPG-θ = 0.1 (b) SUPG-θ = 0.01 (c) SUPG-θ = 0.02 (d) GALJ-θ = 0.01

Figure 3 :

 3 Figure 3: Results without entropy correction, linear schemes

Figure 4 :

 4 Figure 4: Results without entropy correction, non linear schemes

  (a) GALJ-E-θ = 0.01 (b) SUPG-E-θ = 0.01 (c) SUPG-E-θ = 0.02 (d) RDJ-E-0.02

Figure 5 :

 5 Figure 5: Results with with the entropy term

Figure 6 :

 6 Figure 6: Notations for the finite volume schemes. On the left: definition of the control volume for the degree of freedom σ. The vertex σ plays the role of the vertex 1 on the left picture for the triangle K. The control volume C σ associated to σ = 1 is green on the right and corresponds to 1P GR on the left. The vectors n ij are normal to the internal edges scaled by the corresponding edge length

Figure 7 :

 7 Figure 7: Notations for the finite volume On the left: definition of the control volume for the degree of freedom σ. The vertex σ plays the role of the vertex 1 on the left picture for the triangle K. The control volume C σ associated to σ = 1 is green on the right and corresponds to 1P GR on the left. The vectors n ij are normal to the internal edges scaled by the corresponding edge length

Ψ 4 K

 4 Then we choose the boundary flux:f b σ = ∂K ϕ σ n dγ i.e. after having defined [v h ] = v h |Kv h|K -and chosen one orientation of the internal edges e ∈ E h , we get (22a).Proof of Proposition 3.4. We first show that Φ K,Gal σ (π h (u)) = O h k+d . Since u is regular enough, we have pointwise div f (u) = 0 on K, so that, by consistency of the flux,0 = -K ∇ϕ • f (u) dx + ∂K ϕ fn (u, u) dγ.This shows that -K ∇ϕ • f (u) dx + ∂K ϕ fn (u, u) dγ = O(h k+d+1 ). ∇ϕ σ • f (π h (u))f (u) dx + ∂K ϕ σ fn (π h (u), π h (u)) -fn (u, u) dγ + K ∇ϕ • f (u) dx -∂K ϕ fn (u, u) dγ = |K| × O(h -1 ) × O(h k+1 ) + |∂K| × O(1)O(h k+1 ) + O(h k+d+1 ) = O(h d ) × O(h -1 ) × O(h k+1 ) + O(h d-1 ) × O(1) × O(h k+1 ) + O(h k+d+1 ) = O(h d+k )because the flux is Lipschitz continuous and the mesh is regular.

E

  u h , ϕ = -Ω ∇π h (v) • f (π h (u)) dx + ∂Ω π h (v) fn (π h (u), u b )f (π h (u)) • n dγ + e∈E h e [π h (v)] fn (π h (u), π h (u) -) dγ + K⊂Ω 1 #K σ,σ ∈K (v σv σ ) Φ K σ (π h (u)) -Φ K,Gal σ (π h (u)) + Γ⊂∂Ω 1 #Γ σ,σ ∈Γ (v σv σ )(Φ Γ σ π h (u), u b ) -Φ Gal,Γ σ (π h (u), u b ) .

  e∈E h e [π h (v)] fn (π h (u), π h (u) -) dγ = O(h -d+1 ) × O(h d-1 ) × O(h k+1 ) × O(1) = O(h k+1 ), K⊂Ω 1 #K σ,σ ∈K (v σv σ ) Φ K σ (π h (u)) -Φ K,Gal σ (π h (u)) = O(h -d ) × O(h) × O(h k+d ) = O(h k+1 ),and similarlyΓ⊂∂Ω 1 #Γ σ,σ ∈Γ (v σ -v σ )(Φ Γ σ π h (u), u b )-Φ Gal,Γ σ (π h (u), u b ) = O(h -d+1 )×O(h)×O(h k+d-1 ) = O(h k+1 )

Table 1 :

 1 -5 2.11 3.30 10 -5 2.07 3.23 10 -4 1.81 1.76 10 -2 7905 5.04 10 -6 2.06 7.96 10 -6 2.04 8.75 10 -5 1.87 8.82 10 -3 31297 1.51 10 -6 1.74 2.76 10 -6 1.53 3.13 10 -5 1.5 Error for problem (42), linear approximation. -7 3.01 4.77 10 -7 2.92 6.28 10 -6 2.69 1.76 10 -2 31297 3.12 10 -8 2.88 6.35 10 -8 2.90 8.88 10 -7 2.81

	h	n dof s	L 1	slope	L 2	slope	L ∞	slope
	7.06 10 -2 3.53 10 -2 2017 525 2.12 10 h 9.16 10 -5 n dof s L 1	-slope	1.39 10 -4 L 2	-slope	1.14 10 -3 L ∞	-slope
	1.41 10 -1 7.06 10 -2 2017 525	1.35 10 -5 1.88 10 -6 2.84 3.62 10 -6 2.76 4.08 10 -5 2.39 -2.46 10 -5 -2.15 10 -4 -
	3.53 10 -2 7905	2.33 10

Table 2 :

 2 Error for problem (42), quadratic approximation.

Note that we are computing the divergence (in x) of the function V h , f (V h )

Acknowledgements

The author has been funded in part by the SNSF project 200021 153604 "High fidelity simulation for compressible materials". I would also like to thanks Anne Burbeau (CEA-DEN) for her critical reading of the first draft of this paper. Her input has hopefully helped to improve the readability of this paper. I also warmly thanks the three anonymous referees whose comments and criticisms have, I hope, improved the original version of the paper. The remaining mistakes are mine.

A personal note to end: though the author is curently editor in chief of this Journal, the EES system has been parametrised in such a way that I, of course, have absolutely no access to the identity of the referees, nor to the identity of the editor in charge of the paper. This is, of course, the only paper of this type. Hence I would like to thanks the Elsevier team for making this singular situation possible.

and get:

The normals are given by:

There is not uniqueness, and it is possible to construct different solutions to the problem. In what follows, we show another possible construction. We consider the set-up defined by Figure 8.

C Proof of (22a).

Proof of (22a). We start from [START_REF] Jiang | On a cell entropy inequality for discontinuous Galerkin methods[END_REF] which is multiplied by v σ , and these relations are added for each σ ∈ S. We get:

Permuting the sums on σ and K, then on σ and Γ, we get: We look at the first term, the second is done similarly. We have, introducing Φ K,Gal σ and #K the number of degrees of freedom in K,

Similarly, we have

Adding all the relations and using the assumption 3.1 on the quadrature formulas, we get:


A Laghrib
email: laghrib.amine@gmail.com

A Hakim

S Raghay
email: raghay@fstg-marrakech.ac.ma

M El Rhabi

A robust multi-frame super-resolution based on curvature registration and Second Order Variational regularization

Keywords: robust, super resolution, curvature registration, PDE, MAP estimator, image restoration, regularization. 2000 Mathematics Subject Classification: 35G20, 35F20, 35J05, 49J40, 49M25

Multiframe image super-resolution is a technique to obtain a high-resolution image by fusing a sequence of low-resolution ones. This paper deals with a new approach to robust super resolution based on regularization framework. Since registration is an important step that ensures the success of super resolution algorithms, must choose the most suitable method. We suggest a new algorithm specified at low resolution images with small deformations using fourth-order partial differential equations (PDE) regularization in the last step of super resolution. The deformations are not parametric and differs from one image to another. We use a curvature registration specially because image are slightly deformed. Experimental results show the robustness of the proposed method compared to classical super resolution methods.

Introduction

Multiframe super-resolution (SR) aim to improve the spatial resolution by combining detail of a set of low-resolution (LR) degraded images of the same scene to reconstruct a high resolution (HR) image of that scene. The degradation is modelled by different operators representing motion, blur, subsampling and additive noise. The SR process is possible if there exists a subpixel motion between the LR images. Therefore, each information about LR frame complete details of the original HR image. Robust Super-resolution (RSR) reconstruction is separate in three main steps: registration, interpolation and restoration. Registration is the process of estimating the deformation operator derived directly from the LR frame. In the interpolation step, the LR images are superimposed on to the HR image grid, while in the restoration step we remove noise and blur that is present in the degraded HR image.

The multiframe image super-resolution problem was first proposed by Tsai and Huang in E. Lee's 2003 (Tsai and[START_REF] Tsai | Multiframe image restoration and registration[END_REF], many methods have been proposed to overcome the ill-posedness of this problem using a prior distribution on the image as a regularization (E. [START_REF] Lee | Regularized adaptive high-resolution image reconstruction considering inaccurate subpixel registration[END_REF][START_REF] Milanfar | Super-Resolution Imaging (Digital Imaging and Computer Vision[END_REF][START_REF] Farsiu | Fast and robust multiframe super resolution[END_REF][START_REF] Sanches | A map estimation algorithm with gibbs prior using an iir filter[END_REF]. Regularization method is widely used to solve SR problem viewed as a maximum a posteriori approach. Popular Tikhonov regularizer was applied in SR problem in [START_REF] Nguyen | A note on tikhonov regularization of linear ill-posed problems[END_REF][START_REF] Lee | Regularized adaptive high-resolution image reconstruction considering inaccurate subpixel registration[END_REF], Ng et al. used the total variation (TV) regularizer for video SR, and illustrate their efficient results using PSNR criteria (M. Ng, n.d.). Another interesting regularisation was proposed by Farsiu et al., they introduced the bilateral filter with the L 1 norm and showed that L 1 norm is more robust against noise and misregistration (Sina [START_REF] Farsiu | Robust super-resolution[END_REF][START_REF] Farsiu | Fast and robust multiframe super resolution[END_REF], hence in fact the name of the method RSR.

A key issue that guarantees the success of the SR algorithm is the image registration technique. The motion in LR images may include simple translations (Sina Farsiu and Milanfar, OC-TOBER 2004;[START_REF] Term | Survey: Interpolation methods in medical image processing[END_REF][START_REF] Hiep Luong | An image interpolation scheme for repetitive structures[END_REF], affine [START_REF] Elad | Example-based single document image super-resolution: a global map approach with outlier rejection[END_REF][START_REF] Tsai | Multiframe image restoration and registration[END_REF] as well as non parametric transformations (Thomas M. Lehmann and Spitzer, NOVEM-BER 1999; [START_REF] Modersitzki | Fair: Flexible Algorithms for Image Registration (Fundamentals of Algorithms[END_REF][START_REF] Modersitzki | Numerical Methods for Image Registration[END_REF].

In this paper, we propose a curvature approach [START_REF] Modersitzki | Numerical Methods for Image Registration[END_REF] to register the LR frames estimating the deformations that differ between all LR images. The main of this technique is minimizing the distance between the estimation of the HR image and each upsampled LR image. In addition we propose a second PDE order regularisation in the debluring step [START_REF] Lysaker | fourth-order partial differential equation with applications to medical magnetic resonance images in space and time[END_REF][START_REF] Bergounioux | A second-order model for image denoising[END_REF], since we know they success to recover smoother surfaces.

Numerical results indicate the superiority of the new SR algorithm in noise and misregistration removing compared with TV regularization. To evaluate the robustness of our algorithm we use the PSNR criteria. This paper is organized as follows. Section 2 explains the main concepts of the SR. We justify the use of curvature registration and explain the use of second order PDE regularisation. Simulations and comparisons of our result with other classical methods are presented in Section 3.

Problem formulation

The observed images of a real scene are usually in low resolution due to some degradation operators. In practice, the acquired image is corrupted by noise, blur and decimation. In almost all cases the degradation is generated by inappropriate camera parameters or configuration.

In addition we have the effects of atmospheric turbulence. All these facts corrupt the resolution of the images, and improvement of resolution techniques is therefore required in those cases.

We assume that the LR images are taken under the same environmental conditions using the same sensor.

The unknown and desired high-resolution image X is related to the captured low-resolution images Y k (represented by a column vector of size M ) by this formulation [START_REF] Farsiu | Fast and robust multiframe super resolution[END_REF]:

Y k = DF k HX + E k ∀k = 1, 2, ..., n (2.1)
Where E k represents the additive noise for each image, which is assumed to be zero-mean

Gaussian distributed (with a standard deviation σ k ). H and D represent respectively the blurring operator of size N × N and decimation matrix of size M × N . F k is a geometric warp matrix representing a curvature transformation.

We Also suppose that the noise is independent and identically distributed (i.i.d), following a

Gaussian distribution with a standard deviation σ:

p(E k ) = 1 (2π) N 2 σ N exp{ -E k E k 2σ 2 } (2.2)
Due to the complexity of the problem, we separate it in four steps. This is the same approach used in the RSR algorithm [START_REF] Farsiu | Fast and robust multiframe super resolution[END_REF].

1. Fusing the low-resolution images Y k into a blurred HR version B = HX.

2. Finding the estimation of the HR image X from the blurring one B.

3. Computing the warp matrix F k for each image using curvature registration.

4. Resolving the magnification problem.

We will detail these steps in the next subsections.

The fusing step

The first part of our algorithm is computing the blurred HR version B = HX. Then we use the maximum likelihood estimator (ML) which suggests the choice of B that maximizes the likelihood function (2.6).

B = argmax B {p(Y k /B)} (2.3) = argmin B {-log(p(Y k /B))} = argmin B n k=1 Y k -DF k B 2 L 2
(2.4)

In this case, the matrix (DF k ) (DF k ) is positive definite, the problem is well posed, and has a unique solution:

B = n k=1 ((DF k ) (DF k )) -1 (DF k ) Y k (2.5)
If the problem is ill posed we use TV as a robust regularization, and the minimization problem become:

B = argmin B n k=1 Y k -DF k B 2 L 2 + ∇B 1 (2.6)

Deconvolution and denoising Step

Since X has been known in the presence of white noise, the vector measured Y k is also a Gaussian one. Via Bayes rule, finding the HR image X is equivalent to solve the minimization problem (2.7) using the Maximum a posteriori (MAP).

X M AP = argmax X {p(X/ B)} = argmax X { p( B/X).p(X) p( B) } = argmin X {-log(p( B/X)) -log(p(X))} (2.7)
where p( B/X) represents the likelihood term and p(X) denotes the prior knowledge in the high-resolution image. To solve this problem we need to describe the prior Gibbs function (PDF): p.

The prior Gibbs function

To describe the PDF function we us a second order TV regularisation as we know they robustness to remove noise and misregistration in smoother surfaces.

p(X(x)) = c. exp -α||f (|∇ 2 X(x)|)|| 1 (2.8)
Where f is a linear growth increasing function defined: R -→ R + .

And c is a normalizing constant, guaranteeing that the integral over all x is 1.

|.|: is the euclidean norm of R 4 .

α: is the regularisation parameter verifying 0 < α < 1

The construction of warp matrix F k

We obtain the warp matrix F k for each frame trough a curvature registration algorithm, after a transition from discrete to continuous images using 2-linear interpolation. Let us denote by

Y k (x) the intensity of the kth image of coordinate x ∈ Ω ⊂ R 2
where Ω is defined as the domain of the image. So the expression of the new continuous image Y k (x) [START_REF] Modersitzki | Numerical Methods for Image Registration[END_REF] is given by this formulation :

Y k (x) = k∈{0,1} 2 Y k E(n 1 x 1 ) + k 1 n 1 , E(n 2 x 2 ) + k 2 n 2 × 2 j=1 (-1) k j n j (E(x j n j ) + 1 -k j -x j n j ) (2.9) With: n = (n 1 , n 2 ) is the size of Y k and E(x) is the integer part of x.
In the first step we choose arbitrarily one image Y i from Y k as an image of reference, and we try to find the deformations u k between Y i and the other images, such that:

Y i (x) = Y k (u k (x))
for k = i and ∀x ∈ Ω

An intuitive way to find the deformation u k between the images is to minimize the so-called distance measure D. Since the problem is ill-posed, we have to choose an appropriate regularization R.

The image registration problem is now to find a minimizer u k of the variational problem (2.10) defined by the functional J :

J (u k ) = D(Y i , Y k , u k ) + βR(u k ) for u k ∈ T (2.10)
Where T denotes the set of admissible transformations, and β is the regularisation parameter.

A typical choice for the distance D is the mass-preserving M P measure defined as:

D M P (Y i , Y k , u k ) = Ω (Y k (u k (x))det(∇u k (x)) -Y i (x)) 2 dx (2.11)
In this paper we use a curvature regularisation based on strain tensor. This tensor is generally defined via the displacement v k , and we suppose that

u k (x) = x + v k (x), so u k (x) = I 2 + v k (x). S curv (u k ) = 2 i=1 Ω ( u k i ) 2 dx (2.12)
The registration problem is now well defined in (2.13).

min

u k J curv (u k ) (2.13) With: J curv (u k ) = D M P (Y i , Y k , u k ) + βS curv (u k ) (2.14)
To solve this problem, numerical schemes are required after a discretization of the domain and the objective function J curv . Since we have a curvature registration, we choose a cell-centred grid in the discretization step [START_REF] Modersitzki | Numerical Methods for Image Registration[END_REF]. The objective function depends on the discretization h and defined in (2.15).

J h curv (u h k ) = D h M P (Y h i , Y h k , u h k ) + βS h curv (u h k ) (2.15)
To solve the minimisation problem (2.14) we use the Newton Gauss method (Jorge [START_REF] Nocedal | Numerical Optimization[END_REF], since this problem is nonlinear. We finally find the deformation ûk between the k low resolution images Y k . We can now define easily the matrix of deformation F k easily using ûk .

Resolution of the MAP estimator problem

Since we have defined the prior function p and the operators F k , we rewrite the equation of the MAP estimator:

X M AP = argmin X { x∈Ω HX(x) -B(x) 1 + α||f (|∇ 2 X(x)|)|| 1 } (2.16)
where Ω contains all the pixels on the HR grid X. The norm HX -B 1 is used because it is very robust against outliers [START_REF] Farsiu | Fast and robust multiframe super resolution[END_REF].

To minimise the problem (2.16) we use the classical steepest descent algorithm. We finally have the HR image X as follows.

X n+1 (x) = X n (x) -α {H sing(H X n (x) -B(x)) + α div 2 (g(|∇ 2 X(x)|)∇ 2 X(x))} (2.17)
Where the second order divergence operator div 2 : (R n×m ) 4 -→ R n×m with the adjointness property:

div 2 X.Y = X.∇ 2 Y, ∀X ∈ (R n×m ) 4 , Y ∈ R n×m
And:

g(s) = f (s) s

RESULTS

In this section we evaluate the performance of the proposed algorithm specified at slightly deformed low resolution images. We construct a synthetic LR image to test our algorithm, and compare it with the SR using TV regularization and RSR algorithm with bilateral regularization [START_REF] Farsiu | Fast and robust multiframe super resolution[END_REF]. The peak-signal-to-noise ratio (PSNR) is used to measure the quality of our approach. The SR algorithm is used after an curvature registration, while we use the same data for the RSR resolution and our proposed method. we choose the Lisa (the Simpson's) (1) as the original image with size 256 × 256.

We illustrate in figure (2) four of the N = 30 input low-resolution frames chosen arbitrary and blurring with 5 × 5 Gaussian blur kernel with standard deviation equal to 1.5, and sub-sampling by a factor of 3. In addition we add a noise E k arbitrary in each frame. The parameters chosen for our algorithm are α = 0.1, α = 0.6 and maxiter = 100 iteration for the steepest descent, finally we choose β = 0.1 the curvature registration regularization.

In figure (4d) We show the image obtained using different algorithm compared with our method.

To measure the robustness of the proposed algorithm we use the PSNR in the table (3), we can clearly show the efficiency of our algorithm.

Figure 1 :

 1 Figure 1: The original image of peppers

Figure 4 :

 4 Figure 4: Our result compared with the classical approaches

Acknowledgment

We are grateful to the anonymous referee for the corrections and useful suggestions that have improved this article.


Idriss El Mourabit

Mohammed El Rhabi

Abdelilah Hakim

Orthogonal projection algorithm for first and second order total variation denoising

Keywords: Projection algorithm, Denoising, Total variation, Second order total variation, Partial differential equations

The denoising problem is the process of removing the noise from a degraded image. As we know, the Rodin Osher Fatemi (ROF) denoising model based on total variation is a robust approach for solving the ill-posed problem. To avoid the staircasing effects caused by the first order total variation, the second order one is proposed. In this work, we present an orthogonal projection algorithm for solving the ROF model with first and second order total variation. The efficiency and robustness against noise of the proposed model are illustrated and compared with the classical methods through numerical simulations.

Introduction

The general idea behind variational denoising models is to consider an observed image u 0 as a noiseless version of an image one u. In denoising problems, u is the solution of an ill-posed inverse problem which is presented as follow

u 0 = u + n, (1)
where u 0 is noisy image and n is an additive Gaussian noise.

The ROF denoising model introduced in [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] is one of the most successful variational algorithms, which consists of minimizing the first order total varia-tion. To reconstruct the true image u from u 0 , Rodin Osher Fatemi proposed to minimize the following function inf

u∈BV (Ω) λT V (u) + 1 2 u -u 0 2 L 2 (Ω) (ROF 1 )
where T V (u) represents the total variation of a u ∈ BV (Ω) defined as follows

T V (u) = Ω |Du|dx = sup Ω u div ϕ dx ϕ ∈ C 1 (Ω, R), ϕ ∞ ≤ 1 ,
and BV (Ω) is the space of integrable functions with bounded variations on a bounded domain Ω. This space appears to be a suitable functional space for image analysis, since it contains functions that can be discontinuous across edges. The first term in (ROF 1 ) is a regularization part, controlled by a positive weighting parameter λ, the second one, measures the difference between the original image u and the noised one using the L 2 (Ω) fidelity norm. The existence and uniqueness of the solution of (ROF 1 ) has being studied in [START_REF] Vese | A study in the BV space of a denoising-deblurring variational problem[END_REF][START_REF] Acart | Analysis of bounded variation penalty methods for ill-posed problems[END_REF]. To treat the very noised images, Chambolle proposes in [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] an efficient algorithm to approximate this solution. In fact, the use of the T V in the (ROF 1 ) model favors piecewise constant structures more than smooth structures, this makes staircasing effect. One of the solutions of this problem has been proposed in [START_REF] Bergounioux | A second-order model for image denoising[END_REF],

the authors propose to use the second order total variation in ROF 1 , then the problem writes as follows inf

u∈BV 2 (Ω) λT V 2 (u) + 1 2 u -u 0 2 L 2 (Ω) (ROF 2 )
where BV 2 (Ω) is the space of bounded hessian functions and T V 2 is the second order total variation defined as

T V 2 (u) = Ω | D 2 u | dx (2)
with D 2 u is the second distributional derivative of the function u. In the same reference [START_REF] Bergounioux | A second-order model for image denoising[END_REF] the authors present an other version Chombolle algorithm for resolving the problem (ROF 2 ). On the other hand, the (ROF 1 ) can be resolved with another projection orthogonal algorithm which has been presented in [START_REF] Rhabi | A robust algorithm for convolutive blind source separation in presence of noise[END_REF] in dimension one and is based on a projected gradient method. In this work we generalize this algorithm for minimizing the first and second total variation of an image of dimension two.

The ROF1 model

The theoretical framework

The image restoration problem consists of finding an original image u : Ω → R from a degraded one u 0 and the image u 0 is related to u by the degradation model ( 1) with a Gaussian noise n. According to the maximum likelihood principle u can be approximated by solving the following least-square problem

inf u Ω |u -u 0 | 2 dx (3)
but this problem is ill-posed in the sense of Hadamard and it has been regularized in [START_REF] Tikhonov | Solutions of ill-posed problems[END_REF] by Tickonov and Arsenin by adding a regularization term to [START_REF] Aubert | Mathematical problems in image processing. Partial differential equations and the calculus of variations[END_REF]. The authors proposed to consider the following problem

inf u λ ∇u 2 L 2 (Ω) + 1 2 Ω |u -u 0 | 2 dx. (4)
The associated Euler-Lagrange equation of ( 4) is very strong isotropic and then the edges are not preserved. The Rodin Osher Fatemi (ROF 1 ) based on variation total is an efficient model for removing the noise and preserve the image edges. Solving this problem leads to minimizing the following expression:

inf u∈BV (Ω) λT V (u) + 1 2 Ω |u -u 0 | 2 dx. (ROF 1 )
The following theorem has been proved in [START_REF] Vese | A study in the BV space of a denoising-deblurring variational problem[END_REF][START_REF] Acart | Analysis of bounded variation penalty methods for ill-posed problems[END_REF].

Theorem 1. The problem (ROF 1 ) has an unique solution in BV (Ω) space.

The Euler-Lagrange equation associated to (ROF 1 ) is

λ∂T V (u) + (u -u 0 ) 0 (5)
where ∂T V (u) is the sub-differential of the function T V (u) defined as

w ∈ ∂T V (u) ⇔ T V (v) ≥ T V (u)+ < w, v -u > L 2 (Ω) ∀v
the equation ( 5) can be written

(u 0 -u) ∈ λ∂T V (u)
and according to convex analysis [START_REF] Rockafellar | Convex Analysis[END_REF] we have

u ∈ ∂T V * u 0 -u λ with T V * is the Legendre-Fenchel transform of T V . Equivalently, u 0 λ ∈ (u 0 -u) λ + 1 λ ∂T V * u 0 -u λ we get that w = (u 0 -u) λ is the minimizer of 1 2 w - u 0 λ 2 + u 0 λ T V * (w). (6)
Since T V is a convex and one-homogeneous (T V (λu) = λT V (u) for all u and λ > 0) function and according to convex analysis results [START_REF] Rockafellar | Convex Analysis[END_REF] T V * is the indicator function of a closed convex set K:

T V * (u) =      0 u ∈ K ∞ otherwise
where K is given by

K = u ∈ L 2 (Ω) : u, v L 2 (Ω) ≥ T V (v) ∀v ∈ L 2 (Ω)
from the definition of total variation (2), K can be seen as the closure of the following set

div ξ ξ ∈ C 1 c (Ω, R), ξ ∞ ≤ 1 .
Since the solution w of ( 6) is given by an orthogonal projection of u0 λ on K, and then the solution of (ROF 1 ) can be computed by

u = u 0 -π λK (u 0 ) (7)
where π λK is the projection operator onto K.

Computing the discrete solution with the orthogonal projection algorithm

Form the formula [START_REF] Rhabi | A robust algorithm for convolutive blind source separation in presence of noise[END_REF], resolving total variation problem is equivalent to computing the projection π λK (u 0 ), to this end Chambolle (2004) has been proposed an efficient algorithm for estimating this projection based on a fixed point method. In other way, it could be computed by a projected gradient method, this idea was proposed in [START_REF] Rhabi | A robust algorithm for convolutive blind source separation in presence of noise[END_REF] applied in the case of a 1-D noisy signal, in this section we will generalize this for an image of dimension two. Now we will give some notation, we denote by u i,j , i = 1, ...N, j = 1, ..., M a discrete image and X = R N ×M the set of all discrete images of size N × M and Y = X × X.

The spaces X and Y are both equipped respectively with scalar product < ., . > X and < ., . > Y where

∀u, v ∈ X, < u, v > X = N i=1 M j=1 u i,j v i,j and ∀p = (p 1 , p 2 ), q = (q 1 , q 2 ) ∈ Y, < p, q > Y = N i=1 M j=1 p 1 i,j p 1 i,j + p 2 i,j p 2 i,j
The gradient of an element u written ∇u belongs to Y and could be defined by several manners. One of them consists to set ∇u = ((∇u) 1 , (∇u) 2 ) with

(∇u) 1 i,j =      u i+1,j -u i,j if i < N 0 if i = N (∇u) 2 i,j =      u i,j+1 -u i,j if j < M 0 if j = M (8)
The adjoint operator of -∇ is defined as div :

X 2 → X with, for all p = (p 1 , p 2 ) ∈ X 2 we have ∀w ∈ X, < div p, w >= -< p, ∇w >
When the gradient is given by ( 8) then

(div p) i,j = (div p) 1 i,j + (div p) 2 i,j (9)
where

(div p) 1 i,j =            p 1 i,j -p 1 i-1,j if 1 < i < N p 1 i,j if i = 1 -p 1 i-1,j if i = N (div p) 2 i,j =            p 2 i,j -p 2 i,j-1 if 1 < i < M p 2 i,j if i = 1 -p 2 i,j-1 if i = M
Then the discrete version of T V denoted by J is defined as

J(u) = N i=1 M j=1 |(∇u) i,j |
where |.| is the Euclidean norm.

Then we can write the J as

J(u) = sup {< p, ∇u > Y , p ∈ Y, |p i,j | ≤ 1 i = 1, ..., N and j = 1, ..., M } = sup p∈G < p, u > X
and the set G is defined by

G = {-div p| p ∈ Y, |p i,j | ≤ 1 i = 1,
According to the previous the solution of [START_REF] Tikhonov | Solutions of ill-posed problems[END_REF] is

u = u 0 -π λG (u 0 ) (11)
where π λG (u 0 ) is the orthogonal projection of u 0 on the convex set λG.

Thus, to compute u we are lead to compute the projection operator π λG on the convex set λG, i.e., to solve the problem

min p∈λG F (u) ( 12
)
where

F (p) = p -u 0 2 .
Let ρ > 0, the optimality condition for the problem (12) could be written as follows

p = Π D (p -ρ∇F (p)) (13)
where

D = {p ∈ Y | ∀i, j |p i,j | ≤ 1 }, ∇F (p) = -2λ∇(λ div p -u 0 ) and Π D is
the orthogonal projection on D. It is straightforward to see that

Π D (q) i,j =      q i,j |q i,j | if |q i,j | > 1 q i,j otherwise
the problem can be solved by a classical projection algorithm which give as:

1. p 0 = 0.

2. for n ≥ 0

p n+1 = Π D (p n -ρ∇F (p n ))
where ρ > 0 suitably chosen. Then

p n+1 -p = Π D (p n -ρ∇F (p n )) -Π D (p -ρ∇F (p)) ≤ p n -p + 2ρλ 2 ∇ div (p n -p) ≤ (I + 2ρλ 2 ∇ div)(p n -p) ≤ |I + 2ρλ 2 ∇ div | p n -p
Since the eigenvalues of ∇ div are negative the sequence p n -p is necessary decreasing. Then, it converges when

λ 2 ρ < 1 |µ| = 1 ∇ div 2
where µ is the eigenvalue of ∇div having the largest absolute value.

In order to evaluate the performance of the proposed algorithm in sense of the convergence, we will compare it with the Cambolle algorithm [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] (see numerical implementation section).

The ROF2 model

The computed (ROF 1 ) solution makes some numerical perturbations. In fact, the computed solution turns to be piecewise constant which is called the staircasing effect. In order to resolve this problem it has been proposed to use the second order functional space of bounded variation BV 2 (Ω). This model leads to the minimisation of the following expression:

inf u∈BV 2 (Ω) T V 2 (u) + λ 2 u -u 0 2 L 2 (Ω) (ROF 2 )
where T V 2 is the second total variation define by

T V 2 (u) = sup Ω uH * ϕdx ϕ ∈ C 2 c (Ω, R), ϕ ∞ ≤ 1
and BV 2 (Ω) is space of bonded hessian function defined as

BV 2 (Ω) = u ∈ W 1,1 (Ω), T V 2 (u) < ∞
this space endowed with the norm . BV 2 (Ω) = . W 1,1 (Ω) + T V 2 (.) is a Banach space, for others properties of BV 2 (Ω) see [START_REF] Demengel | Fonctions à hessien borné[END_REF].

In the case of the finite dimension the following theorem has been proved in [START_REF] Bergounioux | A second-order model for image denoising[END_REF].

Theorem 2. The ROF 2 model has an unique solution in BV 2 (Ω) space.

Compute the discrete ROF2 solution

In this section we are going to compute the ROF 2 numerical solution with the proposed orthogonal projection algorithm. We first present some recalls and notations. The we note by X = R N ×M set of all discrete images and Z = X 4 .

Let u be an element of X, the Hessian matrix of u denoted by Hu is identified to a vector of Z and we denote by J 2 the discrete version of the T V 2 defined by

J 2 (u) = N i=1 M j=1 (Hu) i,j R 4 .
Thus, the discretization of the ROF 2 model can be defined as

min u∈X 1 2 |u -u 0 | X + λJ 2 (u) ( 14
)
the associate Euler Lagrange equation of ( 12) is the following

u -u 0 + λ∂J 2 (u) 0 then u -u 0 λ ∈ J 2 (u)
this equivalent [START_REF] Rockafellar | Convex Analysis[END_REF] as

u ∈ J * 2 ( u -u 0 λ )
where J * 2 is the Legendre-Fenchel transform of J 2 . Therefore

u 0 λ ∈ z + 1 λ J * 2 (z) with z = u -u 0 λ solution the following problem min 1 2 z - u 0 λ 2 + u 0 λ J * 2 (z) (15)
thus, the solution of ROF 2 is given as [START_REF] Bergounioux | A second-order model for image denoising[END_REF]

u = u 0 -π λG2 (u 0 )
where π λG2 (u 0 ) orthogonal projection of u 0 on G 2 with

G 2 = {H * q, q ∈ Z, |q i,j R 4 ≤ 1 ∀i, j} .
In [START_REF] Bergounioux | A second-order model for image denoising[END_REF] the authors present a new version of Chambolle algorithm for computing the projection π λG2 (u 0 ), but here we are going to compute it with the proposed orthogonal projection algorithm which can be more efficient in sense of convergence and restoration. Computing the projection π λG2 (u 0 ) is equivalent to resolve the following problem

min z∈G2 z -u 0 2 X or min q∈D2 λH * q -u 0 2 X ( 16
)
where D 2 is the following convex set

D 2 = {q ∈ Z, |q i,j R 4 ≤ 1 ∀i, j} .
We denote G(q) = λH * q -u 0 X then ∇G(q) = 2λH(λH * q -u 0 ).

Therefore, the numerical solution of (15) is given by 1. q 0 = 0.

2. for n ≥ 0

q n+1 = Π D2 (q n -2τ λH(λH * q n -u 0 )) (17)
where τ > 0 is a real number.

Proposition 3. The sequence q n converges to q solution of (15), when n → ∞, if λ and τ satisfy the following condition:

τ λ 2 ≤ 1 HH * Proof. We have q n+1 -q R 4 = Π D2 (q n -2τ λH(λH * q n -u 0 )) -Π D2 (q -2τ λH(λH * q -u 0 )) R 4 ≤ q n -2τ λH(λH * q n -u 0 ) -q + 2τ λH(λH * q -u 0 ) R 4 = q n -q -2τ λ 2 HH * (q n -q) R 4 ≤ q n -q R 4 I -2τ λ 2 HH * this equivalent to q n -q R 4 ≤ q R 4 I -2τ λ 2 HH * n+1 then, q n -q R 4 → 0 if only if I -2τ λ 2 HH * < 1 we conclude that τ λ 2 ≤ 1 HH * .

Numerical implementation

In this section we will present some numerical examples of denoising process to illustrate the difference between the proposed algorithm (13, 17) and both the Chambolle 1 [5] and Chambolle 2 [4]. The comparison is about the speed of convergence.

First of all we consider two examples, the first one is the image of Lena degraded with an additive Gaussian noise of standard deviation σ = 28 and the second one is a color image degraded with the same noise. In the Figures 1 and3 we display these images with corresponding results of each algorithm. As we know that the ROF 1 approach makes staircasing effect, since the resulting image piecewise constant on smooth areas (Fig. 1 and Fig. 3: Chambolle 1 and projection 1), which is not the case with the ROF 2 model the staircasing effect disappears (Fig. 1 and Fig. 3: Chambolle 2 and projection 2). where u is the true image and u * is the denoised image. Table 1 displays the averages of the PSNR and SNR of the four algorithms. We show that the projection 1 and Chambolle 1 have probably the same results, but the proposed algorithm (projection 2) provides the best results in terms of both PSNR and SNR as Chombolle 2. To illustrate the efficiency of the proposed algorithm in sense of convergence, we display in the Figures 2 and4 the two examples.

Conclusion

In this work, we have proposed an algorithm based on projected gradient method for solving the famous Rodin Osher Fatemi model (order one and two).

The numerical implementation is described and experiment results for image denoising have been tested and compared with the classical algorithms. The efficiency of this paper is about the convergence, indeed, the proposed algorithm converges faster than Chambolle algorithms (order one and two), in addition it has a best denoising results using the PSNR and SNR criteria for the second order total variation model compared with the two Chambolle algorithm.

1 2

 12 ..., N and j = 1, ..., M } from the previous notation the ROF 1 is equivalent to solve the following problem min u∈X |u -u 0 | Y + λJ(u)

Figure 1 :

 1 Figure 1: From left to right and from top to bottom : true image, noisy image and denoised images provided by Chambolle 1, projection 1, Chambolle 2 and projection 2.

Figure 2 : 2 and

 22 Figure 2: From left to right and from top to bottom: true image, noisy image and denoised images provided by Chambolle 1, projection 1, Chambolle 2 and projection 2.In addition, we can evaluate the performance of each algorithm, by using the Peak Signal to Noise Ratio (PSNR) and Signal to Noise Ratio (SNR) which is defined by

and 2 .

 2 the comparison of the criterion. The left-curve in the Figure 2 and 4 compare the convergence between the proposed algorithm and Chambolle 1 and the right-curve in the same figures present the comparison between projection 2 and Chambolle 2. We show that the proposed algorithm converge much faster to the solution than Chambolle 1 The parameter λ is selected appropriately for each algorithm, and the gradient descent parameter ρ and τ are respectively ρ = 0.01 and τ = 0.001 for

Figure 3 :

 3 Figure 3: Comparative study of convergence for the image Lena, from left to right: comparison of criterion for Chambolle 1(cham1) and projection (proj1), Chambolle 2(cham2) and projection (proj2).

Figure 4 :

 4 Figure 4: Comparative study of convergence for the image Onion, from left to right: comparison of criterion for Chambolle 1(cham1) and projection (proj1), Chambolle 2(cham2) and projection (proj2).

Table 1 :

 1 Comparison of measures PSNR and SNR.

	measures Chambolle1 Projection1 Chambolle2 Projection2
	PSNR	29.65	29.65	29.61	29.92
	21 cm Lena
	SNR	23.97	23.96	23.95	24.25
	PSNR	29.41	29.38	29.18	29.51
	21 cm Onion
	SNR	22.35	22.32	22.17	22.49

Preprint submitted to Annals of the University of Craiova, Math. & and Comp.Science Series