# Applications of approximate gradient schemes for nonlinear parabolic equations

Robert Eymard, Angela Handlovicova, Raphaèle Herbin, Karol Mikula, O.

Stašová

► To cite this version:

Robert Eymard, Angela Handlovicova, Raphaèle Herbin, Karol Mikula, O. Stašová. Applications of approximate gradient schemes for nonlinear parabolic equations. Applications of Mathematics, 2015, 10.1007/s10492-015-0088-4. http://dx.doi.org/10.1007/s10492-015-0074-7

HAL Id: hal-00744767

https://hal.science/hal-00744767

Submitted on 23 Oct 2012

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Applications of approximate gradient schemes for nonlinear parabolic equations

R. Eymard,\* A. Handlovičová,† R. Herbin‡, K. Mikula† and O. Stašová†

## Abstract

We propose the use of gradient schemes for the approximation of possibly non-local regularizations of the Perona-Malik equations. We then prove the convergence of this method thanks to compactness arguments, and we apply it to a particular gradient scheme on rectangular meshes, whose advantage is to lead to nine-point stencil matrices. Numerical examples provide indications about the accuracy of the method.

**Keywords:** regularized Perona-Malik equation, gradient schemes.

## 1 Introduction

A large class of image processing methods is based on the use of approximate solutions to equations of the following type

$$u_t - \operatorname{div} (G(u, x, t)\nabla u) = r(x, t), \text{ for a.e. } (x, t) \in \Omega \times (0, T) \quad (1)$$

with initial condition:

$$u(x, 0) = u_{\text{ini}}(x), \text{ for a.e. } x \in \Omega, \quad (2)$$

and homogeneous Neumann boundary condition:

$$G(u, x, t)\nabla u(x, t) \cdot \mathbf{n}_{\partial\Omega}(x) = 0, \text{ for a.e. } (x, t) \in \partial\Omega \times \mathbb{R}_{+}, \quad (3)$$

where the following hypotheses, called Hypotheses (H) in this paper, are considered:

* $\Omega$ is an open bounded polyhedron in $\mathbb{R}^d$, $d \in \mathbb{N}^*$, with boundary $\partial\Omega$,

* $T > 0$, $u_{\text{ini}} \in L^2(\Omega)$, $r \in L^2(\Omega \times (0, T))$,

* the possibly non local function $G$ is such that:

$$\begin{align*}
G: L^2(\Omega) \times \Omega \times (0, T) &\to \mathcal{L}(\mathbb{R}^d, \mathbb{R}^d), \\
G(\cdot, x, t) &\text{ is continuous for a.e. } (x, t) \in \Omega \times (0, T), \\
G(u, \cdot, \cdot) &\text{ is measurable for all } u \in L^2(\Omega), \\
G(u, x, t) &\text{ is self-adjoint with eigenvalues in } (\lambda, \bar{\lambda}) \\
&\text{ for all } u \in L^2(\Omega) \text{ and for a.e. } (x, t) \in \Omega \times (0, T),
\end{align*} \quad (4)$$

denoting by $\mathcal{L}(\mathbb{R}^d, \mathbb{R}^d)$ the set of linear mappings from $\mathbb{R}^d$ to $\mathbb{R}^d$.

<footnote>\*Université Paris-Est, France, robert.eymard@univ-mlv.fr</footnote>

<footnote>†Slovak University of Technology, angela,mikula,storova@math.sk</footnote>

<footnote>‡Aix-Marseille Université, France, raphaele.herbin@cmi.univ-mrs.fr</footnote>
In image processing applications, $u_{\text{ini}}$ represents an original noisy image, the solution $u(x, t)$ represents its filtering which depends on the scale parameter $t$. The space dimension $d$ is equal to 2 for 2D image filtering, 3 for 3D image or 2D+time movie filtering and 4 for 3D+time filtering of spatio-temporal image sequences. The image processing methods based on approximations of equation (1) differ by the definition of the function $G$. In this paper, we consider cases where $G$ arises from some regularization of the Perona-Malik equation [17], which reads

$$ \partial_t u - \operatorname{div}(g(|\nabla u|)\nabla u) = 0, \qquad (5) $$

where

$$ g(s) = \frac{1}{1 + Ks^2}, \quad \forall s \in \mathbb{R}^{+}, \qquad (6) $$

for a given $K > 0$. Recall that the mapping $s \mapsto sg(s)$ is not monotonously increasing on $\mathbb{R}^{+}$, and therefore the Perona-Malik equation is an ill-posed parabolic problem on general initial data.

The convergence of a numerical scheme for the one-dimensional original Perona-Malik problem (5) is proved in [4]. The analysis of a finite element discretization of a modified Perona-Malik equation is performed in [3]; in this latter work, the function $g$ depends on the $x$ and $y$ derivatives of $u$, rather than on the norm of the full gradient, which is the case considered here.

We first consider a regularization of the example due to Catté, Lions, Morel and Coll [6]; this examples reads

$$ \partial_t u - \operatorname{div}(G(u(\cdot, t), \cdot)\nabla u) = 0 \qquad (7) $$

with

$$ G(u, x) = \max(g\left(\left|\int_{\Omega} u(y)\nabla G_{\sigma}(x-y)dx\right|\right), \alpha), \quad \forall u \in L^{1}(\Omega), \qquad (8) $$

where $\alpha > 0$ is a given small value (in [6], although the parameter $\alpha$ is not introduced, a similar bound by below is obtained) and $G_{\sigma} \in C^{\infty}(\mathbb{R}^d)$ is a smoothing kernel, e.g. the Gauss function or mollifier with a compact support, for which $\int_{\mathbb{R}^d} G_{\sigma}(x)dx = 1$. Thanks to the convolution properties, the nonlinearity in the diffusion term depends on the unknown function $u$, contrary to the original Perona-Malik equation (without convolution) where it depends on the gradient of the solution. Note that, since the function $s \to \max(g(s), \alpha)$ is Lipschitz continuous with some constant $L_g$, we get that, for all $u, v \in L^1(\Omega)$,

$$ \begin{aligned} |G(u, x) - G(v, x)| &\le L_g \left|\int_{\Omega} u(y)\nabla G_{\sigma}(x-y)dy\right| - \left|\int_{\Omega} v(y)\nabla G_{\sigma}(x-y)dx\right| \\ &\le L_g \left|\int_{\Omega} (u(y) - v(y))\nabla G_{\sigma}(x-y)dy\right| \le L_g \|\nabla G_{\sigma}\|_{\infty} \|u - v\|_{L^1(\Omega)}, \end{aligned} $$

which shows that (8) enters into the framework defined by Hypotheses (H). In the case of the regularized model, the convergence of classical finite volume schemes was proved in [16], [7].

Another interesting image processing model with the structure of equation (1) is the so-called nonlinear tensor anisotropic diffusion introduced by Weickert [18]. In that case, the linear mapping $G(u, x, t)$ represents the so-called diffusion tensor depending on the eigenvalues and eigenvectors of the (regularized) structure tensor with matrix

$$ J_{\rho}(\nabla u_{\sigma}) = G_{\rho} * (\nabla u_{\sigma} \nabla u_{\sigma}^{T}), \qquad (9) $$

where $f * g$ denotes the convolution between the two functions $f$ and $g$, $u_{\sigma}$ is defined by

$$ u_{\sigma}(x, t) = (G_{\sigma} * u(\cdot, t))(x) \qquad (10) $$

and $G_{\sigma}$ and $G_{\rho}$ are Gaussian kernels. In computer vision, the matrix $J_{\rho} = \begin{pmatrix} a & b \\ b & c \end{pmatrix}$, which is symmetric and positive semidefinite, is also known as the interest operator or second moment matrix. We may write
$a = G_{\rho} * (\partial_1 G_{\sigma} * u)^2$, $b = G_{\rho} * ((\partial_1 G_{\sigma} * u)(\partial_2 G_{\sigma} * u))$ and $c = G_{\rho} * (\partial_2 G_{\sigma} * u)^2$. The orthogonal set of eigenvectors $(v, w)$ of $J_{\rho}$ corresponding to its eigenvalues $(\mu_1, \mu_2)$, $\mu_1 \ge \mu_2$, is such that the orientation of the eigenvector $w$, which corresponds to the smaller eigenvalue $\mu_2$, gives the so-called coherence orientation. This orientation has the lowest fluctuations in image intensity. The diffusion tensor $G$ in equation (1) is then designed to steer a smoothing process such that the filtering is strong along the coherence direction $w$ and increasing with the coherence defined by the difference of the eigenvalues $(\mu_1 - \mu_2)^2$. To that goal, $G$ must possess the same eigenvectors $v = (v_1, v_2)$ and $w = (-v_2, v_1)$ as the structure tensor $J_{\rho}(\nabla u_{\sigma})$ and the eigenvalues of $G$ can be chosen as follows

