Bruno Despr És
email: despres@ann.jussieu.fr

NON LINEAR FINITE VOLUME SCHEMES FOR THE HEAT EQUATION IN 1D

Keywords: 2000 Mathematics Subject Classification. 65M06, 65M12 Finite volume schemes, heat equation, non linear correction

published or not. The documents may come     Non linear finite volume schemes for the heat equation in 1D.

Introduction

Finite volume schemes are very convenient for complex applications [START_REF] Eymard | Finite Volume MEthods, Handbook of Numerical Analysis[END_REF][START_REF] Sheng | Monotone Finite volume schemes of non-equilibrium radiation diffusion equations of distorted meshes[END_REF][START_REF] Aavatsmark | Discretization on unstructured grids for inhomogeneous, anisotropic media. Part I: Derivation of the methods[END_REF][START_REF] Buet | Design of asymptotic preserving finite volume schemes for the hyperbolic heat equation on unstructured meshes Numerische Mathematik[END_REF][START_REF] Herbin | Benchmark on discretization schemes for anisotropic diffusion problems on general grids[END_REF][START_REF] Boyer | Finite volume method for 2D linear and nonlinear elliptic problems with discontinuities[END_REF]. A very active field of research is nowadays the development of linear and non linear finite volume methods for the heat equation [START_REF] Boyer | Finite volume method for 2D linear and nonlinear elliptic problems with discontinuities[END_REF][START_REF] Cances | Monotone coercive cell-centered finite volume schemes for anisotropic diffusion equations[END_REF][START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general nonconforming meshes. SUSHI: a scheme using stabilization and hybrid interfaces[END_REF][START_REF] Lipnikov | Local flux mimetic finite difference methods[END_REF][START_REF] Sheng | Monotone Finite volume schemes of non-equilibrium radiation diffusion equations of distorted meshes[END_REF]. A major inspiration is the work of Le Potier [START_REF] Potier | Correction non linéaire et principe du maximum pour la discrétisation d'opérateurs de diffusion avec des schémas volumes finis centrés sur les mailles[END_REF][START_REF] Droniou | Construction and convergence study of local-maximum-principle preserving schemes for elliptic equations[END_REF] who has designed ingenious non linear correction terms to guarantee the maximum principle in any dimension. But to our knowledge only partial convergence results are available in the literature [START_REF] Cances | Monotone coercive cell-centered finite volume schemes for anisotropic diffusion equations[END_REF]. In our mind this is fundamentally related to the lack of an a priori estimate which controls the norm of some discrete derivatives of such schemes. In this paper we will show how to derive such an a priori estimate for some non linear schemes for the discretization of the model non stationary heat equation in dimension one (1.1) ∂ t u -∂ xx u = 0, t > 0, x ∈ R, u(0, x) = u 0 (x), x ∈ R.

The long term ambition of this work is to establish, for the heat equation, the same kind of setting that was fruitful in the recent past for non linear schemes for the advection equation [START_REF] Sweby | High-resolution schemes using flux limiters for hyperbolic conservation-laws[END_REF][START_REF] Roe | Characteristic-based schemes for the Euler equations[END_REF][START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF]. Indeed a beautiful and comprehensive theory based on TVD or TVB schemes is available in 1D for the advection equation. Such TVD schemes are based on a control of the discrete L 1 norm of the first derivative for TVD and TVB schemes [START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Després | Lax theorem and Finite Volume schemes[END_REF][START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF]. In this sense a natural question is to generalize this methodology to 1D non linear schemes for the heat equation, trying to establish a systematic method to obtain bounds on some discrete derivatives, and after that to prove convergence. We will show how to obtain these estimates using an approach proposed in [START_REF] Després | Lax theorem and Finite Volume schemes[END_REF] which is fundamentally Fourier based. In order words we show how to control some non linear terms with Fourier linear techniques: this apparent paradox is the cornerstone of this work. The structure of the main stability estimate is the following: if a certain function of the CFL number ν = ∆t ∆x 2 which is the product of three terms is less than one, that is if

(1.2) Q × α × E(ν) ≤ 1 -ε
where Q > 0 is a natural measure of the size of the non linear correction, α ≥ 1 is a constant characteristic of the Le Potier method (at the end of the analysis one always takes α = 1) and E(ν) is a complicated function of the CFL number, then one obtains a simple control of the discrete derivatives of order k on which the non linear correction is based. The function E(ν) is a series and depends on the type of linear scheme that one considers: typically

(1.3) E(ν) = ν l≤L (1 -4ν) l 2 k + L+1≤l 2l 2l + k l k ν(2l + k) k 2
where L = k 2 1 4ν -1 is a threshold value justified in the core of the paper. This idea comes from [START_REF] Després | Lax theorem and Finite Volume schemes[END_REF][START_REF] Després | Convergence of non-linear finite volume schemes for linear transport[END_REF]. An important part of the paper will be devoted to obtain sharp estimates for E(ν) and related functions and to show that the product (1.2) is indeed smaller than 1.

The net result is this work is the design of a new scheme for the heat equation. This finite scheme scheme (3.19) is based on the Le Potier's trick and has enhanced approximation properties. It is explicit, non linear, preserves total mass, is maximum preserving under standard CFL condition: we are able to prove it converges at order 3 towards smooth solutions; the only restriction of our convergence theory is the CFL number which is for the moment slightly more restrictive than the usual one. Numerical results show this order of convergence is optimal. A second scheme with similar properties is designed in the appendix.

The organization is as follows. We first present the family of high order 1D schemes that we desire to analyze: we call theses schemes Le Potier or modified (Le Potier's) schemes. We explain how to use the Le Potier's trick to modify the schemes and to insure the maximum principle. After we will derive the fundamental and new a priori estimates (1.2-1.3) on which the convergence results of the final section are based. Some numerical results are used to confirm the theoretical analysis.

Basic linear schemes

Let ∆x > 0 is the mesh size of our finite difference or finite volume discretization method. We will consider square integrable numerical profiles v = (v j ) j∈Z such that

(2.1) v =   ∆x j∈Z v 2 j   1 2
.

It is convenient to define the space of square integrable numerical profiles

(2.2) l 2 = v ∈ R Z ; v < ∞
equipped with the natural scalar product

(2.3) (u, v) = ∆x j∈Z u j v j .
We start from the finite volume form of a linear explicit discrete scheme

(2.4) ∆x u j -u j ∆t -f j+ 1 2 -f j-1 2 = 0, ∀j ∈ Z,
where u j is the current discrete solution in cell j and f j+ 1 2 is the explicit numerical flux evaluated between cells j and j + 1. Taking

f j+ 1 2 = u j+1 -u j ∆x ,
one gets the classical three points linear scheme

u j -u j ∆t - u j+1 -2u j + u j-1 ∆x 2 = 0, j ∈ Z.
This scheme is of order 1 in time and 2 in space. The Courant number is ν = ∆t ∆x 2 . The three points scheme is stable in l 2 (actually it is stable in all discrete Lebesgue spaces) under CFL condition

(2.5) 2ν ≤ 1.
This monotone scheme is the fundamental brick for the heat equation with the lowest order of approximation in dimension one. The following examples display enhanced approximation properties for the heat equation. If one desires to establish a parallel with linear schemes for the advection equation, the fundamental brick is the upwind scheme while the Lax-Wendroff scheme which is second order with enhanced approximation properties.

Example 1.

The second scheme that we consider is of order 1 in time and 4 in space. It is based on the observation that

u(x j+1 ) -2u(x j ) + u(x j-1 ) ∆x 2 = ∂ xx u(x j ) + ∆x 2 12 ∂ xxxx u(x j ) + O(∆x 4 )
for smooth functions. Therefore

u(x j+2 ) -2u(x j ) + u(x j-2 ) 4∆x 2 = ∂ xx u(x j ) + ∆x 2 3 ∂ xxxx u(x j ) + O(∆x 4 ).
A linear combination yields

4 3 u(x j+1 ) -2u(x j ) + u(x j-1 ) ∆x 2 - 1 3 u(x j+2 ) -2u(x j ) + u(x j-2 ) 4∆x 2 = ∂ xx u(x j ) + O(∆x 4 ).
A similar trick is used in [START_REF] Cohen | Higher-Order Numerical Methods for Transient Wave Equations[END_REF]. That's why we will consider the scheme

(2.6) u j -u j ∆t - 4 3 u j+1 -2u j + u j-1 ∆x 2 + 1 3 u j+2 -2u j + u j-2 4∆x 2 = 0
which is first order in time and fourth order in space. The corresponding flux is

f j+ 1 2 = 4 3 u j+1 -u j ∆x - u j+2 + u j+1 -u j -u j-1 12∆x .
It is convenient for further developments to use a more compact variational formulation. We define the bilinear form (2.7)

a(u, v) = 4 3 ∆x j (u j+1 -u j ) (v j+1 -v j ) ∆x 2 - 1 3 ∆x j (u j+2 -u j ) (v j+2 -v j ) 4∆x 2 .
With these notations the scheme (2.6) restricted to profiles in u, u ∈ l 2 can be rewritten under the variational form

(2.8) u -u ∆t , v + a(u, v) = 0, ∀v ∈ l 2 .
A classical result is the following.

Proposition 2.1. The symmetric bilinear form (2.7) is non negative.

Proof. One has ∆x j (u j+2u j )

2 ≤ 4∆x j (u j+1 -u j ) 2 from which deduce that a(u, u) ≥ 4∆x 3 j (u j+1 -u j ) 2 ∆x 2 - ∆x 3 j (u j+1 -u j ) 2 ∆x 2 = ∆x j (u j+1 -u j ) 2 ∆x 2 ≥ 0.
An elementary upper bound is

(2.9) a(u, u) ≤ 16 3∆x 2 u 2 .
Lemma 2.2. This scheme is stable in l 2 under the CFL condition

(2.10) 16 3 ν ≤ 1.
Proof. Take the test function v = u in (2.8). One gets

1 2 u 2 - 1 2 u 2 + 1 2 u -u 2 = -∆ta(u, u) = - 1 2 ∆t a(u, u) - 1 2 ∆t a(u, u) + 1 2 ∆t a(u -u, u -u).
Since the bilinear form a is non negative one has that (2.11)

u 2 -u 2 ≤ ∆t a(u -u, u -u) -u -u 2 ≤ 16ν 3 -1 u -u 2 .
The CFL condition guarantees the non positivity of the right hand side. It ends the proof.

Example 2.

It is of course tempting to use the modified equation to design a scheme with enhanced consistency in time. Consider a smooth solution of the heat equation. One has

u(t n+1 -u(t n ) ∆t = ∂ t u(t n ) + ∆t 2 ∂ tt u(t n ) + O(∆t 2 ) = ∂ t u(t n ) + ∆t 2 ∂ xxxx u(t n ) + O(∆t 2 )
from which we deduce that the scheme

u j -u j ∆t - 4 3 u j+1 -2u j + u j-1 ∆x 2 + u j+2 -2u j + u j-2 12∆x 2 - ∆t 2 u j+2 -4u j+1 + 6u j -4u j-1 + u j-2 12∆x 4 = 0
is of order 2 in time and 4 in space. We refer to [START_REF] Yu | The method of differential approximation[END_REF] page 43 where the same scheme is introduced mainly for theoretical purposes. The scheme can also be rewritten under the form

u j -u j ∆t - u j+1 -2u j + u j-1 ∆x 2 + (1 -6ν) u j+2 -4u j+1 + 6u j -4u j-1 + u j-2 ∆x 2 = 0.
It admits a variational reformulation: compute u ∈ l 2 such that (2.12)

u -u ∆t , v + b(u, v) = 0, ∀v ∈ l 2
where the bilinear form b is a correction of a (defined in (2.7))

b(u, v) = a(u, v) - ν 2 a(u, v) with a(u, v) = ∆x j (u j+2 -2u j+1 + u j ) (v j+2 -2v j+1 + v j ) ∆x 4 .
We notice that

0 ≤ a(u, u) ≤ 4∆x j (u j-1 -u j ) ∆x 2 . Therefore b(u, u) ≥ (1 -2ν) ∆x j (u j-1 -u j ) 2 ∆x 2 ≥ 0
under condition (2.5).

Lemma 2.3. The scheme (2.12) is stable in l 2 under the CFL condition (2.10).

Proof. Since b(u, u) ≤ a(u, u), the same proof as the one of lemma 2.2 holds. The final stability condition is the more restrictive one between (2.5) and (2.10). The proof is ended.

Maximum principle and Le Potier's trick

The schemes (2.6)-(2.8) and (2.12) are high order. Unfortunately they do not preserve the maximum principle for all ν. For example the explicit formulation of (2.6) is (3.1)

u j = 1 - 5ν 2 u j + 4ν 3 (u j+1 + u j-1 ) - ν 12 (u j+2 + u j-2 ) .
Since the coefficients of the extreme parts are negative, the scheme cannot be maximum preserving: this is independent of the CFL condition. C. Le Potier has proposed in [START_REF] Potier | Correction non linéaire et principe du maximum pour la discrétisation d'opérateurs de diffusion avec des schémas volumes finis centrés sur les mailles[END_REF][START_REF] Droniou | Construction and convergence study of local-maximum-principle preserving schemes for elliptic equations[END_REF] a way to introduce a non linear modification so as to recover the maximum principle in any dimension. It starts from the scheme rewritten under the form

u j -u j = νDiv j , Div j = ∆x f j+ 1 2 -f j-1 2 .
The neighbors of a cell correspond by definition to non zero coefficient in the explicit stencil: for our example (3.1) it corresponds to

(3.2) V(j) = {j -2, j -1, j + 1, j + 2}
Let us define

(3.3) Σ j = l∈V(j) |u l -u j | .
We note that Σ j ≥ 0 by definition. With this notation, the application of the Le Potier's trick to the discrete heat equation is to consider the non linear scheme (3.4)

u j = u j + νDiv j + αν l∈V(j) a jl (u l -u j )
where α ≥ 0 is a coefficient to be prescribed and the non linear part of the scheme is given by

a jl = a lj = |Div j | Σ j + |Div l | Σ l which is a non linear coefficient.
Remark 3.1. One sees that a jl is a bounded quantity by construction.

If one compares with TVD schemes, this term is very close to some extended slope indicator. By convention a jl can be extended by zero a jl = 0 for l ∈ V(j) ⇐⇒ j ∈ V(l), so that the non linear scheme (3.4) can be rewritten under the more compact form (3.5) u j = u j + νDiv j + αν l∈Z a jl (u lu j )

3.1. Stability. The non linear scheme (3.4) inherits the properties of the stationary Le Potier's schemes [START_REF] Potier | Correction non linéaire et principe du maximum pour la discrétisation d'opérateurs de diffusion avec des schémas volumes finis centrés sur les mailles[END_REF][START_REF] Droniou | Construction and convergence study of local-maximum-principle preserving schemes for elliptic equations[END_REF].

Proposition 3.2. The scheme (3.4) admits the finite volume formulation

∆x u j -u j ∆t -f tot j+ 1 2 -f tot j-1 2 = 0, ∀j ∈ Z,
where the total flux is the sum of the linear flux plus a non linear correction, that is f tot

j+ 1 2 = f j+ 1 2 + f cor j+ 1 2
for all j with

(3.6) f cor j+ 1 2 = α ∆x l≤j+ 1 2 ≤m
a ml (u mu l ).

Remark 3.3. If l and m are two indices sufficiently far one to the other, then a ml = 0. So only a finite number of terms enter in (3.6).

Proof. The linear part of the flux f j+ 1 2 does not yield any difficulty. Concerning the non linear part one has by construction

f cor j+ 1 2 = α ∆x j+1≤m a mj (u m -u j ) + α ∆x l≤j-1 and j+1≤m a ml (u m -u l ).
The first sum is for l = j, the second sum is all other terms. One also has

f cor j-1 2 = α ∆x l≤j-1 2 ≤m a ml (u m -u l ) = α ∆x l≤j-1 a jl (u j -u l ) + α ∆x l≤j-1 and j+1≤m a ml (u m -u l )
where the first sum if for m = j and the second sum is the rest. The difference is therefore

f cor j+ 1 2 -f cor j-1 2 = α ∆x j+1≤m a mj (u m -u j ) - α ∆x l≤j-1 a jl (u j -u l ) = α ∆x l∈Z a jl (u l -u j ).
It ends the proof.

Lemma 3.4. The scheme (3.4) is stable in l 2 under CFL condition

16 3 + 4α max j l a jl ν ≤ 1.
Proof. We define the bilinear form

(3.7) a tot (u, v) = a(u, v) + α ∆x a cor (u, v)
which is the sum of the classical bilinear form a(u, v) that corresponds to the linear part of the scheme and which is non negative, and of the additional bilinear form

a cor (u, v) = - α ∆x j l a jl (u l -u j ) v j = α∆x l j a jl u l -u j ∆x v l -v j ∆x
which corresponds to the non linear coefficients a jl . The equivalent of the stability inequality (2.11) is

u 2 -u 2 ≤ ∆t a + α ∆x a cor (u -u, u -u) -u -u 2 ≤ 16ν 3 + 4αν max j l a jl -1 u -u 2
which proves the claim.

Lemma 3.5. Assume α ≥ 1. Assume the CFL condition

(3.8)   l∈V(j) (α + ε jl ) |Div j | Σ j + α |Div l | Σ l   ν ≤ 1
where ε jl = ±1 is defined below. Then the scheme (3.4) satisfies the maximum principle.

Proof. Using the evident identity

Div j = l∈V(j) Div j |u l -u j | Σ j
one rewrites the explicit scheme (3.5) under the form

u j = u j + ν l∈V(j) Div j |u l -u j | + α |Div j | (u l -u j ) Σ j +αν l∈V(j) |Div l | Σ l (u l -u j )
The key observation is the following

• If Div j (u l -u j ) ≥ 0 then we set ε jl = 1 so that Div j |u l -u j | + α |Div j | (u l -u j ) = (α + ε jl ) |Div j | (u l -u j )
• On the other hand if Div j (u lu j ) < 0 then we set ε jl = -1 so that

Div j |u l -u j | + α |Div j | (u l -u j ) = (α + ε jl ) |Div j | (u l -u j ).
In both cases α + ε jl ≥ 0 since α ≥ 1 by hypothesis. Therefore

u j = u j + ν l∈V(j) (α + ε jl ) |Div j | Σ j + α |Div l | Σ l (u l -u j ) =   1 -ν l∈V(j) (α + ε jl ) |Div j | Σ j + α |Div l | Σ l   u j +ν l∈V(j) (α + ε jl ) |Div j | Σ j + α |Div l | Σ l u l .
It is a convex combination under CFL. It ends the proof.

Let us consider a simple example in order to figure out the practical impact of (3.8) on the stability criterion. The worst guess reduces to

Card V(j) × sup j |Div j | Σ j × sup jl (|α + εjl| + α) ≤ 1.
We consider the example described in (2.6) which corresponds to

Div j = - 1 12 (u j+2 + u j-2 ) + 4 3 (u j+1 + u j-1 ) - 5 2 u j = - 1 12 (u j+2 -u j ) + 4 3 (u j+1 -u j ) + 4 3 (u j-1 -u j ) - 1 12 (u j-2 -u j ).
One has the bound

|Div j | Σ j ≤ 1 12 + 4 3 + 4 3 + 1 12 = 34 12 .
So we obtain

4 × 34 12 × (2α + 1) ν ≤ 1 It is natural to choose the smallest value of α, that is α = 1. One obtains the (sufficient) CFL condition 34ν ≤ 1.
It is quite a stringent condition in terms of time constraint.

Remark 3.6. In order to optimize the time step restriction, that is to be able to take ∆t as large as possible, we will systematically use α = 1 either in the theory or in the numerics. However we keep it most of the analysis to be closer to Le Potier's notations.

A first natural question is to try to diminish the numerical value of the constant in this CFL condition: it will be done with the help of modified schemes developed hereafter.

3.2. Modified schemes. In many cases the initial linear scheme can be decomposed in two parts. A first linear part for which the maximum principle holds and a second linear part which does not naturally respect the maximum principle. We propose to call "modified schemes" such schemes where only the second part is modified by the seminal Le Potier's trick. This procedure which is very natural and has been developed also in [14] is a way to obtain a less severe CFL constraint, still guaranteeing the maximum principle.

3.3.

A first modified scheme. For example the second example (2.6) can be rewritten as (3.9)

u j = u j + νDiv 0 j + νDiv 1 j
where the first part is defined by

Div 0 j = 4 3 (u j+1 -2u j + u j-1 )
and the second part is defined by

Div 1 j = - 1 3 (u j+2 -2u j + u j-2 ) .
The stencil associated to the second part is

V 1 (j) = {j -2, j + 2} .
Since the first part Div 0 j naturally corresponds to a scheme which satisfies the maximum principle, we need to modified only the second part. We obtain the scheme (3.10)

u j = u j + νDiv 0 j + νDiv 1 j + αν l∈V1(j) a 1 jl (u l -u j )
where α ≥ 0 is a coefficient to be prescribed and

a 1 jl = a 1 lj = Div 1 j Σ 1 j + Div 1 l Σ 1 l . The quantity Σ 1 j is defined in accordance by (3.11) Σ 1 j = l∈V1(j)
|u lu j | .

3.4.

A second modified scheme. However the decomposition (3.9) is not the only one. One can consider (3.12)

u j = u j + νDiv 2 j + νDiv 3 j
where the first part is now Div 2 j = u j+1 -2u j + u j-1 and the second part is

Div 3 j = - 1 12 (u j+2 -4u j+1 + 6u j -4u j-1 + u j-2 ) .
The stencil associated to this new second part is the total one

V 3 (j) = V(j) = {j -2, j -1, j + 1, j + 2} .
Since the first part Div 2 j satisfies the maximum principle, we need to modified only the second part. We obtain the scheme (3.13) u j = u j + νDiv 2 j + νDiv 3 j + αν l∈V3(j)

a 3 jl (u l -u j )
where α ≥ 0 is a coefficient to be prescribed and

a 3 jl = a 3 lj = Div 3 j Σ 3 j + Div 3 l Σ 3 l .
The quantity Σ 3 j is defined by

(3.14) Σ 3 j = l∈V 3 (j)
|u lu j | .

3.5.

A third modified scheme. Here we perform a more important modification. We start from the decomposition

u j+2 -4u j+1 + 6u j -4u j-1 + u j-2 = (u j+2 -3u j+1 + 3u j -u j-1 ) -(u j+1 -3u j + 3u j-1 -u j-2 ) We define Σ 4 j = |u j+2 -3u j+1 + 3u j -u j-1 | + |u j+1 -3u j + 3u j-1 -u j-2 | and a 4 jl = a 4 lj = Div 3 j Σ 4 j + Div 3 l Σ 4 l for l ∈ {j + 1, j -1} = V 4 (j) (otherwise a jl = 0
). Let us consider the scheme (3.15) u j = u j + νDiv 2 j + νDiv 3 j +ανa 4 j,j+1 (u j+2 -3u j+1 + 3u ju j-1 )ανa 4 j,j-1 (u j+1 -3u j + 3u j-1u j-2 ). Proposition 3.7. Set α = 1. The scheme (3.15) is of Finite Volume type and satisfies the maximum principle under the standard CFL condition (2.5).

Proof. We focus the second statement of the claim since the first statement is evident from (3.15). Performing the same kind of algebra as before we rewrite the second part of the linear flux as follows

Div 3 j = Div 3 j |u j+2 -3u j+2 + 3u j -u j-1 | Σ 4 j - -Div 3 j |u j+1 -3u j + 3u j-1 -u j-2 | Σ 4 j
Plugging in (3.15) one obtains (3.16)

u j = u j + νDiv 2 j + νw(u j+2 -3u j+1 + 3u j -u j-1 ) -νz(u j+1 -3u j + 3u j-1 -u j-2 ) where w = (α + ε j,j+1 ) Div 3 j Σ 4 j + α Div 3 j+1 Σ 4 j+1 and z = (α -ε j,j-1 ) Div 3 j Σ 4 j + α Div 3 j-1 Σ 4 j-1
.

Here

ε j,j+1 = +1 if Div 3 j (u j+2 -3u j+1 + 3u j -u j-1 ) > 0, -1 if Div 3 j (u j+2 -3u j+1 + 3u j -u j-1 ) ≤ 0,
and

ε j,j-1 = +1 if Div 3 j (u j+1 -3u j + 3u j-1 -u j-2 ) > 0, -1 if Div 3 j (u j+1 -3u j + 3u j-1 -u j-2
) ≤ 0. By construction Div 3 j ≤ 1 12 Σ 4 j for all j, so w, z ≤ 1 4 . The expansion of formula (3.16) is (3.17)

u j = (1 -2ν + 3w + 3z)u j +ν (1 -3w -z) u j+1 + ν (1 -w -3z) u j-1 + νwu j+2 + νzu j-2
from which the result is deduced since all weights are non negative under CFL 2ν ≤ 1: that is 1 -2ν + 3w + 3z ≥ 0, 1 -3wz ≥ 0, 1w -3z, νw ≥ 0 and νz ≥ 0.

3.6. Second order in time. The second order in time scheme (2.12) admits the explicit form (3.18)

u j = u j + νDiv 2 j + (1 -6ν)νDiv 3 j
which is very close to (3.12). So it is easy to modify (3.15) which becomes

(3.19) u j = u j + νDiv 2 j + (1 -6ν)νDiv 3 j +α|1 -6ν|νa 4 j,j+2 (u j+2 -3u j+2 + 3u j -u j-1 ) -α|1 -6ν|νa 4 j,j+1 (u j+1 -3u j + 3u j-1 -u j-2
). Proposition 3.8. Assume 6ν ≤ 1. The scheme (3.19) satisfies the maximum principle under the standard CFL condition (2.5).

Proof. Change w (resp. z) in |1 -6ν|w (resp. |1 -6ν|z) in (3.17). u j = u j + νDiv j + ανg j where Div j is a given linear stencil and

g j = l a jl (u l -u j )
is the non linear correction such that a jl vanishes for |j -l| large enough. We introduce natural notations.

• T is the translation operator, that is

(T u) j = u j+1 , u ∈ l 2 .
• D is the difference operator, that is

D = T -I.
• We note for convenience the operator A ∈ L(l 2 ) such that (Au) j = -Div j , ∀u ∈ l 2 .

A quantity that plays an important role in the analysis developed in this work is the following. Definition 3.9. Let q ∈ N ⋆ be an integer naturally defined by the scheme. Throughout this paper, the continuity constant of g with respect to D q u will be denoted by

Q > 0, that is g ≤ Q D q u .
For example let us consider the scheme (3.15) for which the source g = (g j ) ∈ l 2 is defined by

g j = a 4 j,j+2 (u j+2 -3u j+2 + 3u j -u j-1 ) -a 4 j,j-1 (u j+1 -3u j + 3u j-1 -u j-2 )
Proposition 3.10. The source of the scheme (3.15)

is such that g ≤ 1+ √ 2 12 D 4 u , that is Q = 1+ √ 2 12
and q = 4.

Proof. We first notice that A 3 u ≤ 1 12 D 4 u where we have used natural notations compatible with (3.15). It remains to compute the continuity constant of g with respect to A 3 u. One has the decomposition g j = h j + k j + l j with

h j = (u j+2 -3u j+1 + 3u j -u j-1 ) -(u j+1 -3u j + 3u j-1 -u j-2 ) Σ 4 j Div 3 j , k j = u j+2 -3u j+1 + 3u j -u j-1 Σ 4 j+1 Div 3 j+1 and l j = u j+1 -3u j + 3u j-1 -u j-2 Σ 4 j-1 Div 3 j-1 .
Since |h j | ≤ |Div 3 j |, then h ≤ A 3 u . One the other hand one has that |k j | ≤ α j+1 |Div 3 j+1 | and |l j | ≤ β j-1 |Div 3 j-1 | where the coefficients

α j+1 = |u j+2 -3u j+1 + 3u j -u j-1 | Σ 4 j+1 and β j-1 = |u j+1 -3u j + 3u j-1 -u j-2 | Σ 4 j-1
satisfy 0 ≤ α j , β j and α j + β j = 1 ∀j.

So k + l 2 ≤ ∆x j α j+1 |Div 3 j+1 | + β j-1 |Div 3 j-1 | 2 ≤ ∆x j 2α 2 j+1 |Div 3 j+1 | 2 + ∆x j 2β 2 j-1 |Div 3 j-1 | 2 ≤ 2∆x j α 2 j + β 2 j |Div 3 j | 2 ≤ 2∆x j |Div 3 j | 2 .
Therefore k + l ≤ √ 2 A 3 u from which the result is deduced.

Proposition 3.11. The source of the scheme (3.19) is such that

g ≤ 1 + √ 2 12 |1 -6ν| D 4 u , that is Q = 1+ √ 2
12 |1 -6ν| and q = 4. Proof. Evident from previous proposition and definition (3.19). The |1 -6ν| is because we use the second order in time scheme.

Estimates

With the above notations and definition, any of the previous schemes can be rewritten as u = (I -νA)u + ανg. The Duhamel's formula can be used to express the solution at time step n in function of the initial solution and the correction (4.1)

u n = (I -νA) n u 0 + αν n-1 p=0 (I -νA) n-1-p g p .
It shows that the solution at time step n∆t is the sum of the standard discrete solution plus a contribution due the corrections. Our goal is to show that the global correction is small in some norm. Before we need some a priori estimates on the g p s. These estimates will be obtained using a control of the non linear part of the Duhamel's formula by the linear part. By application of the operator D q to the Duhamel formula (4.1), one gets the identity (4.2)

D q u n = (I -νA) n D q u 0 + αν n-1 p=0 (I -νA) n-1-p D q g p .
Using the continuity of g p with respect to D q u p , it yields the estimate

D q u n ≤ (I -νA) n D q u 0 + Qαν n-1 p=0 (I -νA) n-1-p D q D q u p .
A fundamental property is the following.

Proposition 4.1. Assume l 2 contractivity of the linear part of the scheme. Assume there exists ε > 0 such that

(4.3) Qαν ∞ l=0 (I -νA) l D q ≤ 1 -ε.
Then one has the stability estimate

sup n∈N D q u n ≤ 1 ε D q u 0 ∀u 0 ∈ l 2 .
Remark 4.2. The dissipativity of the scheme insures that lim l→∞ (I -νA) l D q = 0 this property is been used in [START_REF] Després | Lax theorem and Finite Volume schemes[END_REF][START_REF] Després | Convergence of non-linear finite volume schemes for linear transport[END_REF] to study non linear schemes for the transport equation. In this study, the dissipativity of the semi-group (I -νA) l is much stronger since it is a fundamental property of the heat equation. This dissipativity is essential to obtain a control of the non linear part of the Duhamel's formula.

Proof. The l 2 contractivity is in our case equivalent to the stability in l 2 with a stability constant equal to 1, that is (I -νA) l ≤ 1 for all l ∈ N. Let us define Z N = sup n≤N D q u n which satisfies the estimate

Z N ≤ D q u 0 + Qαν N -1 p=0 (I -∆tA) n-1-p D q Z N -1 .
It turns into Z N ≤ D q u 0 + (1ε)Z N which shows that Z N ≤ 1 ε D q u 0 . Since it is true for all N , it shows the claim.

For convenience we define the generic function (4.4)

F (ν) = Qν ∞ l=0 (I -νA) l D q .
Since α = 1 is used to guarantee the maximum property of the schemes, the criterion that we study is ultimately F (ν) < 1. We will now study the function F for the different operators A and different ks. The estimates developed below are used to prove this estimate. The interested reader can jump first to figure 3 where some numerical evaluations of these functions are displayed and then go back to the theory of the next section. The main message of figure 3 is that the condition F (ν) is naturally satisfied. 4.1. Basic operator. Some of the justifications of the bounds used hereafter for general operators A are greatly simplified if one can first prove (4.3) with (4.5)

A 1 = -T + 2I -T -1 .
For a given k, the Fourier symbol of (I -νA 1 )

l D k is λ l (θ) = 1 + ν(e iθ -2 + e -iθ l e iθ -1 k , θ ∈ R.
One has that

|λ l (θ)| = 1 -4ν sin 2 θ 2 l 2 sin θ 2 k .
We assume (4.6) 4ν ≤ 1 which is twice more restrictive than (2.5). It will simplify a lot the analysis and is not a real restriction. We perform the change of variable y = 2 sin θ 2 . Let us define the function

f ν l (y) = (1 -νy 2 ) l y k so that (4.7) (I -νA 1 ) l D k = max 0≤y≤2 f ν l (y). Let us set µ l (ν) = max 0≤y≤2 f ν l (y) l ∈ N.
We finally define

E 1 (ν) = ν ∞ l=0 µ l (ν). Proposition 4.3. Assume k < 2. Then E 1 (ν) = +∞.
Proof. Note that µ l , E 1 , . . . can be studied for real positive k as well, not only integer values. By definition µ l (ν) ≥ f ν l (z) for any z ∈ [0, 2]. So

E 1 (ν) ≥ ν ∞ l=0 (1 -νz 2 ) l z k = z k-2 ∀z ∈ [0, 2].
If k < 2 the right hand side is singular for z = 0. It ends the proof.

Remark 4.4. It will be showed hereafter that k = 2 is also singular. So it must be understood that k > 2 in all the rest of the paper.

Elementary properties of f ν l and µ l (ν) are the following: one has

f ν l+1 (y) ≤ f ν l (y), 0 ≤ y ≤ 2 and d dν f ν l (y) ≤ 0, 0 ≤ y ≤ 2.
Therefore µ l+1 (ν) ≤ µ l (ν) ∀l ∈ N. Next we study the function f l . One has

d dy f ν l (y) = -2lνy(1 -νy 2 ) l-1 y k + k(1 -νy 2 ) l y k-1 = (1 -νy 2 ) l-1 y k-1 -2lνy 2 + k -kνy 2 .
We define

y 2 l = k ν(2l + k)
so that f l increases from y = 0 to y l , then decreases from y l to 1

ν 1 2 . We note that 0 ≤ • • • ≤ y 2 l+1 ≤ y 2 l ≤ • • • ≤ y 2 0 = 1 ν and that 4 ≤ 1
ν due to the CFL condition (4.6). So • Either y l < 2 and

µ l (ν) = f ν l (y l ) = 2l 2l + k l k ν(2l + k) k 2 . • Or y l ≥ 2 and µ l (ν) = f ν l (2) = (1 -4ν) l 2 k . The transition is for L(ν) such that k ν(2(L(ν) + 1) + k) < 4 ≤ k ν(2L(ν) + k) that is (4.8) k 2 1 4ν -1 -1 < L(ν) ≤ k 2 1 4ν -1 which means that L(ν) = k 2 1 4ν -1 . Proposition 4.5. Assuming k > 2, one has the bound E 1 (ν) ≤ h(ν) with (4.9) h(ν) = 2 k-2 1 -(1 -4ν) k 8ν + λ k (ν) 2 k-2 k k -2 + ν2 k .
where

λ k (ν) = 2 k 2 1 4ν -1 + 2 2 k 2 1 4ν -1 + 2 + k [ k 2 ( 1 4ν -1)]+1
.

Remark 4.6. It can be checked this estimate is sharp for small ν. In particular one has that the limit value in 0 + is given by E 1 (ν) = h(0) and (4.10)

h(0) = 2 k 4 + e -k 2 2 k-1 k -2 .
Indeed the first term admits the limit value (1 -4ν) k 8ν tends to e -k 2 when ν vanishes. Concerning the second term one noticed that λ k (ν) =

2(L+1) 2(L+1)+k L+1 with L → ∞ in the regime ν → 0 + . Therefore λ k (0 + ) = e -k
2 by continuity: it yields (4.10). This numerical value can be checked in the numerical figure 2.

Remark 4.7. The same proof shows that E 1 (ν) = +∞ for k = 2, because Z 2 (see below) diverges.

Proof. We obtain the formula for the sum of norms (4.7)

(4.11) ∞ l=0 µ l (ν) = l≤L (1 -4ν) l 2 k =Z1 + L+1≤l 2l 2l + k l k ν(2l + k) k 2 =Z2 .

One has the identity νZ

1 = 2 k-2 1 -(1 -4ν) L+1 . Since L + 1 = k 2 1 4ν -1 + 1 ≤ k 2 1 4ν -1 + 1 ≤ k 8ν - k -2 2 ≤ k 8ν . one has that νZ 1 ≤ 2 k-2 1 -(1 -4ν) k 8ν
.

The second term can be analyzed as a staircase Riemann approximation of a convergent integral. Indeed let us define x l = ν l k , so that νZ 2 can be rewritten as

νZ 2 = k    ν k x l ≤ ν k (L+1) 1 1 + ν 2x l kx l ν 1 2x l + ν k 2    Let γ = k 2 > 0: the function z → 1 1+z γ z is increasing. Since ν 2x l ≤ ν 2( ν k (L+1)) , one obtains that 1 1 + ν 2x l kx l ν ≤ 1 1 + ν 2 ν k (L+1) k ν k (L+1) ν = λ k (ν). So νZ 2 ≤ kλ k (ν)G with G = ν k x l ≤ ν k (L+1) 1 2x l + ν k 2
in which we recognize the staircase approximation of the Riemann's integral of the

decreasing function x → 1 2x+ν k 2 . It is convenient to isolate the first term, that is (4.12) G = ν k 1 2x L+1 + ν k 2 + ν k ν k (L+2)≤x l 1 2x l + ν k 2 =H .
The first term can be bounded using

2x L+1 + ν = 2 ν k k 2 1 4ν -1 + 1 + ν ≥ 2 ν k k 2 1 4ν -1 + 1 + ν ≥ 1 4 which yields ν k 1 2x L+1 + ν k 2 ≤ ν2 k k .
The rest can be bounded using ν k

1 2x l +ν k 2 ≤ x l x l -ν k dx (2x+ν) k 2 . Therefore H ≤ ∞ ν k (L+1) dx (2x+ν) k 2 . Since ν k (L + 1) = ν k k 2 1 4ν -1 + 1 ≥ ν k k 2 1 4ν -1 = 1 8 - ν 2 one gets that H ≤ ∞ 1 8 -ν 2 dx (2x + ν) k 2 = 4 k 2 -1 k -2 = 2 k-2 k -2 . That is νZ 2 ≤ kλ k (ν) ν 2 k k + 2 k-1 k-2 .
It ends the proof.

4.2.

Fourth order operator. Next we consider the operator (4.13)

A 2 = -T + 2I -T -1 + 1 12 T 2 -4T + 6I -4T -1 + T -2
which corresponds to the scheme (2.8) or (3.1). For a given k, the Fourier symbol of (I -νA 2 )

l D k is λ l (θ) = 1 + ν(e iθ -2 + e -iθ ) - 1 12 ν e 2iθ -4e iθ + 6 -4e -iθ + e -2iθ l × e iθ -1 k .
One has that

|λ l (θ)| = 1 -4ν sin 2 θ 2 - 4 3 ν sin 4 θ 2 l 2 sin θ 2 k .
The method of analysis is very similar to the previous one. We assume for simplicity of the analysis that

(4.14) 16 3 ν ≤ 1.
Let us define the function

g ν l (y) = (1 -νy 2 - 1 12 νy 4 ) l y k , y = 2| sin θ 2 |, so that (4.15) (I -νA 2 ) l D k = max 0≤y≤2 g l (y). Let us set σ l (ν) = max 0≤y≤2 g ν l (y)
together with

E 2 (ν) = ν ∞ l=0 σ l (ν).
Proposition 4.8. Assume the CFL condition (4.14). Then

• σ l (ν) ≤ µ l (ν) for all l: it yields E 2 ≤ E 1 and E 2 ∈ L ∞ 0, 3 16 .
• σ l+1 (ν) ≤ σ l (ν) for all l.

Proof. Evident since g ν l (y) ≤ f ν l (y) for 16 3 ν ≤ 1.

The extremal point of g l is z l such that g ′ l (z l ) = 0. That is

(2l + k)z 2 l + k 12 + l 3 z 4 l = k ν .
One sees of course that z l ≤ y l . The solution is

z 2 l = -(2l + k)z + (2l + k) 2 + 4 k 12 + l 3 k ν 2 k 12 + l 3 ,
or also in the conjugate form

z 2 l = 2 k ν (2l + k)z + (2l + k) 2 + 4 k 12 + l 3 k ν .
The value of σ l (ν) is as follows

• Either z l < 2 so σ l (ν) = 1 -νz 2 l - 4 3 νz 4 l l z k l .
• Or n l ≥ 2 and

σ l (ν) = g ν l (2) = 1 - 16 3 ν l 2 k .
The transition is the largest M such that

(2l + k)z 2 M + k 12 + l 3 z 4 M < k ν . Therefore (4.16) E 2 (ν) = l≤M 1 - 16 3 ν l 2 k + M +1≤l 1 -νz 2 l - 4 3 νz 4 l l z k l .
4.3. Full fourth order operator. Finally we consider the Fourier symbol of (I -νA 3 ) l D k of the operator (4.17)

A 3 = -T + 2I -T -1 + 1 -6ν 2 T 2 -4T + 6I -4T -1 + T -2
which corresponds to the second order in time and fourth order in space linear scheme. The scheme is full fourth order in the sense that second order in time corresponds to fourth order in space under CFL condition 6ν ≤ 1. We define the function

E 3 (ν) = ν|1 -6ν| ∞ l=0 (I -νA 3 ) l D k .
Using the same method as before, it can be checked that

E 3 (ν) = ν|1 -6ν| l≤N 1 -4 + 4 3 (1 -6ν) ν l 2 k +ν|1 -6ν| N +1≤l 1 -νw 2 l - 4 3 ν(1 -6ν)w 4 l l w k l
where the extremal point of the function

y → h ν l (y) = 1 -νy 2 -(1-6ν) 12 νy 4 l y k is w 2 l = 2 k ν (2l + k)z + (2l + k) 2 + 4 k 12 + l 3 k ν (1 -6ν)
.

The transition is the largest N such that

(2l + k)w 2 N + k 12 + l 3 (1 -6ν)w 4 N < k ν . Since h ν l (y) ≤ f ν l (y) one has that (4.18) E 3 (ν) ≤ |1 -6ν|E 1 (ν).
4.4. Application to the scheme (3.15). In the following we apply the various inequalities to the case k = 4 and

Q = 1+ √ 2
12 |1 -6ν| which allowed a complete analysis of the scheme (3.19) (see also proposition 3.11).

Our first task is to show that F 3 (ν) = QE 3 (ν) is such that F 3 < 1 so that the stability bound D 4 u p ≤ C D 4 u 0 ∀p holds from proposition (4.1). A first simple and convenient estimate that can be derived from (4.9) is the following. Proposition 4.9. Assume k = 4 and ν ≤ 16. Then the function h defined in (4.9) is such that h ≤ p where

p(ν) = 4 1 -(1 -4ν) 1 2ν + (1 -4ν) 1 2ν -2 (8 + 16ν).
Proof. Indeed our hypotheses imply that k 8ν = 1 2ν on the one hand. So

2 k-2 1 -(1 -4ν) k 8ν = 4 1 -(1 -4ν) 1 2ν
.

One the other hand one has k

2 1 4ν -1 ≥ k 2 1 4ν -1 -1: since the function λ k is monotone decreasing, it yields λ k (ν) ≤ k 1 4ν -1 k 1 4ν -1 + k k 2 ( 1 4ν -1)
.

With k = 4 it turns into λ 4 (ν) ≤ (1 -4ν) 1 2ν -2 from the result is deduced.

Proposition 4.10. One has p ≤ q with

q(ν) = 4 + 4e -2 (1 + 24ν)(1 + 12ν), ν ≤ 1 6 .
Proof. One first has that On the interval [0, 1 6 ], one has that e 8ν ≤ 1+24ν which is easy to show by convexity.

(1 -4ν) 1 2ν -2 = e 2 y (1-y) ln(1-y) , y = 4ν. Since (1 -y) ln(1 -y) = (1 -y) -y - 1 2 y 2 - 1 3 y 3 + . . . = -y + 1 - 1 2 y 2 + 1 2 - 2 3 y 3 + • • • ≤ -y + y 2 , y = 4ν ≤ 1, then (1 -4ν)
Proposition 4.11. On the interval [0, 1 6 ], the function (1 -6ν)q(ν) is monotone decreasing from 4 + 4e -2 until 0.

Proof. For such a simple function, we consider that a plot (with gnuplot) is enough. We refer to figure 1. Lemma 4.12. Assume 6ν ≤ 1. The function F 3 associated to the scheme (3.15) is such that

F 3 (ν) ≤ 1 + √ 2 12
4 + 4e -2 ≈ 0.91 . . .

As a consequence the fourth order discrete derivative is bounded uniformly with respect to n and v: there exists C > 0 such that 12 (1 -6ν)q(ν) in the range ν ∈ [0, 1 6 ]. The function is monotone decreasing.

D 4 u p ≤ C D 4 u 0 p ∈ N.

Proof. By construction

F 3 (ν) = QE 3 (ν) ≤ Q(1 -6ν)E 1 (ν) ≤ Q(1 -6ν)q(ν) ≤ Qq(0) = 1 + √ 2 12 4 + 4e -2 .
Numerical application show that

1 + √ 2 12 4 + 4e -2 ≈ 0.913647 • • • < 1.
Therefore F 3 < 1 over the range ν ∈ 0, 1 6 . See also the plot in figure 3.

Convergence

Now we consider any scheme (4.1) such that an estimate on the norm of g p can be obtained. Let us go back to the Duhamel formula (4.1) and consider the reminder

R = ν n-1 p=0 (I -νA) n-1-p g p , A = 1 1 , A 2 or A 3 .
Due to the fact that corrections compatible with the Le Potier's technique may be rewritten under a finite volume form as stated in proposition 3.2, there exists

s p ∈ l 2 such that g p = Ds p , s p ≤ C 2 D 4 u . Therefore one has R ≤ C 2 ν n-1 p=0 (I -νA) n-1-p D D 4 u 0 ≤ C 2       ν n-1 p=0 (I -νA 1 ) n-1-p D =Fn       D 4 u 0 .
Proposition 5.1. Let T > 0. Using the previous assumptions, there exists

C 3 > 0 such that F n ≤ C 3 ∆x , n∆t ≤ T.
Proof. The term between parenthesis is very similar to the series analyzed previously, but here k = 1 so that the infinite series cannot converge, see proposition 4.3. In order to bound this term we consider the truncated series

F n ≤ ν T ∆t l=0 µ l (ν).
Using (4.11) one gets

F n ≤ ν l≤L (1 -4ν) l 2 =Z1 +ν T ∆t L+1≤n 2l 2l + 1 l 1 ν(2l + 1) 1 2 =Z2 , L = 1 2 1 4ν -1 .
As before the first term is bounded since an elementary summation yields νZ 1 ≤ ν 2 4ν = 1 2 . The second term is bounded by exactly the same method as in the proof of (4.9). The only but essential difference is the sum H which becomes now

H ′ = ν ν(L+2)≤x l ≤ν T ∆t 1 2x l + ν 1 2
We obtain Remark 5.2. Inequality (5.1) can be obtained directly from the fundamental inequality derived in [START_REF] Després | Lax theorem and Finite Volume schemes[END_REF][START_REF] Després | Convergence of non-linear finite volume schemes for linear transport[END_REF].

H ′ ≤ ν T ∆t 1 8 -ν 2 dx (2x+ν)
Let Π ∆x be the point-wise projector of a smooth function onto the mesh grid. Take the initial data as u 0 = Π ∆x u(0). Theorem 5.3. Consider the scheme (3.19). Assume u 0 ∈ H 4 (R). Assume the CFL condition 6ν ≤ 1. Let T > 0. Then there exists a constant C > 0 such that

(5.2) u n -Π ∆x u(n∆t) ≤ C∆x 3 |u 0 | H 4 (R) , ∀n, n∆t ≤ T.
Proof. It is an easy matter to show that the regularity assumption implies

D 4 u 0 ≤ c ′′ |u 0 | H 4 (R) ∆x 4 . Therefore R ≤ K|u 0 | H 4 (R) ∆x 4-1 = K|u 0 | H 4 (R) ∆x 3 .
On the other hand the linear part of the scheme is fourth order in space and second order in time. Therefore

(I -νA) n u 0 -Π ∆x u(n∆t) ≤ C∆x 4 |u 0 | H 4 (R) .
Since u n = (I -νA) n u 0 + R, the triangular inequality shows the result.

This theorem can be adapted to take into account the others Le Potier or modified schemes considered in this work. In particular the curve F 2 in figure 3 is under threshold 1. So the schemes (3.15) also converges at order 3.

Remark 5.4. A case of interest is the scheme (3.13) for which one can prove that Q = 1 4 and q = 4. A proof is given in the appendix. Considering the curves F 2 and F 3 of figure 3 that are normalized for Q = 1+ √ 2

12 , one can apply a correction factor λ = 4 1+ √ 2 ≈ 1.24 . . . . One gets that F 3 = λF 3 < 1 for all ν, and that F 2 = λF 2 < 1 for ν < .15 approximatively. On this range it yields a control of the fourth order discrete derivative, and therefore a proof of convergence at order three if one considers the use of the non linear correction (3.13) with the 1 -6ν term.

Remark 5.5. The curve F 3 is smaller than the limit value in a range [0, C[ where the constant C > 1 6 can be identified with numerical experiments. We infer that the scheme converges at order 3 in this larger range also.

Remark 5.6. The main open problem is the adaptation of such theorem of convergence in dimension two and greater. This is fully open problem.

Numerical tests

We perform simple numerical tests to assess the properties of the numerical schemes developed in this work. l=0 σ l (ν). We observe that these quantities are pretty constant, less than 6 for ν ≤ 1 4 . The computed value E 1 (0) is very close to the exact value ≈ 4 + 4e -2 ≈ 4.5413 . . . . Next in figure 3 we plot

F i = 1 + √ 2 12 E i , i = 1, 2, 3.
We observe that F 2 < 1 and F 3 < 1 for ν ≤ 1 4 . On the other hand F 1 < 1 is true only for approximatively ν ≤ .18. By inspection of the graphics, it is clear that the function F 3 is decreasing and bounded as stated in lemma 4.12.

(ν) = 1+ √ 2 12 E 1 (ν) is for theoretical understanding. F 2 (ν) = 1+ √ 2
12 E 2 (ν) corresponds to the scheme (3.15).

F 3 (ν) = 1+ √ 2
12 E 3 (ν) corresponds to the scheme (3.19). What is important is to be under the threshold 1 (in bold). 6.2. Stability test. We consider the numerical solution of the heat equation on a 10 cells mesh. The initial data is a discrete Dirac mass. Such initial profiles are very convenient to illustrate the maximum principle. We observe in figure 4 that the two fourth order in space linear schemes do not preserve the maximum principle. The three points scheme and the new third order non linear scheme preserve it.

Discrete solution at first iteration Zoom

Figure 4. Numerical solution calculated by the four different schemes: one iteration. One sees that the three points schemes and the modified scheme satisfy the maximum principle. The two fourth order in space schemes do not.

6.3. Accuracy test. We solve the heat equation on the interval [0, 1] with periodic boundary conditions. The initial data is u 0 (x) = cos(2πx) so that the exact solution is u(t) = e -4π 2 t u 0 .

We measure in table 1 and figure 5 the relative error in l 2 norm at time t = 0.1 in function of the number of cells. As predicted by the theory, the modified scheme based on a the second order in time and fourth order in space scheme converges at order 3. 1. Error in function of the mesh size for the four different schemes: S 1 = second order in space and first order in time; S 2 = fourth order in space and first order in time; S 3 = fourth order in space and second order in time; S 4 = S 3 +non linear correction. Appendix A. Constant Q for the scheme (3.13)

cells

The control of the fourth order discrete difference used in the scheme (3.13) can be analyzed as in propositions 3.10 and 3.11. This scheme is also a modified (Le Potier's) scheme, but closer to the initial Le Potier's scheme than (3.15) or (3.19).

For convenience we start from the fourth order in space and second order in time linear scheme, plus the non linear correction (A.1) u j = u j + νDiv 2 j + (1 -6ν)νDiv 3 j + ν|1 -6ν| l∈V3(j)

a 3 jl (u l -u j )
where all terms are defined in (3.13): see also below. The main point is to determine the continuity constant of the non linear correction g = (g j ) where is sufficiently small to be sure that the control of the fourth order difference holds. It yields the convergence at order 3 of this scheme. Notice however that the time step to achieve the maximum principle is a priori more stringent than for (3.19).

g j = l∈V3(j)

3. 7 .

 7 General formulation. Motivated by the previous examples, we will study the following family of modified schemes(3.20)

1 2ν - 2 ≤ 1 2ν -2 8 + 16ν - 4 ( 1 -4ν) 2 = 4 + ( 1 -4ν) 1 2ν -2 4 + 1 2ν - 2

 1218412411412 e -2+2y = e -2+8ν . Therefore p(ν) = 4 + (1 -4ν) 48ν -64ν 2 ≤ 4 + (1 -4ν) (4 + 48ν) ≤ 4 + 4e -2 e 8ν (1 + 12ν)

Remark 4 . 13 .Figure 1 .

 4131 Figure 1. Plot of 1+ √ 2

1 2 (

 12 = I). For small ∆t the integral diverges like (5.1) I ≈ C ν T ∆t ≤ C ′ (T ) ∆x thanks to the CFL condition. It ends the proof.

Figure 2 . 6 . 1 . 10 6 l=0 6 l=0

 2611066 Figure 2. The constant E 1,2,3 for three different set of parameters.

Figure 3 .

 3 Figure 3. The functions F 1,2,3 for three different set of parameters. F 1 (ν) = 1+

Figure 5 .

 5 Figure 5. Error of the different schemes with respect to the number of cells[START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general nonconforming meshes. SUSHI: a scheme using stabilization and hybrid interfaces[END_REF][START_REF] Sweby | High-resolution schemes using flux limiters for hyperbolic conservation-laws[END_REF] 40, 80, 160). Log-scaled. Same as table 1.

3 j 3 j 3 j 2 = d j- 1 3 j- 2 . 3 j= 1 2 + |b j | 2 + |c j | 2 + |d j | 2 3 j 2 ≤ 4 Div 3 2 . 4 ) 4 D 4 u

 33321323122223242444 = l∈V 3 (j) |u lu j | , V 3 (j) = {j + 2, j + 1, j -1, j -2}. One has that g j = h j + k j+2 + l j+1 + m j-1 + n j-2 where the first term h = (h j ) ish j = l∈V 3 (j) (u lu j ) Σ Div =⇒ |h j | ≤ Div 3 j =⇒ h ≤ Div 3and the other terms are             k j+2 = a j+2 Div 3 j+2 , a j+2 = Div 3 j-2 , d j-2 = uj-2-uj ΣOne has the relation|a j | + |b j | + |c j | + |d j | = |u ju j-2 | + |u ju j-1 | + |u ju j+1 | + |u ju j+2 | Σ ∀j. Since k + l + m + n 2 ≤ ∆x j |a j+2 | Div 3 j+2 + |b j+1 | Div 3 j+1 + |c j-1 | Div 3 j-1 + |d j-2 | Div 3 Div Therefore g ≤ h + k + l + m + n ≤ (1 + √ Div 3 = 3 Div 3 ≤1by definition of Div 3 , see(3.12). We now refer remark 5.4 where it is proved that this numerical value Q = |1-6ν|

4

 4

Acknowledgements The author deeply wishes to thank Christophe Le Potier and Clément Cances for valuable comments.