Approximation of nonlinear parabolic equations using a family of conformal and non-conformal schemes

We consider a family of space discretisations for the approximation of nonlinear parabolic equations, 
such as the regularised mean curvature flow level set equation, using semi-implicit or fully implicit time schemes. 
The approximate solution provided by such a scheme is shown to converge 
thanks to compactness and monotony arguments. Numerical examples show the 
accuracy of the method.


Introduction
Nonlinear parabolic equations are involved in different physical or engineering frameworks. For example, the porous medium equation u t − ∆u m = 0, the Stefan problem u t − ∆ϕ(u) = 0 arise in the framework of fluid flows within porous media. Important improvements in the approximation of their solutions have been obtained, using finite volume methods. Indeed, such methods are well suited to the conservative form of these equations. More surprising is the success of finite volume methods for the approximation of some nonlinear problems, under the more general form u t − F (u, ∇u, D 2 u) = 0. For example, in [15], a few algorithms are proposed for the approximation of motion by mean curvature equation, including finite volume methods, whereas the equation, namely u t − |∇u|div (∇u/|∇u|) = 0, is not in the divergence form. In such cases, finite difference methods have more intensively been used. The mathematical framework which is under consideration for the analysis of the convergence of these finite difference schemes relies on the notion of viscosity solution and monotonous scheme. Such a monotonous behaviour does not seem straightforward in the framework of finite volume schemes. Indeed, in a recent paper [7], we study a finite volume method for the approximation of the motion by mean curvature equation in a regularised sense. The principles, used in [7] for the mathematical analysis of the convergence of the finite volume scheme, completely differ from that of the viscosity solutions [5,3], and do not allow for handling the case of the non-regularised motion by mean curvature equation (nevertheless, this case is handled in some numerical examples provided at the end of this paper). A regularised sense, as detailed below, must be used for the proof of convergence of the method.
In the present paper, our aim is to propose a more general framework of approximation methods for some nonlinear parabolic equations in non-conservative form.
In the spirit of [7], where we prove the convergence of a finite volume scheme for the approximation of a weak solution of (1)-(2)-(3) in the sense of Definition 1.1, we develop in this paper a series of new features: 1. We consider a more general framework for the space discretisations, including conformal and nonconformal finite element methods and finite volume methods inspired by multipoint flux approximation [1].
2. In [7], the discrete norms involved in the scheme as arguments of functions µ and ν do not correspond to the exact L 2 norm of an approximate gradient (this imposes to separately prove the strong convergence of this approximate norm, and of an approximate gradient), whereas we consider in this paper a family of schemes such that exact norms of the discrete gradients are used, which allows to directly prove the strong convergence of the gradient from the convergence of its norm. Hence we can more easily consider in this paper the framework of a function ν(u, ∇u) instead of ν(u, |∇u|), since we prove the strong convergence of the discrete gradient used in the discretisation of ν(u, ∇u).
3. We present numerical schemes which can resume to 9-point stencil finite volume scheme (see Section 3.2, where the local elimination of interface unknowns is possible).
4. The proof that the discrete gradient and its norm are strongly convergent relies on Hypothesis (H5) instead of Leray-Lions method [11] (see (39)).
The main result of this paper, i.e. the strong convergence of the discrete schemes to a solution of (5), is proved thanks to the following property. Let F be the function defined by Then, for any sufficiently regular function u, it holds Therefore, assuming that this function u is solution of (1) with f = 0 for the sake of simplicity, we get, by taking v = u t in (5), that ∇u ∈ C 0 ([0, T ]; L 2 (Ω)) and The discrete equivalent of this property is shown in Lemma 4.1 for the fully-implicit scheme (using that x → xµ(x) is strictly increasing), and in Lemma 4.3 for the semi-implicit scheme (using that µ is decreasing). Note that the hypothesis that x → xµ(x) is strictly increasing is used in both schemes for the proof of the strong convergence of the discrete approximate of the gradient. Unfortunately, although it is possible to extend some of these properties to the case µ(x) = 1/x, the convergence study provided in this paper does not hold in this framework.
Remark 1.2 Note that, thanks to the convergence result proved in this paper, we also prove the existence of a weak solution u of (1)-(2)- (3) in the sense of Definition 1.1, which satisfies, for all T > 0: 1. u ∈ L 2 (0, T ; H 1 0 (Ω)) and u t ∈ L 2 (Ω×]0, T [) (hence u ∈ C 0 (0, T ; L 2 (Ω))), This paper is organised as follows. In Section 2, we present a family of discretisation tools, examples of which (case of rectangular or simplicial meshes) are given in Section 3. Then in Section 4, we show some estimates that are used on one hand in the proof of the existence of at least one solution to the fully implicit scheme, and of the existence and uniqueness of the solution to the semi-implicit scheme, on the other hand, in the convergence proof provided in Section 5. Finally, numerical results are given in Section 6.