$$ \kappa_1 = \alpha, \quad \alpha \in (0, 1), \quad \alpha \ll 1, \qquad (11) $$

$$ \kappa_2 = \begin{cases} \alpha, & \text{if } \mu_1 = \mu_2, \\ \alpha + (1 - \alpha) \exp\left(\frac{-C}{(\mu_1 - \mu_2)^2}\right), & \text{otherwise.} \end{cases} \quad C > 0 \quad \text{otherwise.} $$

So, the matrix $M_G$ of the linear operator $G(u, x, t)$ is finally defined by

$$ M_G = ABA^{-1}, \quad \text{where} \quad A = \begin{pmatrix} v_1 & -v_2 \\ v_2 & v_1 \end{pmatrix} \quad \text{and} \quad B = \begin{pmatrix} \kappa_1 & 0 \\ 0 & \kappa_2 \end{pmatrix}. \qquad (12) $$

Thanks to this construction by convolutions, the diffusion matrix (nonlinearly) depends on the solution $u$ and it satisfies smoothness, symmetry and uniform positive definiteness properties. It is then possible to show that it enters into Hypotheses (H) by similar computations to those done above and in [9]. The so-called diamond-cell finite volume schemes for the nonlinear tensor anisotropic diffusion were suggested and analyzed in [9, 8].

Another type of regularization of the classical Perona-Malik approach is obtained by considering the gradient information from delayed time $t - \bar{t}$, for a given $\bar{t} > 0$. We call this model the time-delayed Perona-Malik equation: consider (1) with $u_{\text{ini}} \in H^1(\Omega)$, and define $u(x, t) = u_{\text{ini}}(x)$ for $x \in \Omega$ and $t < 0$. Then, for any $k \in \mathbb{N}$, we define the function $G$ in the time interval $]k\bar{t}, (k+1)\bar{t}]$, by

$$ G(x, t) = \max(g(|\nabla u(x, t - \bar{t})|), \alpha), \qquad (13) $$

where $\alpha > 0$ is again a small parameter. Then problem (1) boils down to a standard linear parabolic problem on $]k\bar{t}, (k+1)\bar{t}]$ (which is included in the framework of Hypotheses (H)), and, as shown in the numerical part of this paper, it leads to an efficient approximation of the Perona-Malik equation. An integral average on a given shifted time interval is used in [2] for regularization purposes. Our regularization can be understood as a discrete approximation of the integral average, e.g. by a mid-point rule. Note that the Amann theory holds for the original integral regularization without numerical integration while our theory and convergence proofs are valid for the gradient approximation shifted backward in time.

This paper is organized as follows. After providing the sense for a weak solution in Section 2, we define the class of gradient schemes [13, 12], to approximate the models based on equation (1). We then prove their convergence in Section 3. Then in Section 4, we study a particular scheme among this class, which is particularly efficient for image processing. Finally, numerical results are given in Section 5, confirming this efficiency, and are followed by a short conclusion.

## 2 The weak formulation

**Definition 2.1 (Weak solution to (1)-(2)-(3))** Under Hypotheses (H), a function $u$ is a weak solution of (1)-(2)-(3) if, for all $T > 0$,

$$ 1. u \in L^2(0, T; H^1(\Omega)), $$
2. the following holds

$$ \begin{aligned} & \int_0^T \int_{\Omega} (-u(x,t)\varphi_t(t)w(x) + \varphi(t)G(u,x,t)\nabla u(x,t) \cdot \nabla w(x)) \,dxdt - \int_{\Omega} u_{\text{ini}}(x)w(x)\varphi(0)dx \\ & = \int_0^T \int_{\Omega} r(x,t)w(x)\varphi(t)dxdt, \quad \forall w \in H^1(\Omega), \quad \forall \varphi \in C_c^\infty([0,T]), \end{aligned} \quad (14) $$

where we denote by $C_c^\infty([0, T])$ the set of functions of $C_c^\infty([0, \infty, T])$ restricted to $[0, T]$.

We have the following standard result [5].

**Theorem 2.1 (Properties of a weak solution u to (1)-(2)-(3)) Under Hypotheses (H), the function u is a weak solution to (1)-(2)-(3) in the sense of Definition 2.1 if and only if**

1. $u \in L^2(0, T; H^1(\Omega))$, $u_t \in L^2(0, T; H^1(\Omega))'$ (defining the standard continuous embedding of $L^2(\Omega)$ in $H^1(\Omega)'$), and therefore $u \in C^0([0, T]; L^2(\Omega))$,

2. $u(x, 0) = u_{\text{ini}}(x)$ for a.e. $x \in \Omega$,

3. the following holds

$$ \begin{aligned} & \int_0^T \left( \langle u_t(t), v(t) \rangle_{H^1(\Omega)', H^1(\Omega)} + \int_{\Omega} G(u, x, t) \nabla u(x, t) \cdot \nabla v(x, t) dx \right) dt \\ & = \int_0^T \int_{\Omega} r(x, t)v(x, t) dxdt, \quad \forall v \in L^2(0, T; H^1(\Omega)). \end{aligned} \quad (15) $$

Then u satisfies the following consequence of the above properties:

$$ \begin{aligned} & \frac{1}{2} \int_{\Omega} (u(x, t_0)^2 - u_{\text{ini}}(x)^2) dx + \int_{0}^{t_0} \int_{\Omega} G(u, x, t) \nabla u(x, t) \cdot \nabla u(x, t) dx dt \\ & = \int_{0}^{t_0} \int_{\Omega} r(x, t) u(x, t) dx dt, \quad \forall t_0 \in [0, T]. \end{aligned} \quad (16) $$

# 3 Approximate gradient schemes for parabolic equations

## 3.1 Definition of the scheme

**Definition 3.1 (Approximate gradient discretization and gradient scheme)** Let $\Omega$ be a bounded open domain of $\mathbb{R}^d$, with $d \in \mathbb{N}^*$. An approximate gradient discretization $\mathcal{D}$ is defined by $\mathcal{D} = (X_{\mathcal{D}}, \Pi_{\mathcal{D}}, \nabla_{\mathcal{D}})$, where:

1. the set of discrete unknowns $X_{\mathcal{D}}$ is a finite dimensional vector space on $\mathbb{R}$,

2. the mapping $\Pi_{\mathcal{D}} : X_{\mathcal{D}} \to L^2(\Omega)$ is the reconstruction of the approximate function,

3. the mapping $\nabla_{\mathcal{D}} : X_{\mathcal{D}} \to L^2(\Omega)^d$ is the reconstruction of the gradient of the function,

4. $\|u\|_{\mathcal{D}} = (\|\Pi_{\mathcal{D}}u\|_{L^2(\Omega)}^2 + \|\nabla_{\mathcal{D}}u\|_{L^2(\Omega)^d}^2)^{1/2}$ is a norm on $X_{\mathcal{D}}$.

Then the **compactness** of the discretization is measured through the function $T_{\mathcal{D}} : \mathbb{R}^d \to \mathbb{R}^+$, defined by

$$ T_{\mathcal{D}}(\xi) = \max_{v \in X_{\mathcal{D}} \setminus \{0\}} \frac{\|\Pi_{\mathcal{D}}v(\cdot + \xi) - \Pi_{\mathcal{D}}v\|_{L^2(\mathbb{R}^d)}}{\|v\|_{\mathcal{D}}}, \quad \forall \xi \in \mathbb{R}^d, \quad (17) $$
where $\Pi_D v$ is set to zero outside of $\Omega$. Note that $\lim_{|\xi| \to 0} T_D(\xi) = 0$.

The strong consistency of the discretization is measured through the interpolation error function $S_D: H_0^1(\Omega) \to [0, +\infty)$, defined by

$$S_D(\varphi) = \min_{v \in X_D} \left( \|\Pi_D v - \varphi\|_{L^2(\Omega)}^2 + \|\nabla_D v - \nabla \varphi\|_{L^2(\Omega)^d}^2 \right)^{\frac{1}{2}}, \quad \forall \varphi \in H^1(\Omega), \quad (18)$$

The limit conformity of the discretization is measured through the conformity error function $W_D: H_{\text{div},0}(\Omega) \to [0, +\infty)$ (where $H_{\text{div},0}(\Omega)$ denotes the set of all elements of $H_{\text{div}}(\Omega)$ with zero normal trace), defined by

$$W_D(\varphi) = \max_{u \in X_D \setminus \{0\}} \frac{1}{\|u\|_D} \left| \int_{\Omega} (\nabla_D u(x) \cdot \varphi(x) + \Pi_D u(x) \text{div} \, \varphi(x)) \, dx \right|, \quad \forall \varphi \in H_{\text{div},0}(\Omega). \quad (19)$$

**Definition 3.2 (Space-time discretization)** Let $\Omega$ be an open subset of $\mathbb{R}^d$, with $d \in \mathbb{N}^*$ and let $T > 0$ be given. We say that $(D, \tau)$ is a space-time gradient discretization of $\Omega \times (0, T)$ if

* $D = (X_D, \Pi_D, \nabla_D)$ is an approximate gradient discretization of $\Omega$ in the sense of Definition 3.1,

* there exists $N_T \in \mathbb{N}$ with $T = N_T \tau$, where $\tau > 0$ is the time step.

We then define $X_{D,\tau} = \{(u^n)_{n=1,...,N_T}, u^n \in X_D\}$, and we define the mappings $\Pi_{D,\tau}: X_{D,\tau} \to L^2(\Omega \times (0,T))$ and $\nabla_{D,\tau}: X_{D,\tau} \to L^2(\Omega \times (0,T))^d$ by

$$\Pi_{D,\tau} u(x,t) = \Pi_D u^n(x), \quad \text{for a.e. } x \in \Omega, \forall t \in (n-1)\tau, n\tau, \forall n = 1, \dots, N_T, \quad (20)$$

and

$$\nabla_{D,\tau} u(x,t) = \nabla_D u^n(x), \quad \text{for a.e. } x \in \Omega, \forall t \in (n-1)\tau, n\tau, \forall n = 1, \dots, N_T. \quad (21)$$

Let $(D, \tau)$ be a space-time discretization of $\Omega \times (0, T)$. We define the following scheme for the discretization of Problem (1):

$$\begin{align}
u \in X_{D,\tau}, \quad D_{\tau} u(x,t) &= \frac{1}{\tau} (\Pi_D u^1(x) - u_{\text{ini}}(x)), \quad \text{for a.e. } x \in \Omega, \forall t \in [0, \tau], \tag{22} \\
D_{\tau} u(x,t) &= \frac{1}{\tau} (\Pi_D u^n(x) - \Pi_D u^{n-1}(x)), \quad \text{for a.e. } x \in \Omega, \forall t \in (n-1)\tau, n\tau, \forall n = 2, \dots, N_T, \nonumber
\end{align}$$

and

$$\int_0^T \int_{\Omega} (D_{\tau} u \Pi_{D,\tau} v + G_{D,\tau} (\Pi_{D,\tau} u, x, t) \nabla_{D,\tau} u \cdot \nabla_{D,\tau} v) \, dxdt = \int_0^T \int_{\Omega} r \Pi_{D,\tau} v \, dxdt, \quad \forall v \in X_{D,\tau}, \quad (23)$$

where

$$\begin{align}
G_{D,\tau}: L^2(\Omega) \times \Omega \times (0,T) &\to \mathcal{L}(\mathbb{R}^d, \mathbb{R}^d), \text{ there exists } C_{D,\tau} > 0 \text{ with} \nonumber \\
&\|G_{D,\tau}(v, x, t) - G(v, x, t)\|_{\mathcal{L}(\mathbb{R}^d, \mathbb{R}^d)} \le C_{D,\tau} \|v\|_{L^2(\Omega)}, \forall v \in L^2(\Omega), \text{ for a.e. } (x, t) \in \Omega \times (0, T), \nonumber \\
G_{D,\tau}(u, \cdot, \cdot) \text{ is measurable for all } u \in L^2(\Omega \times (0, T)), \nonumber \\
G_{D,\tau}(u, x, t) \text{ is self-adjoint with eigenvalues in } (\lambda, \bar{\lambda}) \nonumber \\
&\text{ for all } u \in L^2(\Omega) \text{ and for a.e. } (x, t) \in \Omega \times (0, T). \tag{24}
\end{align}$$

We then denote $\Pi_{D,\tau} u(\cdot, 0) = u_{\text{ini}}$, hence defining $\Pi_{D,\tau} u(\cdot, t)$ for all $t \in [0, T]$.

**Remark 3.1** Note that Assumption (24) holds in particular for $G_{D,\tau} = G$, which occurs in some of the numerical applications.
## 3.2 Properties of the scheme

**Lemma 3.1** ($L^\infty(0,T;L^2(\Omega))$ and $L^2(0,T;H^1(\Omega))$ estimates and existence of a discrete solution)
Under Hypotheses (H), let $(\mathcal{D}, \tau)$ be a space-time gradient discretization of $\Omega \times (0,T)$ in the sense of Definition 3.2. Then there exists at least one solution to Scheme (22)-(23), which moreover satisfies that there exists a constant $C_1 > 0$ such that:

$$ \|\Pi_{\mathcal{D}} u^m\|_{L^2(\Omega)}^2 \le C_1 (\|u_{\text{ini}}\|_{L^2(\Omega)}^2 + T \|r\|_{L^2(\Omega \times (0,T))}^2), \quad \forall m = 1, \dots, N_T, \qquad (25) $$

and

$$ \|\nabla_{\mathcal{D},\tau} u\|_{L^2(\Omega \times (0,T))}^2 \le \frac{C_1}{\lambda} (\|u_{\text{ini}}\|_{L^2(\Omega)}^2 + T \|r\|_{L^2(\Omega \times (0,T))}^2). \qquad (26) $$

**PROOF.** Let $m = 1, \dots, N_T$ and let us set $v^n = u^n$ for $n = 1, \dots, m$ and $v^n = 0$ for $n = m + 1, \dots, N_T$ in (23). We obtain, thanks to the equality $a(a-b) = \frac{1}{2}(a^2 + (a-b)^2 - b^2)$,

$$ \begin{aligned} & \int_{\Omega} \left( \frac{1}{2} (\Pi_{\mathcal{D}} u^m(x)^2 - u_{\text{ini}}(x)^2) + \lambda \int_{0}^{m\tau} |\nabla_{\mathcal{D},\tau} u(x,t)|^2 dt \right) dx \\ & \le \int_{0}^{m\tau} \int_{\Omega} r(x,t) \Pi_{\mathcal{D},\tau} u(x,t) dx dt. \end{aligned} \qquad (27) $$

Applying the Young inequality to the right hand side provides

$$ \begin{aligned} & \int_{\Omega} \left( \frac{1}{2} (\Pi_{\mathcal{D}} u^m(x)^2 - u_{\text{ini}}(x)^2) + \lambda \int_{0}^{m\tau} |\nabla_{\mathcal{D},\tau} u(x,t)|^2 dt \right) dx \\ & \le T \int_{0}^{m\tau} \int_{\Omega} r(x,t)^2 dx dt + \frac{1}{4T} \int_{0}^{m\tau} \int_{\Omega} \Pi_{\mathcal{D},\tau} u(x,t)^2 dx dt. \end{aligned} $$

We now give a discrete version of the Gronwall lemma. Defining $a_m = \int_{\Omega} \Pi_{\mathcal{D}} u^m(x)^2 dx$, $b = \int_{\Omega} u_{\text{ini}}(x)^2 dx + 2T \int_{0}^{T} \int_{\Omega} r(x,t)^2 dx dt$, we get

$$ a_m \le b + \frac{1}{2T} d_m \quad \text{with} \quad d_m = \tau \sum_{n=1}^{m} a_n, $$

which can be written, setting $c = \tau/2T \le 1/2$ and $d_0 = 0$,

$$ d_m - d_{m-1} \le \tau b + cd_m. $$

This gives $d_m \le (d_{m-1} + \tau b)/(1-c)$, and therefore $d_m \le \tau \sum_{n=1}^{m} b/(1-c)^{m+1-n}$. Using the inequality $mc \le 1/2$, which leads to $-\log(1-c) \le \log(2m/(2m-1))$, we get

$$ \log \frac{1}{(1-c)^m} \le m \log(\frac{2m}{2m-1}) \le m(\frac{2m}{2m-1} - 1) = \frac{m}{2m-1} \le 1, $$

and therefore $1/(1-c)^m \le e$. We have

$$ \tau \sum_{n=1}^{m} \frac{1}{(1-c)^{m+1-n}} \le \tau \frac{\frac{1}{(1-c)^{m+1}} - 1}{\frac{1}{1-c} - 1} \le \tau \frac{1}{\tau/2T} \frac{1}{(1-c)^m} \le 2Te, $$

which gives $d_m \le 2Teb$. Hence we conclude $a_m \le (1+2e)b$, which concludes (25) and (26). The topological degree argument applied to numerical schemes (see e.g. [14]) allows to conclude to the existence of at least one solution to the scheme. $\square$

Note that the minimal regularity of the initial data corresponding to the framework of image processing, and the generality of the operator G prevent from obtaining more regular estimates than those of Lemma 3.1. In particular, multiplying the scheme with discrete time derivative does not seem to yield a better estimate.

We may now state the convergence theorem.
**Theorem 3.1** Let Hypotheses (H) be fulfilled. Let $(D_m, \tau_m)_{m \in \mathbb{N}}$ be a sequence of space-time discretizations of $\Omega \times (0, T)$ in the sense of Definition 3.2 such that $C_{D_m, \tau_m}$ (introduced in (24)) and $\tau_m > 0$ tend to 0 as $m \to \infty$. We assume that there exists a function $T: \mathbb{R}^d \to \mathbb{R}^+$ with $\lim_{|\xi| \to 0} T(\xi) = 0$ such that $T_{D_m} \le T$ for all $m \in \mathbb{N}$, and that

$$ \forall \varphi \in H^1(\Omega), \lim_{m \to \infty} S_{D_m}(\varphi) = 0, \qquad (28) $$

and

$$ \forall \varphi \in H_{\text{div},0}(\Omega), \lim_{m \to \infty} W_{D_m}(\varphi) = 0. \qquad (29) $$

Let, for all $m \in \mathbb{N}$, $u_m \in X_{D_m, \tau_m}$ be such that (22)-(23) hold.

Then there exists a weak solution $u$ to (1)-(2)-(3) in the sense of Definition 2.1 such that, up to a subsequence, $\Pi_{D_m, \tau_m} u_m$ tends to $u$ in $L^2(\Omega \times (0, T))$, and $\nabla_{D_m} u_{D_m, \tau_m}$ tends to $\nabla u$ for the weak topology of $L^2(\Omega \times (0, T))^d$. In the case that $u$ is unique, then the whole sequence converges in the same sense.

**PROOF.** Thanks to (25)-(26), we first have the existence of some $C_2 > 0$, independent of $m$, such that

$$ \|\Pi_{D_m, \tau_m} u_m(\cdot + \xi, \cdot) - \Pi_{D_m, \tau_m} u_m\|_{L^2(\mathbb{R}^d \times (0,T))} \le C_2 T(\xi). $$

Let us now prove an estimate on the time translates. Let $\eta \in (0, T)$. We drop the index $m$ for the simplicity of notation. We consider

$$ A(\eta) = \|\Pi_{D,\tau} u(\cdot, \cdot + \eta) - \Pi_{D,\tau} u\|_{L^2(\Omega \times (0,T-\eta))}^2. $$

Denoting by $\nu(t)$ the integer $n$ such that $t \in ((n-1)\tau, n\tau]$, we have $A(\eta) = \int_0^{T-\eta} \sum_{n=\nu(t)}^{\nu(t+\eta)} \tau(A_1^n - A_2^n)dt$,

with

$$ A_1^n = \frac{1}{\tau} \int_{\Omega} (\Pi_{D,\tau} u^n - \Pi_{D,\tau} u^{n-1}) \Pi_D u^{\nu(t+\eta)} dx $$

and

$$ A_2^n = \frac{1}{\tau} \int_{\Omega} (\Pi_{D,\tau} u^n - \Pi_{D,\tau} u^{n-1}) \Pi_D u^{\nu(t)} dx. $$

Using (23) with $v^n = u^{\nu(t+\eta)}$ and $v^p = 0$ for $p \ne n$, we get

$$ A_1^n = \frac{1}{\tau} \int_{(n-1)\tau}^{n\tau} \int_{\Omega} \left( r(x, s) \Pi_{D,\tau} u^{\nu(t+\eta)}(x) - G_{D,\tau} (\Pi_{D,\tau} u, x, s) \nabla_D u^n(x) \cdot \nabla_D u^{\nu(t+\eta)}(x) \right) dx ds. $$

Hence we get, using Young's inequality, that

$$ A_1^n \le \frac{1}{2} (A_{11}^n + A_{12}^{\nu(t+\eta)} + A_{13}^n + A_{13}^{\nu(t+\eta)}), $$

with

$$ A_{11}^n = \frac{1}{\tau} \int_{(n-1)\tau}^{n\tau} \int_{\Omega} r(x, s)^2 dx ds, \quad A_{12}^n = \int_{\Omega} (\Pi_{D,\tau} u^n)^2 dx, \quad A_{13}^n = \bar{\lambda} \int_{\Omega} (\nabla_D u^n)^2 dx. $$

A similar inequality holds for $A_2^n$. Hence, applying Lemma 6.1 stated in the appendix, and (25)-(26), we get the existence of $C_3$, which does not depend on $m$, such that

$$ A(\eta) \le C_3 \eta. $$

Thanks to (25) which is a $L^\infty(0, T; L^2(\Omega))$ bound on the approximate solution, it is easy to extend the time translates from $(0, T - \eta)$ to $\mathbb{R}$, for the approximate solution set to 0 outside of $\Omega \times (0, T)$. This proves that the time translates of $\Pi_{D_m, \tau_m} u_m$ uniformly tend to 0. Then, thanks to Kolmogorov's theorem, we may extract a subsequence such that $\Pi_{D_m, \tau_m} u_m$ converges in $L^2(\Omega \times (0, T))$ to some function
$u \in L^2(\Omega \times (0, T))$. Since $\nabla_{\mathcal{D}_m, \tau_m} u_m$ is bounded in $L^2(\Omega \times (0, T))^d$, it weakly converges in $L^2(\Omega \times (0, T))^d$, up to a subsequence, to some function $G$. For given $\varphi \in H_{\text{div},0}(\Omega)$ and $\varphi \in C_c^\infty((0, T))$, we can write, using (19),

$$ \left| \int_0^T \int_{\Omega} \varphi(t) (\nabla_{\mathcal{D}_m, \tau_m} u_m \cdot \varphi + \Pi_{\mathcal{D}_m, \tau_m} u_m \, \text{div} \, \varphi) \, dxdt \right| \le \int_0^T W_{\mathcal{D}_m}(\varphi) \varphi(t) \|u_m(t)\|_{\mathcal{D}_m} \, dt. $$

Letting $m \to \infty$ in the above inequality (which is possible thanks to bounds of the right hand side due to (25)-(26)) and using (29), we get

$$ \forall \varphi \in C_c^\infty((0, T)), \quad \forall \varphi \in H_{\text{div},0}(\Omega), \quad \int_0^T \int_{\Omega} \varphi(t) (G \cdot \varphi + u \, \text{div} \, \varphi) \, dxdt = 0. $$

This shows that $u \in L^2(0, T; H^1(\Omega))$ and that $G(x, t) = \nabla u(x, t)$ for a.e. $(x, t) \in \Omega \times (0, T)$.
Let us now prove that $u$ is a weak solution to (1)-(2)-(3) in the sense of Definition 2.1.
Let $\varphi \in C_c^\infty([0, T])$, and $w \in H^1(\Omega)$ be given. We denote by

$$ w_m = \underset{v \in X_{\mathcal{D}_m}}{\operatorname{argmin}} \left( \|\Pi_{\mathcal{D}_m} v - w\|_{L^2(\Omega)}^2 + \|\nabla_{\mathcal{D}_m} v - \nabla w\|_{L^2(\Omega)^d}^2 \right)^{\frac{1}{2}}. $$

Using (28), we then have that $\Pi_{\mathcal{D}_m} w_m$ converges to $w$ in $L^2(\Omega)$ and that $\nabla_{\mathcal{D}_m} w_m$ converges to $\nabla w$ in $L^2(\Omega)^d$. We then set $v^n = \varphi((n-1)\tau)w_m$ in (23), and we denote by $\varphi_m(t)$ the value $\varphi((n-1)\tau)$ for all $t \in [(n-1)\tau, n\tau]$. We get $T_1^m + T_2^m = T_3^m$, with

$$ T_1^m = \int_0^T \int_{\Omega} D_{\tau_m} u \, \Pi_{\mathcal{D}_m} w_m \varphi_m(t) \, dxdt, $$

$$ T_2^m = \int_0^T \varphi_m(t) \int_{\Omega} G_{\mathcal{D}_m, \tau_m} (\Pi_{\mathcal{D}_m, \tau_m} u, x, t) \nabla_{\mathcal{D}_m, \tau_m} u \cdot \nabla_{\mathcal{D}_m} w_m \, dxdt, $$

and

$$ T_3^m = \int_0^T \varphi_m(t) \int_{\Omega} r \, \Pi_{\mathcal{D}_m} w_m \, dxdt. $$

We can rewrite $T_1^m$ as

$$ T_1^m = - \int_{\Omega} u_{\text{ini}} \Pi_{\mathcal{D}_m} w_m \varphi(0) \, dx - \sum_{n=1}^{N_T} \int_{\Omega} \Pi_{\mathcal{D}_m} u^n \, \Pi_{\mathcal{D}_m} w_m (\varphi(n\tau) - \varphi((n-1)\tau)) \, dxdt, $$

which gives

$$ T_1^m = - \int_{\Omega} u_{\text{ini}} \Pi_{\mathcal{D}_m} w_m \varphi(0) \, dx - \int_0^T \int_{\Omega} \Pi_{\mathcal{D}_m, \tau_m} u \, \Pi_{\mathcal{D}_m} w_m \varphi_t(t) \, dxdt. $$

Hence we get that

$$ \lim_{m \to \infty} T_1^m = - \int_{\Omega} u_{\text{ini}} w \varphi(0) \, dx - \int_0^T \int_{\Omega} u \, w_m \varphi_t(t) \, dxdt. $$

Thanks to Assumption (24), we have, for a.e. $(x, t) \in \Omega \times (0, T)$:

$$ |G_{\mathcal{D}_m, \tau_m} (\Pi_{\mathcal{D}_m, \tau_m} u, x, t) - G(u, x, t)| \le C_{\mathcal{D}_m, \tau_m} \|\Pi_{\mathcal{D}_m, \tau_m} u(\cdot, t)\|_{L^2(\Omega)} + |G(\Pi_{\mathcal{D}_m, \tau_m} u, x, t) - G(u, x, t)|. $$

Thanks to the assumption that $C_{\mathcal{D}_m, \tau_m}$ tends to 0 as $m \to \infty$, and thanks to the convergence of $\Pi_{\mathcal{D}_m, \tau_m} u(\cdot, t)$ to $u(\cdot, t)$ in $L^2(\Omega)$ for a.e. $t \in (0, T)$, we get that $G_{\mathcal{D}_m, \tau_m}(\Pi_{\mathcal{D}_m, \tau_m} u, x, t)$ converges to $G(u, x, t)$ for a.e. $(x, t) \in \Omega \times (0, T)$. Noticing that:

$$ \int_{\Omega} G_{\mathcal{D}_m, \tau_m} (\Pi_{\mathcal{D}_m, \tau_m} u, x, t) \nabla_{\mathcal{D}_m, \tau_m} u \cdot \nabla_{\mathcal{D}_m} w_m \, dx = \\ \int_{\Omega} (G_{\mathcal{D}_m, \tau_m} (\Pi_{\mathcal{D}_m, \tau_m} u, x, t) - G(u, x, t)) \nabla_{\mathcal{D}_m, \tau_m} u \cdot \nabla_{\mathcal{D}_m} w_m \, dx + \int_{\Omega} G(u, x, t) \nabla_{\mathcal{D}_m, \tau_m} u \cdot \nabla_{\mathcal{D}_m} w_m \, dx, $$
we get by dominated convergence for the first term and weak-strong convergence for the second term, that:

$$ \lim_{m \to \infty} T_2^m = \int_0^T \varphi(t) \int_{\Omega} G(u, x, t) \nabla u \cdot \nabla w \, dx \, dt. $$

Finally, we also have:

$$ \lim_{m \to \infty} T_3^m = \int_0^T \varphi(t) \int_{\Omega} r w \, dx \, dt. $$

This achieves the proof that $u$ is a weak solution to (1)-(2)-(3) in the sense of Definition 2.1. $\square$

**Theorem 3.2** Under hypotheses of Theorem 3.1, let $(\mathcal{D}_m, \tau_m)_{m \in \mathbb{N}}$ be a subsequence satisfying the properties stated in the conclusions of Theorem 3.1. Then $\Pi_{\mathcal{D}_m, \tau_m} u_m(\cdot, t)$ tends to $u(t)$ in $L^2(\Omega)$ for all $t \in [0, T]$, and $\nabla_{\mathcal{D}_m} u_{\mathcal{D}_m, \tau_m}$ tends to $\nabla u$ in $L^2(\Omega \times (0, T))^d$. In the case that $u$ is unique, then the whole sequence converges in the same sense.

**PROOF.** Let $w \in H^1(\Omega)$ and let

$$ w_m = \underset{v \in X_{\mathcal{D}_m}}{\operatorname{argmin}} \left( \|\Pi_{\mathcal{D}_m} v - w\|_{L^2(\Omega)}^2 + \|\nabla_{\mathcal{D}_m} v - \nabla w\|_{L^2(\Omega)^d}^2 \right)^{\frac{1}{2}}. $$

Let us consider, for $s \le t \in [0, T]$,

$$ B(s,t) = \int_{\Omega} (\Pi_{\mathcal{D}_m, \tau_m} u_m(x,t) - \Pi_{\mathcal{D}_m, \tau_m} u_m(x,s)) \Pi_{\mathcal{D}_m} w_m(x) dx. $$

We get

$$ \begin{aligned} B(s,t) &= \int_{(\nu(t)-1)\tau}^{\nu(t+s)\tau} \int_{\Omega} D_{\tau_m} u_m \Pi_{\mathcal{D}_m} w_m(x) dx \\ &= \int_{(\nu(t)-1)\tau}^{\nu(t+s)\tau} \int_{\Omega} (r \Pi_{\mathcal{D}_m} w_m - G_{\mathcal{D}_m, \tau_m} (\Pi_{\mathcal{D}_m, \tau_m} u_m, x, t) \nabla_{\mathcal{D}_m, \tau_m} u \cdot \nabla_{\mathcal{D}_m} w_m) dx dt. \end{aligned} $$

Thanks to the Cauchy-Schwarz inequality, we get the existence of $C_4$, independent of $m$ and $w$, such that

$$ B(s,t) \le (\nu(t+s)\tau - (\nu(t)-1)\tau)^{1/2} C_4 \|w_m\|_{\mathcal{D}_m} \le (t-s+2\tau_m)^{1/2} C_4 \|w_m\|_{\mathcal{D}_m}. $$

This continuity property is sufficient to apply Theorem 6.1 (given in the Appendix), proving that, for all $t \in [0, T]$, $\Pi_{\mathcal{D}_m, \tau_m} u_m(t)$ tends to $u(t)$ for the weak topology of $L^2(\Omega)$. In the same way as in the proof of Lemma 3.1, we have the property, for a given $t_0 \in [0, T]$ (dropping some indices $m$),

$$ \begin{aligned} & \int_{\Omega} \left( \frac{1}{2} (\Pi_{\mathcal{D},\tau} u(x, t_0)^2 - u_{\text{ini}}(x)^2) + \int_{0}^{\nu(t_0)\tau} G_{\mathcal{D},\tau} (\Pi_{\mathcal{D},\tau} u, x, t) \nabla_{\mathcal{D},\tau} u \cdot \nabla_{\mathcal{D},\tau} u dt \right) dx \\ & \le \int_{0}^{\nu(t_0)\tau} \int_{\Omega} r(x,t) \Pi_{\mathcal{D},\tau} u(x,t) dx dt. \end{aligned} \quad (30) $$

Passing to the limit in the above inequality, we get

$$ \begin{aligned} & \limsup_{m \to \infty} \left( \frac{1}{2} \int_{\Omega} \Pi_{\mathcal{D},\tau} u(x, t_0)^2 + \int_{0}^{\nu(t_0)\tau} G_{\mathcal{D},\tau} (\Pi_{\mathcal{D},\tau} u, x, t) \nabla_{\mathcal{D},\tau} u \cdot \nabla_{\mathcal{D},\tau} u dt \right) dx \\ & \le \frac{1}{2} \int_{\Omega} u_{\text{ini}}(x)^2 dx + \int_{0}^{t_0} \int_{\Omega} r(x,t) u(x,t) dx dt. \end{aligned} $$

Since $u$ satisfies (16), we have

$$ \frac{1}{2} \int_{\Omega} u(x, t_0)^2 dx + \int_{0}^{t_0} \int_{\Omega} G(u, x, t) \nabla u \cdot \nabla u dx dt = \frac{1}{2} \int_{\Omega} (u_{\text{ini}})^2 dx + \int_{0}^{t_0} \int_{\Omega} r u dx dt. $$
The weak convergence of $\Pi_{\mathcal{D},\tau}u(\cdot, t_0)$ to $u(\cdot, t_0)$ and of $\nabla_{\mathcal{D},\tau}u$ to $\nabla u$ implies that

$$ \liminf_{m \to \infty} \frac{1}{2} \int_{\Omega} \Pi_{\mathcal{D},\tau} u(x, t_0)^2 \geq \frac{1}{2} \int_{\Omega} u(x, t_0)^2 dx, \\ \liminf_{m \to \infty} \int_{0}^{\nu(t_0)\tau} G_{\mathcal{D},\tau}(\Pi_{\mathcal{D},\tau} u, x, t) \nabla_{\mathcal{D},\tau} u \cdot \nabla_{\mathcal{D},\tau} u dt dx \geq \int_{0}^{t_0} \int_{\Omega} G(u, x, t) \nabla u \cdot \nabla u dx dt. $$

Therefore we obtain

$$ \lim_{m \to \infty} \frac{1}{2} \int_{\Omega} \Pi_{\mathcal{D},\tau} u(x, t_0)^2 = \frac{1}{2} \int_{\Omega} u(x, t_0)^2 dx, \\ \lim_{m \to \infty} \int_{0}^{\nu(t_0)\tau} G_{\mathcal{D},\tau}(\Pi_{\mathcal{D},\tau} u, x, t) \nabla_{\mathcal{D},\tau} u \cdot \nabla_{\mathcal{D},\tau} u dt dx = \int_{0}^{t_0} \int_{\Omega} G(u, x, t) \nabla u \cdot \nabla u dx dt. $$

hence concluding the proof of the strong convergences. □

## 4 A particular gradient scheme

In order to describe the scheme, we now introduce some notations for the space discretization.

<figure><img src="image_2.png" /><figcaption>Figure 1: Notations for the meshes</figcaption></figure>

1. A rectangular discretization of $\Omega$ is defined by the increasing sequences $a_i = x_0^{(i)} < x_1^{(i)} < \dots < x_{n(i)}^{(i)} = b_i$, $i = 1, \dots, d$.

2. We denote by

$$ \mathcal{M} = \{ ]x_{i(1)}^{(1)}, x_{i(1)+1}^{(1)}[ \times \dots \times ]x_{i(d)}^{(d)}, x_{i(d)+1}^{(d)}, \dots, 0 \le i^{(1)} < n^{(1)}, \dots, 0 \le i^{(d)} < n^{(d)} \} $$

the set of the control volumes. The elements of $\mathcal{M}$ are denoted $p, q, \dots$. We denote by $x_p$ the center of $p$. For any $p \in \mathcal{M}$, let $\partial p = \bar{p} \setminus 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_{\mathcal{D}}$ denote the maximum value of $(h_p)_{p \in \mathcal{M}}$.

3. We denote by $\mathcal{E}_p$ the set of all the faces of $p \in \mathcal{M}$, by $\mathcal{E}$ the union of all $\mathcal{E}_p$, and for all $\sigma \in \mathcal{E}$, we denote by $|\sigma|$ its $(d-1)$-dimensional measure. For any $\sigma \in \mathcal{E}$, we define the set $\mathcal{M}_{\sigma} = \{p \in \mathcal{M}, \sigma \in \mathcal{E}_p\}$ (which has therefore one or two elements), we denote by $\mathcal{E}_p$ the set of the faces of $p \in \mathcal{M}$ (it has $2d$ elements) and by $x_{\sigma}$ the cent-re of $\sigma$. We then denote by $d_{p\sigma} = |x_{\sigma} - x_p|$ the orthogonal distance between $x_p$ and $\sigma \in \mathcal{E}_p$ and by $n_{p,\sigma}$ the normal vector to $\sigma$, outward to $p$.
4. We denote by $\mathcal{V}_p$ the set of all the vertices of $p \in \mathcal{M}$ (it has $2^d$ elements), by $\mathcal{V}$ the union of all $\mathcal{V}_p$, $p \in \mathcal{M}$. For $y \in \mathcal{V}_p$, we denote by $K_{p,y}$ the rectangle whose faces are parallel to those of $p$, and whose the set of vertices contains $x_p$ and $y$. We denote by $\mathcal{V}_\sigma$ the set of all vertices of $\sigma \in \mathcal{E}$ (it has $2^{d-1}$ elements), and by $\mathcal{E}_{p,y}$ the set of all $\sigma \in \mathcal{E}_p$ such that $y \in \mathcal{V}_\sigma$ (it has $d$ elements).

5. We define the set $X_D$ of all $u = ((u_p)_{p \in \mathcal{M}}, (u_{\sigma,y})_{\sigma \in \mathcal{E}, y \in \mathcal{V}_\sigma})$.

6. We denote, for all $u \in H_D$, by $\Pi_D u \in L^2(\Omega)$ the function defined by the constant value $u_p$ a.e. in $p \in \mathcal{M}$.

7. For $u \in X_D$, $p \in \mathcal{M}$ and $y \in \mathcal{V}_p$, we denote by

$$ \nabla_{p,y} u = \frac{2}{|p|} \sum_{\sigma \in \mathcal{E}_{p,y}} |\sigma|(u_{\sigma,y} - u_p)\mathbf{n}_{p,\sigma} = \sum_{\sigma \in \mathcal{E}_{p,y}} \frac{u_{\sigma,y} - u_p}{d_{p\sigma}}\mathbf{n}_{p,\sigma}, \quad (31) $$

and by $\nabla_D u$ the function defined a.e. on $\Omega$ by $\nabla_{p,y} u$ on $K_{p,y}$.

We then have the following result.

**Lemma 4.1** Let $\Omega =]a_1, b_1[ \times \dots \times ]a_d, b_d[$ be an open rectangle in $\mathbb{R}^d$. Let $\mathcal{D} = (X_D, \Pi_D, \nabla_D)$ be a rectangular discretization as described above. Then $\mathcal{D}$ is an approximate gradient discretization in the sense of Definition 3.1, such that, under the regularity condition $(x_{i(j)}^{(j)} - x_{i(j)})/(x_{i(k)}^{(k)} - x_{i(k)}) \le C$, we get the existence of $T: \mathbb{R}^d \to \mathbb{R}^+$ with $\lim_{|\xi| \to 0} T(\xi) = 0$ such that $T_D \le T$ independently of $h_D$, and that

$$ \forall \varphi \in H^1(\Omega), \quad \lim_{h_D \to 0} S_D(\varphi) = 0, \quad (32) $$

and

$$ \forall \varphi \in H_{\mathrm{div},0}(\Omega), \quad \lim_{h_D \to 0} W_D(\varphi) = 0. \quad (33) $$

**PROOF.** Let us recall the result, proved in [11]: for such a rectangular discretization $\mathcal{D}$, the expression $\|u\|_{\mathcal{D}}$, defined by

$$ \|u\|_{\mathcal{D}}^2 = \sum_{p \in \mathcal{M}} |p| u_p^2 + \sum_{p \in \mathcal{M}} \sum_{\sigma \in \mathcal{E}_p} \sum_{y \in \mathcal{V}_\sigma} \frac{|\sigma|}{d_{p\sigma}} (u_{\sigma,y} - u_p)^2, \quad \forall u \in X_D, $$

is a norm on $X_D$ such that (32) holds, where $C$ only depends on the bound on $\theta$. The limit conformity property (33) is proved in the same way as in [11]. We then remark that $\|u\|_{\mathcal{D}}^2$ controls the semi-norm $|\cdot|_{1,\tau}$ as defined in [10, Definition 10.2, p. 795]. We can therefore follow the proof of [10, Theorem 10.3, p. 810] using the discrete trace inequality [10, Lemma 10.5, p. 807] in the case of the homogeneous Neumann boundary conditions, which proves the existence of $T$ as given in this statement (hence proceeding in the same way as in [9]). $\square$

**Remark 4.1** The equations obtained, for a given $y \in \mathcal{V}$, defining $v \in X_D$ for a given $\sigma \in \mathcal{E}_y$ by $v_{\sigma,y} = 1$ and all other degrees of freedom null, constitute a local invertible linear system, allowing for expressing all $(u_{\sigma,y})_{\sigma \in \mathcal{E}_y}$ with respect to all $(u_{p,y})_{p \in \mathcal{M}}$. This leads to a nine-point stencil on rectangular meshes in $2D$, $27$-point stencil in $3D$ (this property is the basis of the MPFA O-scheme [1]).

# 5 Numerical experiments

## Example 1.

This example is devoted to a 1D illustration of the time-delayed regularization of the Perona-Malik equation, as described by (13) in the introduction of this paper. We consider the case $\Omega =]0, 1[$, $u_{\mathrm{ini}}(x) = 2x$, and $\alpha = 1/101$ in (13). We then apply the scheme presented in Section 4 (it then resumes to a
<figure><img src="image_3.png" /><figcaption>Figure 2: $\bar{t} = .2$, $\tau = .01$ (left), $\bar{t} = .5$, $\tau = .01$ (right)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: $\bar{t} = 1$, $\tau = .01$ (left), $\bar{t} = 5$, $\tau = .1$ (right)</figcaption></figure>

standard 3-point finite volume scheme), using various values of $\bar{t}$. We then show in Figure 2 the results computed at the final time $T = 50$, for $h = 1/500$ for $\bar{t} = 0.2$ and $\bar{t} = 0.5$ and $\tau = .01$. We see that, although the problem is regularized, the numerical solution shows discontinuities resulting from the ill-posedness of the Perona-Malik equation. Turning to the values $\bar{t} = 1$ and $\bar{t} = 5$, with $\tau = .1$, we observe in Figure 3 that the regularity of the numerical solution is significantly increased.

**Example 2.** In this example we present results of numerical computations for the regularized Perona-Malik equation in $2D$, where the regularization is again of the type (13). We assume that $\Omega = [0, 1] \times [0, 1]$, and we chose the right hand side $r$ such that $u(x, y, t) = C((x^2 + y^2)/2 - (x^3 + y^3)/3)t$ is an exact strong solution of the problem (1)-(2)-(3). Two test cases, with different time delays $\bar{t}$, final time $T$ and coefficient $C$, are presented. The errors in $L_2((0, T), L_2(\Omega))$ and $L_\infty((0, T), L_2(\Omega))$ for the solution $u$ (denoted by $E_2, E_\infty$) and for its gradient (denoted by $EG_2, EG_\infty$) together with the experimental order of convergence ($EOC$) are presented in Table 1 in the case $\bar{t} = 0.625$, $T = 0.625$ and $C = 1$, and in Table 2 in the case $\bar{t} = 0.625$, $T = 6.25$ and $C = \frac{\bar{t}}{T} = 0.1$. These numerical results indicate that the method is of second order accurate, both in solution and gradient approximation.

**Example 3.** Image filtering by the time-delayed Perona-Malik model. In this experiment, we show the results obtained when filtering a noisy image, using the time-delay regularization of the Perona-Malik equation as given by (13). The original image (Figure 4 left top) is damaged by 40 percent
<table><caption>Table 1: Example 2, $\bar{t} = 0.0625$, $T = 0.625$, $C = 1$.</caption><thead><tr><th>n</th><th>τ</th><th>E<sub>2</sub></th><th>EOC</th><th>E<sub>∞</sub></th><th>EOC</th><th>EG<sub>2</sub></th><th>EOC</th><th>EG<sub>∞</sub></th><th>EOC</th></tr></thead><tbody><tr><td>4</td><td>0.0625</td><td>4.771e-4</td><td>-</td><td>1.022e-3</td><td>-</td><td>7.184e-3</td><td>-</td><td>1.450e-2</td><td>-</td></tr><tr><td>8</td><td>0.015625</td><td>1.172e-4</td><td>1.429</td><td>2.692e-4</td><td>1.925</td><td>1.707e-3</td><td>2.073</td><td>3.615e-3</td><td>2.004</td></tr><tr><td>16</td><td>0.00390625</td><td>2.913e-5</td><td>2.604</td><td>6.812e-5</td><td>1.982</td><td>4.213e-4</td><td>2.019</td><td>9.031e-4</td><td>2.001</td></tr><tr><td>32</td><td>0.0009765625</td><td>7.270e-6</td><td>2.002</td><td>1.708e-5</td><td>1.996</td><td>1.050e-4</td><td>2.004</td><td>2.257e-4</td><td>2.000</td></tr><tr><td>64</td><td>0.000244140625</td><td>1.815e-6</td><td>2.001</td><td>4.273e-6</td><td>1.999</td><td>2.624e-5</td><td>2.000</td><td>5.643e-5</td><td>1.999</td></tr></tbody></table>

<table><caption>Table 2: Example 2, $\bar{t} = 0.625$, $T = 6.25$, $C = \frac{\bar{t}}{T} = 0.1$</caption><thead><tr><th>n</th><th>τ</th><th>E<sub>2</sub></th><th>EOC</th><th>E<sub>∞</sub></th><th>EOC</th><th>EG<sub>2</sub></th><th>EOC</th><th>EG<sub>∞</sub></th><th>EOC</th></tr></thead><tbody><tr><td>4</td><td>0.0625</td><td>1.778e-3</td><td>-</td><td>1.209e-3</td><td>-</td><td>2.032e-2</td><td>-</td><td>1.400e-2</td><td>-</td></tr><tr><td>8</td><td>0.015625</td><td>4.594e-4</td><td>1.952</td><td>3.173e-4</td><td>1.930</td><td>5.034e-3</td><td>2.031</td><td>3.478e-3</td><td>2.001</td></tr><tr><td>16</td><td>0.00390625</td><td>1.158e-4</td><td>1.988</td><td>8.021e-5</td><td>1.984</td><td>1.255e-3</td><td>2.004</td><td>8.680e-4</td><td>2.002</td></tr><tr><td>32</td><td>0.0009765625</td><td>2.898e-5</td><td>1.998</td><td>2.011e-5</td><td>1.996</td><td>3.137e-4</td><td>2.000</td><td>2.169e-4</td><td>2.001</td></tr><tr><td>64</td><td>0.000244140625</td><td>7.248e-6</td><td>1.999</td><td>5.030e-6</td><td>1.999</td><td>7.842e-5</td><td>2.000</td><td>5.422e-5</td><td>2.000</td></tr></tbody></table>

additive noise (Figure 4 right top). In the bottom row of Figure 4 we present the 1st, 10th and 25th denoising steps which show a correct reconstruction of the original figure. The following parameters are used in the computations: $n^{(1)} = n^{(2)} = 200$, $h = 0.1$, $\tau = 0.1$, $\bar{t} = 0.1$. In the nonlinear function $g$ defined by (6), we use the value of parameter $K = 100$. We observe that the denoised image does not show any significant alteration after 25 steps, compared to the initial one.

**Example 4.** Image filtering by the nonlinear tensor diffusion. In this example we present the image denoising by the so-called coherence enhancing tensor diffusion [18], that means the case, where the diffusion coefficient is given by a matrix. The mathematical model is given in the introduction of this paper by (9)-(10) and (11)-(12). In this experiment we use the following parameters: $n^{(1)} = n^{(2)} = 250$, $h = 0.001$, $\bar{t} = 0.000001$, $\tau = 0.000001$. The chosen parameters in the convolution operators are defined by $\sigma = 0.0001$ and $\rho = 0.01$, and the parameters in (11) are $\alpha = 0.0001$, $C = 1$. The scheme uses explicit values for $G$, which leads to define a semi-implicit scheme for the nonlinear tensor diffusion model as in [9, 8]. In Figure 5, one can see the original image with three crackling lines, and the results after 2 and 5 time steps respectively, showing that the coherence of the line structures has been improved.

# 6 Conclusion

Gradient schemes have recently been shown to be efficient for the discretization of elliptic problems and steady state coupled problems. The present paper shows that the mathematical framework of the Gradient Scheme also applies for the discrete analysis of parabolic problems involved in some practical applications. The particular gradient scheme used in this paper for image processing applications leads to a 9-point stencil finite volume type scheme, which is unconditionally coercive and convergent. Its use in further applications will be studied in future works.

# Appendix

**Lemma 6.1** Let $T > 0$ be given, and let $\tau = T/N_T$, for $N_T \in \mathbb{N}^*$. For all $t \in [0, T]$, we denote by $\nu(t)$ the element $n \in \mathbb{N}$ such that $t \in ((n-1)\tau, n\tau])$. Let $(a^n)_{n=1,...,N_T}$ be a family of non negative real values. Then