The family of discrete schemes
We now introduce the tools used for prescribing the space discretisation.
We can notice that the operators used in Definition 2.1 are quite general, and provided by a large variety of discretisation schemes.
Then ∇ D is the natural gradient in H 1 0 (Ω) and Π D is equal to identity. This example will not be further considered in this paper, since we prefer focusing on non-conformal methods including finite volume ones.
Let us now turn to space-time discretisations. We now define two numerical schemes. The fully implicit scheme is defined by and The semi-implicit scheme is defined by (9), and by where u is the function defined by u n−1 (x) for a.e. (x, t) ∈ Ω×](n − 1)τ, nτ [. In order to be able to prove convergence properties, we now give the framework which is considered for a sequence of space discretisations. and 2. the following consistency property holds 3. the following compactness property holds: for all sequence (u m ) m∈N with u m ∈ H Dm such that there exists C > 0 with u m Dm ≤ C for all m ∈ N, then there exists u in L 2 (Ω) such that, up to a sub-sequence, Π Dm u m converges to u in L 2 (Ω), 4. for all sequence (u m ) m∈N with u m ∈ H Dm such that there exists C > 0 with u m Dm ≤ C for all m ∈ N, and such that there exists u ∈ L 2 (Ω) such that Π Dm u m converges to u in L 2 (Ω), then ∇ Dm u m converges to ∇u for the weak topology of L 2 (R) d , prolonging by 0 all functions outside Ω.

Remark 2.2
It results from the above definition that, if a sequence (u m ) m∈N with u m ∈ H Dm is such that there exists C > 0 with u m Dm ≤ C for all m ∈ N, and that there exists u in L 2 (Ω) such that u m converges to u in L 2 (Ω), then u ∈ H 1 0 (Ω). The next section is devoted to the presentation of precise examples of space discretisations, and to the detailed expression of Schemes (10) and (11) in these cases.

Examples of non-conformal space discretisations
Since the main applications which are considered are devoted to image processing, we first focus on non conformal rectangular finite elements on rectangular domains, and then on non conformal simplicial finite elements on polygonal domains. All these non conformal finite element methods can also be seen as finite volume methods.

A first scheme on rectangular domains
We consider the particular case where A space discretisation in the sense of Definition 2.1 is defined by the following way.

We denote by
the set of the control volumes. The elements of M are denoted p, q, . . .. We denote by x p the centre of p. For any p ∈ M, let ∂p = p \ p be the boundary of p; let |p| > 0 denote the measure of p and let h p denote the diameter of p and h D denote the maximum value of (h p ) p∈M .
3. We denote by E the set of all the faces of p ∈ M, and for all σ ∈ E, we denote by |σ| its (d − 1)dimensional measure. For any σ ∈ E, we define the set M σ = {p ∈ M, σ ∈ E p } (which has therefore one or two elements), we denote by E p the set of the faces of p ∈ M (it has 2 d elements) and by x σ the centre of σ. We then denote by d pσ = |x σ − x p | the orthogonal distance between x p and σ ∈ E p and by n p,σ the normal vector to σ, outward to p.

4.
We denote by V p the set of all the vertexes of p ∈ M (it has 2 d elements), by V the union of all V p , p ∈ M. For y ∈ V, we denote by K p,y the rectangle whose faces are parallel to those of p, and whose the set of vertexes contains x p and y. We denote by V σ the set of all vertexes of σ ∈ E (it has 2 d−1 elements), and by E p,y the set of all σ ∈ E p such that y ∈ V σ (it has d elements).

We define the set H
6. We denote, for all u ∈ H D , by Π D u ∈ L 2 (Ω) the function defined by the constant value u p a.e. in p ∈ M.

We denote, for all
8. For u ∈ H D , p ∈ M and y ∈ V p , we denote by and by ∇ D u the function defined a.e. on Ω by ∇ p,y u on K p,y .
It is then possible to show, using the results of [6] that a sequence of space discretisations defined as above is an admissible sequence of space discretisations in the sense of Definition 2.3, letting h Dm tending to 0 whereas min p∈M min σ∈Ep dpσ hp remains uniformly bounded. The proof of this result is a consequence of the fact that the discrete norm is related to the one which is used in the finite volume setting (see [6]). Let us write the schemes (10) and (11) in this case. We first choose for test function v ∈ H D,τ , the function such that v n p = 1 for a given p ∈ M and n = 1, . . . , N T , and all other components equal to 0. We get where we set with m = n for (10) and m = n − 1 for (11), and We then choose for test function v ∈ H D,τ , the function such that v n σ = 1 for a given interior face σ common to both control volumes p, q ∈ M and n = 1, . . . , N T , and all other components equal to 0. We obtain µ n The above expression allows, in the case of Scheme (11), for eliminating u n σ with respect to u n p and u n q . It is then easy to derive an L ∞ estimate in this case, which resumes to L ∞ stability if f = 0.

A second scheme on rectangular domains
We again consider the particular case where 5. We define the set H D of all u = ((u p ) p∈M , (u σ,y ) σ∈E,y∈Vσ ), with u σ,y = 0 for σ ⊂ ∂Ω, y ∈ V σ and n = 1, . . . , N T .
6. We denote, for all u ∈ H D , by Π D u ∈ L 2 (Ω) the function defined by the constant value u p a.e. in p ∈ M.

We denote, for all
and by ∇ D u the function defined a.e. on Ω by ∇ p,y u on K p,y . Let us write the schemes (10) and (11) in this case. For a given p ∈ M and n = 1, . . . , N T , we get where we set with m = n for (10) and m = n − 1 for (11), and For a given interior σ common to p, q ∈ M, y ∈ V σ and n = 1, . . . , N T , we have Again, the above expression allows to eliminate u n σ,y in the case of Scheme (11), and an L ∞ estimate is derived.

A scheme applying on simplicial meshes
This scheme has a few common points with the scheme presented in Section 3.2, although we now consider that Ω be an open bounded polyhedron in R d . A space discretisation in the sense of Definition 2.1 is now defined by the following method. such that Ω = p∈M p. The elements of M are denoted p, q, . . .. We denote by x p the centre of gravity of p. For any p ∈ M, let ∂p = p \ p be the boundary of p; let |p| > 0 denote the measure of p and let h p denote the diameter of p and h D denote the maximum value of (h p ) p∈M .
2. We denote by E the set of all the faces of p ∈ M, and for all σ ∈ E, we denote by |σ| its (d − 1)dimensional measure. For any σ ∈ E, we denote by M σ = {p ∈ M, σ ∈ E p }. We assume that the mesh is conformal, in the sense that M σ has exactly one element and then σ ⊂ ∂Ω or M σ has two elements and σ ⊂ Ω. We then denote by E p the faces of p ∈ M (it has d + 1 elements) and by x σ the centre of gravity of σ. We then denote by d pσ = |x σ − x p | the orthogonal distance between x p and σ ∈ E p and by n p,σ the normal vector to σ, outward to p.
3. We denote by V p the set of all the vertexes of p ∈ M (it has d + 1 elements), by V the union of all V p , p ∈ M. For y ∈ V, we denote by K p,y the polyhedron, defined as the set of all x ∈ p such that the barycentric coordinates ( , such that s y ′ ≥ 0 and y ′ ∈Vp s y ′ = 1). We denote by V σ the set of all vertexes of σ ∈ E (it has d elements), and by E p,y the set of all σ ∈ E p such that y ∈ V σ (it has d elements). We then denote, for σ ∈ E and y ∈ V σ , by x σ,y the point of σ defined by the barycentric coordinates (in σ) (s y ′ ) y ′ ∈Vσ , such that s y ′ = 1/(d + 1) for all y ′ ∈ V σ \ {y} (therefore s y = 2/(d + 1)).

5.
We denote, for all u ∈ H D , by Π D u ∈ L 2 (Ω) the function defined by the constant value u p a.e. in p ∈ M.
6. We denote, for all v ∈ H 1 0 (Ω), by P D v ∈ H D the element defined by (P D v) p = 1 |p| p v(x)dx for all p ∈ M. We denote by σ y the subset of all x ∈ σ such that the barycentric coordinates (s y ′ ) y ′ ∈Vσ of x in σ satisfy s y > 1/2, and we set (P D v) σ,y = d−1 (d+1)|σ| σ v(x)ds(x) + 2 (d+1)|σy| σy v(x)ds(x) for all σ ∈ E and y ∈ V σ (hence computing a second order approximation at point x σ,y ). 7. For u ∈ H D , p ∈ M and y ∈ V p , we denote by and by ∇ D u the function defined a.e. on Ω by ∇ p,y u on K p,y .
Remark 3.2 Note that the measure of K p,y is |p|/(d + 1) (this is easily shown, considering the affine transformation which sends p to a tetrahedron with all edges equal).
It can then be shown (see [2,12,13]) that a sequence of space discretisations defined as above, such that h Dm tends to 0, under a regularity property of the mesh, is an admissible sequence of space discretisations in the sense of Definition 2.3. Let us write the schemes (10) and (11) For a given interior σ common to p, q ∈ M, y ∈ V σ and n = 1, . . . , N T , we have µ(|∇ p,y u m |)∇ p,y u n · n p,σ = µ(|∇ q,y u m |)∇ q,y u n · n p,σ , and in this case no L ∞ estimate can be easily derived.

Properties of the schemes
Before focusing on the estimates satisfied by the approximate solutions, we first present a few properties which are useful in the convergence study.  (9) and (10). Then it holds:

Estimates and existence of a solution to the fully implicit scheme
Proof.

Convergence
Thanks to the estimates proved in the above section, we are now in position for proving the convergence of the scheme, using the monotonicity properties of the operators.

Convergence properties for the fully implicit scheme
We consider u ∈ H D,τ satisfying (9) and (10). We define Note that u D,τ is the solution of We then have the following convergence lemma.

it holds
Proof. We then remark that the sequence u m is bounded in L ∞ (0, T ; H Dm ), which provides, thanks to compactness property assumed in Definition 2.3, to the L 2 (Ω×]0, T [) bound on D τm u m and to an adaptation of Ascoli's theorem similar to that done in [7], that there existsū ∈ L ∞ (0, T ; H 1 0 (Ω)) ∩ C 0 (0, T ; L 2 (Ω)), withū t ∈ L 2 (Ω×]0, T [) such that, up to a sub-sequence, Π Dm,τm u m converges in L ∞ (0, T ; L 2 (Ω)) toū as m → ∞. We then get that D τm u m weakly converges in L 2 (Ω×]0, T [) toū t as m → ∞. The proof that u(., 0) = u 0 results from the initialisation of the scheme (9) and from Property (14). One of the difficulties is to respectively identifyḠ andw with µ(|∇ū|)∇ū and ν(ū, ∇ū). This will be done in further lemmas, thanks to the property (29) stated in the present lemma, that we have now to prove. Note that in the proof below, we drop some indexes m for the simplicity of notation. Let ϕ ∈ L 2 (0, T ; H 1 0 (Ω)) be given. Letting v = P D,τ ϕ in (28), and passing to the limit, we get Hence, setting ϕ =ū in (30), we get Passing to the limit in (28) with v = u m (the right hand side converges thanks to weak/strong convergence), we then get (29).

Convergence properties for the semi-implicit scheme
We consider u ∈ H D,τ , given by (9) and (11). We define and G D,τ defined by (27). Note that u D,τ is the solution of We then have the following convergence lemma.

relation (29) holds.
Proof. The proof mainly follows the same steps as that of Lemma 5.1. Let us focus on the points which are specific. Writing tends to 0 as m → ∞ thanks to (35) and properties of function µ. The same holds for T 0 Ω ( G Dm,τm (x, t)− G Dm,τm (x, t)) · ∇ Dm u m (x, t)dxdt, which proves (34).

Strong convergence of ∇ D u
The problem is now to show the convergence in L 2 (Ω×]0, T [) of ∇ Dm u m to ∇ū. This will result from property (29) (which holds for both the fully implicit and the semi-implicit schemes), and from the properties of function µ. Indeed, this property is the key point of the proof of the following lemma which uses Minty's trick.
Then the following holds and thereforeḠ Proof. In order to pass to the limit in T m (W ), we write T m (W ) = T Hence we get (37), which is sufficient to prove next Lemma 5.4. Nevertheless, let us apply Minty's trick (which remains available in the framework of non strictly monotonous operators): we set W = ∇ū − λψ, with λ > 0 and ψ ∈ C ∞ c (Ω×]0, T [) d in (37). We get, dividing by λ, We can let λ −→ 0 in the above inequality, using Lebesgue's dominated convergence theorem. We then get Since this also holds for −ψ, we get HenceḠ − µ(|∇ū|∇ū) = 0 a.e. in Ω×]0, T [, which achieves the proof of (38).
We now have the following lemma. Proof. We first remark that, for all V, W ∈ L 2 (Ω×]0, T [) d , it holds Indeed, thanks to the Cauchy-Schwarz inequality, we get  We can now conclude the convergence of the scheme. Proof. We first apply Lemmas 5.1 or 5.2. We get (28) or (33). We apply Lemma 5.4. We thus get thatw = f − ν(ū, ∇ū)ū t a.e. in Ω×]0, T [, which, in addition to (38), concludes the proof.

Numerical experiments
In this section we present several numerical examples to illustrate the properties of the proposed numerical schemes. They are devoted to the solution of regularised mean curvature flow and to the motion of 2D curves by curvature in level set formulation. In all examples we use both the semi-implicit and fully implicit schemes. We compute the errors and experimental order of convergence (EOC) for the whole level set function and also for the level set representing moving curve. Let us emphasise that for all the tests also proposed in [7], the order of the obtained errors and EOC are very close, using similar space and time steps. In the tables below n is number of finite volumes along each boundary side and n 2 is a total number of finite volumes. We consider square domain Ω to the equation u t |∇u| 2 + 0.5 − div ∇u  Table 1 for first numerical scheme and in Table 3 for the second one. For the fully implicit scheme and T OL = 10 −12 , we need about 20-30 Gauss-Seidel iterations and about 40-50 nonlinear iterations. The results are presented in Table 2 and Table 4. As one can see, the fully implicit schemes are about 10 times more accurate than the semi-implicit ones in this smooth example, but clearly they are more computationally complex and thus resulting in longer CPU times. All the schemes have EOC = 2 in solution error and EOC is slightly better than 1 in the gradient.     [14] u(x, y, t) = min{ to the level set equation   In this example, the solution contains flat regions and so we have to use the Evans-Spruck type regularisation [5]. Moreover, a singular circular curve with gradient jump is present, so we cannot expect second order accuracy in numerical solution. However, as we see from the Tables 5-8, the numerical schemes converge also in this singular case and naturally, EOC is equal (or close to) 1 for the solution error. In order to mimic convergence of numerical solution to (43) we use the coupling a = h 2 . One can also observe that the error obtained using the fully implicit schemes is only slightly better than that provided by the semi-implicit ones. This may lead to the conclusion that concerning both precision and CPU time aspects, one should use semi-implicit schemes as a reasonable compromise (cf. [10,4]   Example 3. In this example we compute the displacement of the unit circle by its curvature, and we compare the numerical results with the exact solution. The exact radius r(t) of a shrinking circle can be analytically expressed by The initial condition is given by u 0 (x, y) = −1 + x 2 + y 2 , which represents the distance function to the initial unit circle. Since we use zero Neumann boundary conditions in this example, the initial level set function is deformed, see Figure 3, but the error on the interface decreases with respect to the space and time steps, indicating the convergence of the method, as it can be seen in Tables 9 and 10. The measurement of the error is similar to that of [9]. For every discrete time step k = 0, 1, . . . , N , we first compute all the points x k i , i = 1, 2, . . . , P where the piecewise linear representation of the numerical solution becomes equal to zero along the finite element grid lines. Then we compute the distances r k i , i = 1, 2, . . . , P between the origin and the points x k i , i = 1, 2, . . . , P . Finally, these distances are compared to the radius r(kτ ) of the exact evolving circle. Then the formula is used for assessing the error in the L 2 (0, T ; L 2 (S 1 )) norm, denoting by S 1 the unit circle and setting T = N τ . The results for the numerical schemes on rectangular grids in semi-implicit and fully implicit versions are summarised in Tables 9 and 10. One can see that the fully implicit scheme (18) -(20) has the second order behaviour in this case. In Figure 4 we represent the numerical evolution of a circle together with the exact solution, setting n = 80, τ = h 2 and using 400 time steps. We hardly distinguish in this figure the numerical solution and the exact one.

Conclusions
The family of discrete schemes presented in this paper shows very easy implementation properties, and satisfactory accuracy. The adaptation of the viscosity solution sense to this discrete framework remains an open problem.