$$ \int_0^{T-\eta} \sum_{n=\nu(t)}^{\nu(t+\eta)} (\tau a^n) dt = \eta \sum_{n=1}^{N_T} (\tau a^n), \quad \forall \eta \in (0, T), \quad (34) $$
<figure><img src="image_5.png" /><figcaption>Figure 4: Example 3, the original image (left top), the noisy image (right top) and the results after 1, 10, 25 time steps (bottom, from left to right) of filtering.</figcaption></figure>

and

$$ \int_{0}^{T-\eta} \left( \sum_{n=\nu(t)+1}^{\nu(t+\eta)} \tau \right) a^{\nu(t+s)} dt \le \eta \sum_{n=1}^{N_T} (\tau a^n), \quad \forall \eta \in (0, T), \forall s \in (0, \eta). \quad (35) $$

The proof of the above lemma follows that of [10, Proof of Lemma 18.6 p 855].

**Theorem 6.1** Let $\Omega$ be an open bounded subset of $\mathbb{R}^2$, $a < b \in \mathbb{R}$ and $(u_m)_{m \in \mathbb{N}}$ be a sequence of functions from $[a, b]$ to $L^2(\Omega)$, such that there exists $C_1 > 0$ with

$$ \|u_m(t)\|_{L^2(\Omega)} \le C_1, \quad \forall m \in \mathbb{N}, \forall t \in [a, b]. \quad (36) $$

We also assume that there exists $(\tau_m)_{m \in \mathbb{N}}$ with $\tau_m \ge 0$ and $\lim_{m \to \infty} \tau_m = 0$ and a dense subset $R$ of $L^2(\Omega)$ such that, for all $\varphi \in R$, there exists $C_\varphi > 0$ such that

$$ |\langle u_m(t_2) - u_m(t_1), \varphi \rangle_{L^2(\Omega), L^2(\Omega)}| \le C_\varphi (t_2 - t_1 + 2\tau_m)^{1/2}, \quad \forall m \in \mathbb{N}, \forall a \le t_1 \le t_2 \le b. \quad (37) $$

Then there exists $u \in L^\infty(a, b; L^2(\Omega))$ with $u \in C_w([a, b], L^2(\Omega))$ (where we denote by $C_w([a, b], L^2(\Omega))$ the set of functions from $[a, b]$ to $L^2(\Omega)$, continuous for the weak topology of $L^2(\Omega)$) and a subsequence of $(u_m)_{m \in \mathbb{N}}$, again denoted $(u_m)_{m \in \mathbb{N}}$, such that, for all $t \in [a, b]$, $u_m(t)$ converges to $u(t)$ for the weak topology of $L^2(\Omega)$.
<figure><img src="image_6.png" /><figcaption>Figure 5: Example 4, the original image (left), the results after 2 (middle) and 5 (right) time steps.</figcaption></figure>

PROOF. The proof follows that of Ascoli's theorem, and we provide it only for completeness (a similar proof is also provided in [15] in the framework of the weak star topology of Radon measures). Let $(t_p)_{p \in \mathbb{N}}$ be a sequence of real numbers, dense in $[a, b]$. Due to (36), for each $p \in \mathbb{N}$, we may extract from $(u_m(t_p))_{m \in \mathbb{N}}$ a subsequence which is convergent to some element of $L^2(\Omega)$ for the weak topology of $L^2(\Omega)$. Using a diagonal method, we can choose a sub-sequence, again denoted $(u_m)_{m \in \mathbb{N}}$, such that $(u_m(t_p))_{m \in \mathbb{N}}$ is weakly convergent for all $p \in \mathbb{N}$. For any $t \in [a, b]$ and $v \in L^2(\Omega)$, we then prove that the sequence $((u_m(t), v)_{L^2(\Omega), L^2(\Omega)})_{m \in \mathbb{N}}$ is a Cauchy sequence. Indeed, let $\varepsilon > 0$ be given. We first choose $\varphi \in R$ such that $\|\varphi - v\|_{L^2(\Omega)} \le \varepsilon$. Then, we choose $p \in \mathbb{N}$ such that $|t - t_p| \le (\varepsilon/C_\varphi)^2$. Since $((u_m(t_p), \varphi)_{L^2(\Omega), L^2(\Omega)})_{m \in \mathbb{N}}$ is a Cauchy sequence, we choose $n_0 \in \mathbb{N}$ such that, for $k, l \ge n_0$,

$$|\langle u_k(t_p) - u_l(t_p), \varphi \rangle_{L^2(\Omega), L^2(\Omega)}| \le \varepsilon,$$

and such that $\tau_k, \tau_l \le (\varepsilon/C_\varphi)^2$. We then get, using (37),

$$|\langle u_k(t) - u_l(t), \varphi \rangle_{L^2(\Omega), L^2(\Omega)}| \le C_\varphi ( (|t - t_p| + 2\tau_k)^{1/2} + (|t - t_p| + 2\tau_l)^{1/2}) + \varepsilon,$$

which gives

$$|\langle u_k(t) - u_l(t), \varphi \rangle_{L^2(\Omega), L^2(\Omega)}| \le 2 \ 3^{1/2} \varepsilon + \varepsilon.$$

We then get, using (36),

$$|\langle u_k(t) - u_l(t), v \rangle_{L^2(\Omega), L^2(\Omega)}| \le 2C_1\varepsilon + 2 \ 3^{1/2}\varepsilon + \varepsilon.$$

This proves that the sequence $((u_m(t), v)_{L^2(\Omega), L^2(\Omega)})_{m \in \mathbb{N}}$ converges. Since

$$|\langle u_m(t), v \rangle_{L^2(\Omega), L^2(\Omega)}| \le C_1 \|v\|_{L^2(\Omega)},$$

we get the existence of $u(t) \in L^2(\Omega)$ such that $(u_m(t))_{m \in \mathbb{N}}$ converges to $u(t)$ for the weak topology of $L^2(\Omega)$. Then $u \in C_w([a, b], L^2(\Omega))$ is obtained by passing to the limit $m \to \infty$ in (37), and by using the density of $R$ in $L^2(\Omega)$. $\square$

**Acknowledgement.** The work of A.H., K.M. and O.S. was supported by the grants APVV-0184-10 and VEGA 1/0269/09.

## References

[1] I. Aavatsmark, T. Barkve, O. Boe, and T. Mannseth. Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media. J. Comput. Phys., 127(1):2–14, 1996.
[2] H. Amann. Time-delayed Perona-Malik type problems. *AMUC*, LXXVI(1):15–38, 2007.

[3] S. Bartels and A. Prohl. Stable discretization of scalar and constrained vectorial Perona-Malik equation. *Interfaces and Free Boundaries. Mathematical Modelling, Analysis and Computation*, 9(4):431–453, 2007.

[4] G. Bellettini, M. Novaga, M. Paolini, and C. Tornese. Convergence of discrete schemes for the Perona-Malik equation. *Journal of Diff. Equations*, 245:892–924, 2008.

[5] Clément Cancès and Thierry Gallouët. On the time continuity of entropy solutions. *J. Evol. Equ.*, 11(1):43–55, 2011.

[6] F. Catté, P. L. Lions, J. M. Morel, and T. Coll. Image selective smoothing and edge detection by nonlinear diffusion. *SIAM J. Num. Anal.*, 29:182–193, 1992.

[7] A. Handlovičová and Z. Krivá. Error estimates for finite volume scheme for Perona-Malik equation. *AMUC*, LXXIV(1):79–94, 2005.

[8] O. Drblíková A. Handlovičová and K. Mikula. Error estimates of the finite volume scheme for the nonlinear tensor-driven anisotropic diffusion. *Applied Numerical Mathematics*, 59:2548–2570, 2009.

[9] O. Drblíková and K. Mikula. Convergence analysis of finite volume scheme for nonlinear tensor anisotropic diffusion in image processing. *SIAM Journal on Numerical Analysis*, 46(1):37–60, 2007.

[10] R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. In P. G. Ciarlet and J.-L. Lions, editors, *Techniques of Scientific Computing, Part III*, Handbook of Numerical Analysis, VII, pages 713–1020. North-Holland, Amsterdam, 2000.

[11] R. Eymard, T. Gallouët, and R. Herbin. Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes SUSHI: a scheme using stabilization and hybrid interfaces. *IMA J. Numer. Anal.*, 30(4):1009–1043, 2010.

[12] R. Eymard, C. Guichard, and R. Herbin. Small-stencil 3d schemes for diffusive flows in porous media. *M2AN Math. Model. Numer. Anal.*, 46:265–290, 2012.

[13] R. Eymard and R. Herbin. Gradient scheme approximations for diffusion problems. In *Finite volumes for complex applications VI*, volume 1, pages 439–447. Springer, London, 2011. Finite Volumes for Complex Applications VI (FVCA VI), Prague, Czech Republic, June 2011.

[14] R. Eymard, Herbin R., and J.-C. Latché. Convergence analysis of a colocated finite volume scheme for the incompressible Navier-Stokes equations on general 2 or 3d meshes. *SIAM J. Numer. Anal.*, 45(1):1–36, 2007.

[15] Robert Eymard, Sophie Mercier, and Alain Prignet. An implicit finite volume scheme for a scalar hyperbolic problem with measure data related to piecewise deterministic Markov processes. *J. Comput. Appl. Math.*, 222(2):293–323, 2008.

[16] K. Mikula and N. Ramarosy. Semi-implicit finite volume scheme for solving nonlinear diffusion equations in image processing. *Numer. Math.*, 89(3):561–590, 2001.

[17] P. Perona and J. Malik. Scale space an edge detection using anisotropic diffusion. *IEEE transactions on pattern Analysis and machine Intelligence*, 12(7):629–639, 1990.

[18] J. Weickert. Coherence-enhancing diffusion filtering. *Int. J. Comput. Vision*, 31:111–127, 1999.