A class of elliptic systems with discontinuous variable exponents and L1 data for image denoising

Nonlinear Analysis: Real World Applications 50 (2019) 448–468

Contents lists available at ScienceDirect

Nonlinear Analysis: Real World Applications www.elsevier.com/locate/nonrwa

A class of elliptic systems with discontinuous variable exponents and L1 data for image denoising Dazhi Zhang a , Kehan Shi b , Zhichang Guo a ,∗, Boying Wu a a b

Department of Mathematics, Harbin Institute of Technology, Harbin, 150001, China Department of Mathematics, China Jiliang University, Hanzhou, 310018, China

article

info

abstract

Article history: Received 22 March 2018 Received in revised form 26 May 2019 Accepted 26 May 2019 Available online 3 June 2019 Keywords: Elliptic system Discontinuous variable exponent L1 data Image denoising

This paper investigates a class of elliptic systems consisting of the p(x)-Laplacian equation and the Poisson equation for image denoising. Under the assumption }, where p− := ess inf x∈Ω p(x) and N is the dimension of Ω , that p− > max{1, N 3 we prove the existence and uniqueness of weak solutions for the homogeneous Neumann boundary value problem with discontinuous variable exponent p(x) and L1 data. The proof, which is based on Schauder’s fixed-point theorem, also provides an iterative scheme for solving the system numerically. Experimental results illustrate that the proposed system with piecewise constant p(x) performs better than commonly used smooth p(x) in image denoising. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction In this paper, we investigate the existence and uniqueness of weak solutions for the coupled elliptic system ( ) p(x)−2 div g(x)|∇u| ∇u = λv, in Ω , (1) ∆v = f − u, p(x)−2

g(x)|∇u|

∇u · ⃗n = ∇v · ⃗n = 0,

in Ω ,

on ∂Ω ,

(2) (3)

and its application in the reconstruction of images corrupted by Gaussian noise. Here Ω is a bounded domain of RN with Lipschitz boundary ∂Ω , ⃗n is the outward pointing unit normal vector to ∂Ω , N ≥ 2 and λ is a positive free parameter. For a given noisy image f (x), the solution u(x) is the restored image. System (1)–(3) can be regarded as an approximation of the H −1 model [1] for image processing. The authors of [1] proposed to remove Gaussian noise in a grayscale image f by minimizing the energy ∫ λ (4) |∇u|dx + ∥u − f ∥2H −1 (Ω) . 2 Ω ∗

Corresponding author. E-mail addresses: [email protected], [email protected] (Z. Guo).

https://doi.org/10.1016/j.nonrwa.2019.05.012 1468-1218/© 2019 Elsevier Ltd. All rights reserved.

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

449

∫ 2 The H −1 norm is defined as ∥ · ∥2H −1 = Ω |∇∆−1 (·)| dx, where the inverse Laplacian operator ∆−1 is ∫ understood in the following sense: If u − f ∈ L2 (Ω ) and Ω (u − f )dx = 0, then by the classical results (see, for instance, [2]), problem { ∆P = u − f, in Ω , ∇P · ⃗n = 0, on ∂Ω , ∫ admits a unique solution P ∈ H 1 (Ω ) with Ω P dx = 0. We hence define P := ∆−1 (u − f ). ∫ The total variation (TV) regularization Ω |∇u|dx, which was first introduced by Rudin et al. [3], is one of the most popular and successful methodologies for image restoration. It has been studied extensively and has been proven to be an invaluable tool for preserving edges in image restoration. However, it also causes the staircasing effect, where areas of gradual intensity transitions are approximated by piecewise constant regions separated by sharp edges. Blocky and non-natural structures (like a staircase in 1D) can be observed. To overcome this drawback, regularization with variable exponents [4–6] has been proposed to replace the TV regularization. In this case, energy (4) becomes ∫ λ g(x) p(x) |∇u| dx + ∥u − f ∥2H −1 (Ω) , (5) 2 Ω p(x) where g(x) and p(x) are usually taken as g(x) =

1

2,

1 + k1 |(∇Gσ1 ∗ f )(x)|

p(x) = 1 +

1

2,

1 + k2 |(∇Gσ2 ∗ f )(x)|

(6)

with k1 , k2 , σ1 , σ2 being positive constants and Gσ being the Gaussian kernel. By Euler–Lagrange’s formula, minimizing energy (5) is equivalent to solve the fourth order elliptic equation ( ) p(x)−2 div g(x)|∇u| ∇u = λ∆−1 (f − u), (7) p(x)−2

with Neumann boundary conditions g(x)|∇u| ∇u · ⃗n = 0 and ∇(∆−1 (f − u)) · ⃗n = 0 on ∂Ω . By introducing a new function v := ∆−1 (f − u) in (7), we split the fourth order elliptic equation into a second order elliptic system, which is the proposed system (1)–(3). In this paper, we study system (1)–(3) with f (x) ∈ L1 (Ω ). Generally, it is reasonable to study solutions of image processing problems in the space of functions of bounded variation, BV (Ω ), allowing for discontinuities which are necessary for edge reconstruction. However, Chan et al. [7,8] showed that geometrical features are better preserved in L1 (Ω ) than in BV (Ω ) in image decomposition. Besides, the study of elliptic systems with L1 data is more challenging from a mathematical perspective. We also allow the presence of discontinuities in the measurable variable exponent p(x). The reason for choosing p(x) in the way as in (6) is not well explained by researchers. By dropping the continuity assumption for p(x), we are able to construct more general variable exponents for image denoising. For instance, the restored image achieves better quality if p(x) is constructed as an appropriate piecewise constant function. Both aspects distinguish the proposed system (1)–(3) from related works concerning the approximation of the H −1 model [6,9–11]. Since f (x) is merely in L1 (Ω ), it seems reasonable to work with entropy solutions [12,13] or renormalized solutions [14–16] for system (1)–(3), which need less regularity than the usual weak solutions. The notion of renormalized solutions was first introduced by DiPerna and Lions [14] in the study of the Boltzmann equation, and the notion of entropy solutions was first introduced by B´enilan et al. [12] in the study of p-Laplacian type elliptic equations. Notice that the right-hand sides of Eq. (1) and Eq. (2) are coupled to each other, we could not define entropy solutions or renormalized solutions for the system in a similar way to single equations. As a remedy, we shall define weak solutions for system (1)–(3) under the restriction p− > max{1, N3 }, where p− := ess inf x∈Ω p(x). This is feasible since f (x) ∈ L1 (Ω ) appears only on the right-hand side of Eq. (2).

450

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

The variable exponent p(x) in Eq. (1) is allowed to be discontinuous, which causes some difficulties. In ∫ p(x) this case, the integrand of the functional Ω |∇u| dx is not regular [17]. As a consequence, the undesired Lavrentiev phenomenon occurs. Besides, Poincar´e’s inequality and Sobolev’s imbedding theorem are no longer valid for functions in the standard Sobolev space with variable exponents. We shall define weak solutions for Eq. (1) in a non-standard Sobolev space with variable exponents. At last, the difficulty caused by the coupling of system (1)–(3) is solved by Schauder’s fixed-point theorem. That is to say, we propose to replace v on the right-hand side of Eq. (1) by a new function w and decouple the system into two subproblems, which are the classical p(x)-Laplacian equation and the Poisson equation. Then we construct a mapping from w to v and prove that it has a fixed point. The rest of this paper is organized as follows. In Section 2, we recall the definitions of function spaces with variable exponents, define function spaces for our problem and present the main result of this paper. In Section 3, we decouple system (1)–(3) into two subproblems and discuss both of them. The existence and uniqueness of weak solutions for the system is proven in Section 4. Numerical scheme and experiments are given in Section 5 to explain the role of the variable exponent p(x) and to illustrate the image denoising performance of system (1)–(3). 2. Preliminaries and the main result We assume that p(x) and g(x) are measurable functions such that 1 < p− ≤ p(x) ≤ p+ ≤ 2 and 0 < g− ≤ g(x) ≤ g+ ≤ 1, where p− := ess inf p(x), x∈Ω

p+ := ess sup p(x). x∈Ω

In the following, C is a generic constant that takes different values at different positions. We also use the p(x) p(x) and p∗ (x) = NN−p(x) . notations p′ (x) = p(x)−1 p 1,p By L (Ω ) and W (Ω ) we mean the classical Lebesgue space and Sobolev space. H 1 (Ω ) represents the Hilbert space W 1,2 (Ω ). We now give definitions of some function spaces involving variable exponents. Let us start with a brief overview of the generalized Lebesgue spaces Lp(x) (Ω ). A detailed description can be found in [18–20]. By Lp(x) (Ω ) we denote the space consisting of measurable functions f (x) in Ω such that ∫ p(x) ϱp(x) (f ) := |f (x)| dx < ∞. Ω p(x)

The space L

(Ω ) equipped with the norm ∥f ∥p(x),Ω = ∥f ∥p(·) := inf{α > 0 : ϱp(x) (f /α) ≤ 1},

becomes a Banach space. The following two properties of Lp(x) (Ω ) will be used throughout this paper. • It follows directly from the definition that { } { } p+ p+ p− p− min ∥f ∥p(·) , ∥f ∥p(·) ≤ ϱp(x) (f ) ≤ max ∥f ∥p(·) , ∥f ∥p(·) . { min ϱp(x) (f )

1 p−

, ϱp(x) (f )

1 p+

}

{ ≤ ∥f ∥p(·) ≤ max ϱp(x) (f )

• We have the H¨ older inequality ∫ |f (x)g(x)|dx ≤ 2∥f ∥p(·) ∥g∥p′ (·) , Ω ′

for any f (x) ∈ Lp(x) (Ω ) and g(x) ∈ Lp (x) (Ω ).

1 p−

, ϱp(x) (f )

(8) 1 p+

} .

(9)

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

451

Since the variable exponent p(x) may be discontinuous, we define the generalized Sobolev space W 1,p(x) (Ω ) by W 1,p(x) (Ω ) := {u ∈ W 1,p− (Ω ) :

∫

p(x)

|∇u(x)|

dx < ∞},

Ω

equipped with the norm ∥u∥W 1,p(x) (Ω) = ∥u∥1,p(·) := ∥u∥1,p− + ∥∇u∥p(·) . It is easy to check that W 1,p(x) (Ω ) is a separable and reflexive Banach space. At last, since the Neumann 1,p(x) boundary value problem is considered, let W0 (Ω ) be the subset of W 1,p(x) (Ω ) such that ∫ 1,p(x) W0 (Ω ) := {u ∈ W 1,p(x) (Ω ) : udx = 0}. Ω

W01,q (Ω )

The space with a constant q is defined in a similar way. We would like to emphasize that Poincar´e’s inequality and the embedding inequalities fail for functions 1,p(x) in the space W0 (Ω ) if p(x) is discontinuous. But we still have the Poincar´e inequality ⎫ ⎧ (∫ ) pp− ⎬ ∫ ∫ ⎨∫ + p(x) p(x) p p , |∇u| dx, |∇u| dx |∇u| − dx ≤ C max |u| − dx ≤ C ⎭ ⎩ Ω Ω Ω Ω 1,p(x)

for u ∈ W0 (Ω ) due to W 1,p− (Ω ) ⊂ W 1,p(x) (Ω ), Holder’s inequality and (9). In the following, the concept of entropy solutions will be utilized. For this, we need the notion of truncated functions. The truncated function Tk (x) at height k ≥ 0 is taken as ⎧ if x > k, ⎨ k, x, if |x| ≤ k, Tk (x) := min{k, max{x, −k}} = ⎩ −k, if x < −k. It can be observed that |Tk (u)| ≤ k, Tk (−u) = −Tk (u) and ∇Tk (u) = χ{|u|0

for ∀ z ∈ Rn with |z| = R :

n ∑

fi (z)zi ≥ 0.

i=1

Then there exists z ∗ ∈ Rn , such that fi (z ∗ ) = 0 for all i = 1, . . . , n and |z ∗ | ≤ R. We now give the definition of weak solutions for system (1)–(3). Definition 2. A pair of measurable functions (u, v) is called a weak solution of system (1)–(3), if u ∈ W 1,p(x) (Ω ), v ∈ W 1,q (Ω ), 1 ≤ q < NN−1 , Tk (v) ∈ H 1 (Ω ), ∀k > 0, such that ∫ ∫ p(x)−2 g(x)|∇u| ∇u · ∇ϕdx + λ vϕdx = 0, (10) Ω

∫

Ω

∫ ∇v · ∇Tk (v − φ)dx +

Ω

(f − u)Tk (v − φ)dx ≤ 0, Ω

hold for any ϕ ∈ W 1,p(x) (Ω ), φ ∈ H 1 (Ω ) ∩ L∞ (Ω ) and k > 0.

(11)

452

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

Remark 1. If p− > max{1, N3 }, then every term in Definition 2 is well defined. In fact, by the embedding ∗ ∗ q W 1,q (Ω ) ↪→ Lq (Ω ), where q ∗ = NN−q ∈ [1, NN−2 ), we have v ∈ Lq (Ω ). If p− > max{1, N3 }, again by the ∗ embedding inequality, ϕ ∈ W 1,p(x) (Ω ) ⊂ W 1,p− (Ω ) ⊂ Lp− (Ω ), where p∗− ∈ [1, ∞)\( N2 , ∞). This implies that the second term of Eq. (10) is well defined. It is obvious that the rest terms in Eqs. (10) and (11) are well defined for f ∈ L1 (Ω ). Remark 2. It follows from (11) that ∫

∫ ∇v · ∇φdx +

Ω

(f − u)φdx = 0,

(12)

Ω

′

for any φ ∈ W 1,q (Ω ) ∩ L∞ (Ω ). The proof can be found in [22, Theorem 6.5]. Our main result is as follows. Theorem 3. Let Ω be a bounded domain of RN with Lipschitz boundary ∂Ω , N ≥ 2 and λ > 0 be sufficiently small. If p− > max{1, N3 }, then for any f ∈ L1 (Ω ), system (1)–(3) admits a unique weak solution (u, v). ∫ ∫ Remark 3. Let φ = 1 in (12), we have Ω udx = Ω f dx. It implies that the residual f − u satisfies the property of noise or texture of being of zero mean. Remark 4. In image (or video) processing, where N = 2 (or 3), Theorem 3 holds true for any p(x) with p− > 1. Consequently, we can choose p(x) close to 1+ at edges of images, so that system (1)–(3) would not blur images. 3. The decoupled subproblems In this section, we decouple system (1)–(3) into two independent subproblems, which are the p(x)Laplacian equation and the Poisson equation. This is achieved by replacing v in Eq. (1) by a new function w. The a priori estimates, as well as the existence and uniqueness of weak solutions are established for both subproblems. 3.1. The p(x)-Laplacian equation Let us consider the first subproblem div(g(x)|∇u|

p(x)−2

g(x)|∇u|

∇u) − λw = 0, in Ω ,

p(x)−2

∇u · ⃗n = 0, on ∂Ω .

(13) (14)

We have the following result. Theorem 4. Let w ∈ W01,q (Ω ), 1 ≤ q < NN−1 , p− > max{1, N3 }. Then for any λ > 0, problem (13)–(14) 1,p(x) admits a unique weak solution u ∈ W0 (Ω ), such that ∫ ∫ p(x)−2 g(x)|∇u| ∇u · ∇ϕdx + λ wϕdx = 0, (15) Ω

holds for any ϕ ∈ W

1,p(x)

Ω

(Ω ). Moreover, ∥u∥1,p(·) ≤ C∥w∥s(p∗ )′ , −

where s > 0 depends only on p− and p+ .

(16)

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

453

Proof . It follows from Remark 1 that the second term of (15) is well defined. In the case of Dirichlet boundary conditions, the existence and uniqueness of weak solutions for Eq. (13) is a straightforward generalization of the results obtained by Fan and Zhang [23]. Since the Neumann boundary condition is considered here and p(x) is allowed to be discontinuous, we give a complete proof below by the standard Galerkin method. 1,p(x) 1,p(x) (Ω ), 1. The space W0 (Ω ) is separable. Then there exists a sequence of functions {θi }∞ i=1 ⊂ W0 1,p(x) which is linearly independent and is dense in W0 (Ω ). Next, we introduce the Galerkin approximation of problem (13)–(14). Namely, for each n ∈ N, we look for un in the form un (x) =

n ∑

ani θi (x),

i=1

(an1 , an2 , . . . , ann )

where the constant coefficient ⃗a = is obtained by solving the following system of algebraic equations ∫ p(x)−2 Fi (⃗a) := g(x)|∇un | ∇un · ∇θi + λwθi dx = 0, i = 1, 2, . . . , n. (17) Ω ∫ ∫ It follows directly from Ω θi dx = 0 that Ω un dx = 0. 2. It is obvious that ∥⃗a∥ := ∥∇un ∥p(·) is an equivalent norm for the vector ⃗a ∈ Rn . By H¨older’s inequality, (8) and (9), ∫ p(x)−2 g(x)|∇un | ∇un · ∇θi dx ≤ 2∥|∇un |p(x)−2 ∇un ∥p′ (·) · ∥∇θi ∥p(·) Ω

≤ C∥∇un ∥sp(·) · ∥∇θi ∥p(·) = C∥⃗a∥s · ∥∇θi ∥p(·) , which means that Fi (⃗a) is a continuous function of ⃗a. Here s > 0 is a constant depending only on p− and p+ . ∫ By the fact that Ω un dx = 0, H¨ older’s inequality and (9), we have } { 1 1 p+ p− . ∥un ∥p∗− ≤ C∥∇un ∥p− ≤ C max ϱp(x) (∇un ) , ϱp(x) (∇un ) Multiplying (17) by ani and summing over i = 1, 2, . . . , n, we obtain ∫ ∫ n ∑ p(x)−2 Fi (⃗a)ani = g(x)|∇un | ∇un · ∇un dx + λ wun dx Ω

i=1

∫ ≥ g−

Ω

p(x)

|∇un |

dx − λ∥w∥(p∗− )′ ∥un ∥p∗− { } 1 1 ≥ g− ϱp(x) (∇un ) − C∥w∥(p∗− )′ · max ϱp(x) (∇un ) p− , ϱp(x) (∇un ) p+ ( { }) 1 1 p+ −1 p− −1 = ϱp(x) (∇un ) g− − C∥w∥(p∗− )′ max ϱp(x) (∇un ) , ϱp(x) (∇un ) , Ω

For any ⃗a with ∥⃗a∥ := ∥∇un ∥p(·) = R > 0, by utilizing (8), ∫ p(x) ϱp(x) (∇un ) = |∇un | dx ≥ min{Rp− , Rp+ }. Ω

∑n Recalling the assumption that ∥w∥(p∗− )′ is finite, then i=1 Fi (⃗a)ani ≥ 0 for large R. By Lemma 1, we find a vector ⃗a ∈ Rn that solves (17). 3. To pass to the limit n → +∞, we shall establish some estimates. It follows from step 2 that { } 1 1 g− ϱp(x) (∇un ) ≤ C∥w∥(p∗− )′ max ϱp(x) (∇un ) p− , ϱp(x) (∇un ) p+ ,

454

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

which leads us to

} { 1− 1 1− 1 min ϱp(x) (∇un ) p+ , ϱp(x) (∇un ) p− ≤ C∥w∥(p∗− )′ .

Then by (8), { } p+ −1 p− −1 min ∥∇un ∥p(·) , ∥∇un ∥p(·) ≤ C∥w∥(p∗− )′ . Consequently, we finally get ∥un ∥1,p(·) ≤ C∥∇un ∥p− + ∥∇un ∥p(·) ≤ C∥w∥s(p∗ )′ , −

(18)

where s > 0 depends only on p− and p+ . 4. Now we pass to the limit n → +∞ in (17). Since W 1,p(·) (Ω ) is a reflexive space, there exists a ( )N ⃗ ∈ Lp′ (x) (Ω ) , such subsequence of {un } (still denoted by itself) and limit functions u ∈ W 1,p(x) (Ω ), A that un ⇀ u, weakly in W 1,p(x) (Ω ), ( )N p(x)−2 ⃗ weakly in Lp′ (x) (Ω ) . g(x)|∇un | ∇un ⇀ A, ∫ 1,p(x) We also have (16) and u ∈ W0 (Ω ), due to (18) and the fact that Ω un dx = 0. By these convergence 1,p(x) results and the density of {θi }∞ (Ω ), we let n → +∞ in (17) and obtain i=1 in W0 ∫ ∫ ⃗ · ∇ϕdx + λ A wϕdx = 0, (19) Ω

Ω

∫ 1,p(x) for any ϕ ∈ W0 (Ω ). Recalling that Ω wdx = 0, we further have (19) holds for any ϕ ∈ W 1,p(x) (Ω ). To finish the proof of the existence of solutions, we need only to show ∫ ∫ 1,p(x) p(x)−2 ⃗ · ∇ϕdx = A g(x)|∇u| ∇u · ∇ϕdx, ∀ϕ ∈ W0 (Ω ). Ω

Ω

p−2

It relies on the monotonicity of the operator |s| s. Namely, for any ξ, η ∈ RN , ( ) p−2 p−2 2 p p p−2 |ξ| ξ − |η| η (ξ − η) ≥ (p − 1)|ξ − η| (|ξ| + |η| ) p , 1 < p < 2. According to the monotonicity, we have ∫ p(x)−2 p(x)−2 (|∇un | ∇un − |∇ζ| ∇ζ)(∇un − ∇ζ)dx ≥ 0, Ω 1,p(x)

for any ζ ∈ W0

(Ω ). Let us replace θi in (17) by un − ζ to get ∫ ∫ p(x)−2 |∇un | ∇un · ∇(un − ζ)dx + λ w(un − ζ)dx = 0. Ω

Ω

Combining the above two inequalities, we obtain ∫ ∫ p(x)−2 |∇ζ| ∇ζ · ∇(un − ζ)dx + λ w(un − ζ)dx ≤ 0. Ω

Ω

Substituting this inequality into (19) with ϕ = (un − ζ) and letting n → +∞, we have ∫ ⃗ − |∇ζ|p(x)−2 ∇ζ) · ∇(u − ζ)dx ≥ 0. (A Ω

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

455

1,p(x)

Taking ζ = u − αϕ with α ≥ 0 and ϕ ∈ W0 (Ω ) and letting α → 0+ , we obtain ∫ ⃗ − |∇u|p(x)−2 ∇u) · ∇ϕdx ≥ 0 ∀ϕ ∈ W 1,p(x) (Ω ). (A 0 Ω

Taking α ≤ 0, we can obtain an opposite inequality. This finishes the proof of existence of solutions. 5. We are left to prove the uniqueness of solutions. Let u1 and u2 be two solutions of problem (13)–(14), ∫ ∫ which means that u1 and u2 satisfy equation (15) and Ω u1 dx = Ω u2 dx = 0. Then ∫ p(x)−2 p(x)−2 (|∇u1 | ∇u1 − |∇u2 | ∇u2 )∇ϕ = 0, ∀ϕ ∈ W 1,p(x) (Ω ). Ω

Let ϕ = (u1 − u2 ) and utilizing the monotonicity of the operator |s| ∫ 2 |∇u1 − ∇u2 | dx = 0.

p−2

s, we have

Ω

This together with

∫ Ω

u1 dx =

∫ Ω

u2 dx = 0 shows the uniqueness of solutions.

□

3.2. The Poisson equation Now let us consider the second subproblem ∆v − (f − u) = 0, in Ω ,

(20)

∇v · ⃗n = 0, on ∂Ω .

(21)

We have the following result. ∫ ∫ Theorem 5. Let f ∈ L1 (Ω ), u ∈ L1 (Ω ) and Ω udx = Ω f dx. Then problem (20)–(21) admits a unique entropy solution v, such that v ∈ W01,q (Ω ), 1 ≤ q < NN−1 , Tk (v) ∈ H 1 (Ω ), ∀k > 0, and ∫ ∫ ∇v · ∇Tk (v − φ)dx + (f − u)Tk (v − φ)dx ≤ 0, (22) Ω

Ω

for any φ ∈ H 1 (Ω ) ∩ L∞ (Ω ) and k > 0. Moreover, ∥v∥W 1,q (Ω) + ∥Tk (v)∥H 1 (Ω) ≤ C(∥f ∥L1 (Ω) + ∥u∥L1 (Ω) ), for any 1 ≤ q <

N N −1

(23)

and k > 0.

Proof . The proof is mainly based on the ideas developed in [22]. Although some of the arguments are not new, we present a self-contained proof for the sake of clarity and readability. 1. Let us consider the approximate problem of (20)–(21), i.e., ∆vn −(fn − un ) = 0, ∇vn · ⃗n = 0,

in Ω ,

(24)

on ∂Ω ,

(25)

where fn ∈ C0∞ (Ω ) and un ∈ C0∞ (Ω ) are the mollification of f and u such that ∫ ∫ fn → f in L1 (Ω ), un → u in L1 (Ω ), fn dx = un dx, Ω

∥fn ∥L1 (Ω) ≤ ∥f ∥L1 (Ω) ,

∥un ∥L1 (Ω) ≤ ∥u∥L1 (Ω) .

Ω

456

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

Then for any fixed n > 0, problem (24)–(25) admits a unique weak solution vn ∈ H01 (Ω ), such that ∫ ∫ ∇vn · ∇φdx + (fn − un )φdx = 0, ∀φ ∈ H 1 (Ω ). Ω

(26)

Ω

Clearly, φ = Tk (vn ) (k > 0) can be taken as the test function in (26), i.e., ∫ ∫ 2 |∇Tk (vn )| dx = ∇vn · ∇Tk (vn )dx Ω ∫Ω = (−fn + un )Tk (vn )dx ≤ k(∥f ∥L1 + ∥u∥L1 ). Ω

By Poincar´e’s inequality, we have ∫ ∫ 2 2 |∇Tk (vn )| dx ≤ Ck(∥f ∥L1 + ∥u∥L1 ), |Tk (vn ) − Tk (vn )| dx +

(27)

Ω

Ω

∫ 1 T (v )dx is the average value. where Tk (vn ) = |Ω| Ω k n ∫ ∫ 2. By Lebesgue’s dominated convergence theorem, {|vn |>k} vn dx → 0 as k → +∞. Since Ω vn dx = 0, there exists a k > 0, such that ⏐ ⏐ ⏐∫ ⏐∫ ⏐ k ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ ⏐ vn dx⏐ ≤ |Ω |. vn dx⏐ = ⏐ ⏐ ⏐ 4 ⏐ ⏐ {|vn |>k} ⏐ {|vn | k} ≤ 14 |Ω |. Then ⏐ ⏐∫ ⏐ ⏐∫ ∫ ⏐ ⏐ ⏐ 1 1 ⏐⏐ ⏐ ⏐ ⏐ = |Tk (vn )| = T (v )dx v dx + kdx ⏐ ⏐ k n n ⏐ |Ω | ⏐ ⏐ |Ω | ⏐ Ω {|vn |k} ) ( 1 k k k ≤ |Ω | + |Ω | = . |Ω | 4 4 2 It follows the above inequality that |k − k2 | ≤ |k − Tk (vn )|. Then by the classical embedding theorem and (27), ( ) 2N ∫ ( ) 2N 1 N −2 k − k/2 N −2 dx 2 k {|vn |>k} {|vn |>k} ( ) 2N ∫ ( ) 2N 1 N −2 Tk (vn ) − Tk (vn ) N −2 ≤ dx 2 k {|vn |>k} ( ) 2N ∫ 2N 1 N −2 ≤ |Tk (vn ) − Tk (vn )| N −2 dx 2k Ω ( ) 2N ( ) N N −2 N 1 ∥f ∥L1 + ∥u∥L1 N −2 N −2 ≤C (k(∥f ∥L1 + ∥u∥L1 )) =C . k k ∫

meas{|vn | > k} =

1dx =

Now for fixed µ > 0, we have {|∇vn | ≥ µ} ⊂ {|∇vn | ≥ µ, |vn | < k} ∪ {|vn | > k}. Since meas{|∇vn | ≥ µ, |vn | < k} ≤

⏐ ∫ ⏐ ⏐ ∇Tk (vn ) ⏐2 ⏐ ⏐ dx ≤ k(∥f ∥L1 + ∥u∥L1 ) , ⏐ ⏐ µ µ2 Ω

we have meas{|∇vn | ≥ µ} ≤

k(∥f ∥L1 + ∥u∥L1 ) +C µ2

(

∥f ∥L1 + ∥u∥L1 k

)

N N −2

.

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

N −2

457

1

Let k = µ N −1 (∥f ∥L1 + ∥u∥L1 ) N −1 in the above inequality, we have ( meas{|∇vn | ≥ µ} ≤ C Let 1 ≤ q <

N N −1

∥f ∥L1 + ∥u∥L1 µ

)

N N −1

,

be fixed and t > 0, then ∫ ∫ q |∇vn | dx =

∫ q q |∇vn | dx + |∇vn | dx {|∇vn |
Ω

t N

C(q − 1) (∥f ∥L1 + ∥u∥L1 ) N −1 . = tq |Ω | + N N t N −1 −q N −1 − q Choosing t = ∥f ∥L1 + ∥u∥L1 , we obtain ∫ q |∇vn | dx ≤ C(∥f ∥L1 + ∥u∥L1 )q .

(28)

Ω

Now by Poincar´e’s inequality, we conclude from (28) that {vn } is uniformly bounded in W 1,q (Ω ) for any 1 ≤ q < NN−1 . 3. Since problem (24)–(25) is linear, we can repeat step 2 working with vn −vm , and find that (28) becomes ∫ q |∇vn − ∇vm | dx ≤ C(∥fn − fm ∥L1 + ∥un − um ∥L1 )q , Ω

< NN−1 . 1,q

Since {fn } and {un } are Cauchy sequences in L1 (Ω ), it follows that {vn } is a Cauchy for any 1 ≤ q sequence in W (Ω ). Then there exist a subsequence of {vn } (still denoted by itself) and a limit function v ∈ W 1,q (Ω ), such that (23) holds and vn → v, strongly in W 1,q (Ω ), ∇vn → ∇v, a.e. in Ω . ∫ Besides, we have v ∈ W01,q (Ω ), since Ω vn dx = 0. 4. Now we show that v is an entropy solution of problem (20)–(21). According to (27), {Tk (vn )} is a bounded (uniformly in n) sequence in H 1 (Ω ). This together with the above convergence results imply that Tk (v) ∈ H 1 (Ω ) and Tk (vn ) (up to subsequences) converges to Tk (v) weakly in H 1 (Ω ). We now fix k > 0 and φ ∈ H 1 (Ω ) ∩ L∞ (Ω ). Let Tk (vn − φ) be the test function in (26), which means that ∫ ∫ ∇vn · ∇Tk (vn − φ)dx + (fn − un )Tk (vn − φ)dx = 0. (29) Ω

Ω

By Lebesgue’s dominated convergence theorem, we can pass to the limit n → +∞ in the second term of (29) to find ∫ ∫ lim (fn − un )Tk (vn − φ)dx = (f − u)Tk (v − φ)dx. n→+∞

Ω

Ω

The first term of (29) can be rewritten as ∫ ∫ ∇Tk (vn − φ) · ∇Tk (vn − φ)dx + ∇φ∇Tk (vn − φ)dx. Ω

Ω

458

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

The first term is nonnegative. Then by Fatou’s lemma, ∫ ∫ ∇Tk (v − φ) · ∇Tk (v − φ)dx ≤ lim inf ∇Tk (vn − φ)∇Tk (vn − φ)dx. n→+∞

Ω

Ω

While for the second term ∫ ∇φ∇Tk (vn − φ)dx =

lim

n→+∞

∫

Ω

∇φ · ∇Tk (v − φ)dx. Ω

Substituting these results into (29) we finally get ∫ ∫ ∇v · ∇Tk (v − φ)dx + (f − u)Tk (v − φ)dx ≤ 0, Ω

Ω

which means that v is an entropy solution of problem (20)–(21). 5. At last, we show that the solution is unique. In the above proof, we have constructed an entropy solution v, which is the limit of the solution vn to the approximate problem (24)–(25). Let v˜ be another entropy solution of problem (20)–(21). We show that v˜ = v a.e. in Ω . By fn ∈ L∞ (Ω ) and un ∈ L∞ (Ω ), we have vn ∈ H01 (Ω ) ∩ L∞ (Ω ). Then we can choose φ = vn in the definition of entropy solution, which means that ∫ ∫ ∇˜ v · ∇Tk (˜ v − vn )dx + (f − u)Tk (˜ v − vn )dx ≤ 0. Ω

Ω

On the other hand, Tk (˜ v − vn ) can be chosen as the test function in the weak formulation for vn , which means that ∫ ∫ ∇vn · ∇Tk (˜ v − vn )dx + (fn − un )Tk (˜ v − vn )dx = 0. Ω

Ω

By subtracting the above results, we have ∫ ∫ 2 |∇Tk (˜ v − vn )| dx = ∇(˜ v − vn ) · ∇Tk (˜ v − vn )dx Ω ∫ Ω ≤ − (f − fn − u + un )Tk (˜ v − vn )dx ≤ k(∥f − fn ∥L1 + ∥u − un ∥L1 ). Ω

Repeating step 2 for v˜ − vn , f − fn and u − un , we have ∫ q |∇(˜ v − vn )| dx ≤ C(∥f − fn ∥L1 + ∥u − un ∥L1 )q , Ω d d−1 .

Let n → ∞ and recalling that vn → v in W 1,q (Ω ), we obtain ∇˜ v = ∇v. This together with for 1 ≤ q < ∫ ∫ vdx = v ˜ dx = 0 finish the proof of the uniqueness of solutions. □ Ω Ω 4. Proof of the main result In this section, we prove Theorem 3. The existence of weak solutions for system (1)–(3) is proved by Schauder’s fixed-point theorem. For the proof of the uniqueness, we utilize the concept of entropy solutions. 4.1. The uniqueness of weak solutions Let (u1 , v1 ) and (u2 , v2 ) be two weak solutions of system (1)–(3). We shall prove u1 = u2 and v1 = v2 a.e. in Ω . 1. Clearly, for any ϕ ∈ W 1,p(x) (Ω ), Tk (u − ϕ) can be chosen as the test function in (10), i.e., ∫ ∫ p(x)−2 g(x)|∇u| ∇u · ∇Tk (u − ϕ)dx + λ vTk (u − ϕ)dx = 0. (30) Ω

Ω

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

459

Let h > 0 be a constant. We write (30) with respect to solution (u1 , v1 ) with ϕ = Th (u2 ) and with respect to solution (u2 , v2 ) with ϕ = Th (u1 ). Adding up both results, we get ∫ g(x)|∇u1 |p(x)−2 ∇u1 · ∇Tk (u1 − Th (u2 ))dx Ω ∫ + g(x)|∇u2 |p(x)−2 ∇u2 · ∇Tk (u2 − Th (u1 ))dx (31) Ω∫ ∫ = −λ v1 Tk (u1 − Th (u2 ))dx − λ v2 Tk (u2 − Th (u1 ))dx. Ω

Ω

2. In the following, we conclude ∇u1 = ∇u2 by passing to the limit h, k → +∞ in (31) and disregarding some positive terms. For this, define E1 := {|u1 − Th (u2 )| ≤ k, |u2 | ≤ h}, E2 := E1 ∩ {|u1 | ≤ h},

E3 := E1 ∩ {|u1 | > h}.

Then the first term in (31) ∫ p(x)−2 = g(x)|∇u1 | ∇u1 · ∇(u1 − Th (u2 ))dx {|u1 −Th (u2 )|≤k} ∫ p(x)−2 ≥ g(x)|∇u1 | ∇u1 · ∇(u1 − u2 )dx ∫E1 ∫ p(x)−2 p(x)−2 ≥ g(x)|∇u1 | ∇u1 · ∇(u1 − u2 )dx − g(x)|∇u1 | ∇u1 · ∇u2 dx. E2

E3

Considering the last term in the above inequality, we have ∫ | g(x)|∇u1 |p(x)−2 ∇u1 · ∇u2 dx| E3

≤ C∥|∇u1 |p(x)−2 ∇u1 ∥p′ (·),{h<|u1 |
where ε(h) converges to zero as h → +∞. 3. Repeating step 2 for the second term of (31) yields ∫ p(x)−2 g(x)|∇u2 | ∇u2 · ∇(u2 − Th (u1 ))dx {|u2 −Th (u1 )|≤k} ∫ p(x)−2 ≥ g(x)|∇u2 | ∇u2 · ∇(u2 − u1 )dx + ε(h).

(33)

E2

Then combining (31)–(33) and letting h → +∞, we have ∫ p(x)−2 g(x)|∇u1 | ∇u1 · ∇(u1 − u2 )dx ∫{|u1 −u2 |
(34)

460

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

4. We repeat step 1–3 for Eq. (11) and prove similarly that ∫ ∫ ∇v1 · ∇(v1 − v2 )dx − ∇v2 · ∇(v1 − v2 )dx {|v1 −v2 |
(35)

Ω

By Lebesgue’s dominated convergence theorem, ∫ ∫ lim (u1 − u2 )Tk (v1 − v2 )dx − (v1 − v2 )Tk (u1 − u2 )dx = 0. k→+∞

Ω

Ω p−2

Combining this with (34)–(35), recalling the monotonicity of the operator |s| s and letting k → +∞, we have ∇u1 = ∇u2 and ∇v1 = ∇v2 a.e. in Ω . 5. Now we finish the proof of the uniqueness by the definition of weak solutions. By substituting (u1 , v1 ), (u2 , v2 ) and ∇u1 = ∇u2 into (10), we have ∫ (v1 − v2 )ϕdx = 0, Ω

for any ϕ ∈ W 1,p(x) (Ω ). This implies that v1 = v2 a.e. in Ω . Similarly, substituting (u1 , v), (u2 , v) into (11), we have ∫ (u1 − u2 )Tk (v − φ)dx = 0, Ω

for any φ ∈ H 1 (Ω ) ∩ L∞ (Ω ). This implies that u1 = u2 a.e. in Ω . 4.2. The existence of weak solutions Let

N }, N −1 where the constant C0 > 0 will be determined in the following. Then W is a non-empty, convex and weakly 1,p(x) compact subset of W 1,q (Ω ). For any w ∈ W , by Theorem 4, there exists a u ∈ W0 (Ω ), which is the unique solution of subproblem (13)–(14). Besides, W = {w ∈ W01,q (Ω ) : ∥w∥1,q ≤ C0 , 1 ≤ q <

∥u∥L1 ≤ C∥u∥1,p(·) ≤ C1 ∥λw∥s(p∗ )′ ≤ C1 λs C0s . −

Next, according to Theorem 5, there exists a v ∈ W01,q (Ω ), which is the unique solution of subproblem (20)–(21). Furthermore, ∥v∥1,q ≤ C2 (∥f ∥L1 + ∥u∥L1 ) ≤ C2 (∥f ∥L1 + C1 λs C0s ) =: R. Let λ > 0 be sufficiently small. Then there always exists some C0 such that R ≤ C0 , which means that v ∈ W . We define a fixed-point mapping S : W → W with S(w) = v =: vw . In order to utilize Schauder’s fixed-point theorem, we need only to show that S : w → vw is weakly continuous (W → W ). Let {wj } be a sequence that converges weakly to some w ∈ W . We have to prove that S(wj ) = vj converges weakly to S(w) = vw . From (16), (23) and the classical results of compact inclusion in Sobolev spaces, we can extract from {wj }, respectively from {uj } and from {vj }, subsequences (still denoted by themselves) such that for some u and v, wj ⇀ w, weakly in W 1,q (Ω ), uj ⇀ u, weakly in W 1,p(x) (Ω ), a.e. in Ω ,

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

461

vj ⇀ v, weakly in W 1,q (Ω ), a.e. in Ω , Tk (vj ) ⇀ Tk (v), weakly in H 1 (Ω ),

∀k > 0,

∇vj ⇀ ∇v, a.e. in Ω . Notice that vn − vm satisfies ∆(vn − vm ) + un − um = 0 with homogeneous Neumann boundary condition, we conclude from (23) that {vj } is a Cauchy sequence in W 1,q (Ω ), which leads us to the last convergence result listed above. Since uj and vj are solutions of subproblems (13)–(14) and (20)–(21), we have ∫ ∫ p(x)−2 g(x)|∇uj | ∇uj · ∇ϕdx − λ wj ϕdx = 0, (36) Ω

for any ϕ ∈ W 1,p(x) (Ω ), and ∫

Ω

∫ ∇vj · ∇Tk (vj − φ)dx −

Ω

(f − uj )Tk (vj − φ)dx ≤ 0,

(37)

Ω

for any φ ∈ H 1 (Ω )∩L∞ (Ω ). The above convergence results allow us to pass to the limit j → +∞ in (36)–(37) ∫ ∫ and obtain v = vw = S(w). In fact, according to Remark 1, we have limj→∞ Ω wj ϕdx = Ω wϕdx. By the same method we used in the step 5 of the proof of Theorem 4, we deduce from (36) that ∫ ∫ p(x)−2 g(x)|∇u| ∇u · ∇ϕdx − λ wϕdx = 0. Ω

Ω

By the convergences uj → u, vj → v and ∇vj → ∇v a.e. in Ω and Lebesgue’s dominated convergence theorem, we let j → +∞ in (37) to find ∫ ∫ ∇v · ∇Tk (v − φ)dx − (f − u)Tk (v − φ)dx ≤ 0. Ω

Ω

The above results show that v = vw = S(w). At last, since the solution is unique, the whole sequence vj = S(wj ) converges weakly in W to vw = S(w), which means that S is weakly continuous. Consequently, by Schauder’s fixed-point theorem, there exists a 1,p(x) v ∈ W such that v = S(v). Then (u, v) ∈ (W0 (Ω ), W01,q (Ω )) is a weak solution of system (1)–(3), where u is the solution of subproblem (13)–(14) with w replaced by v. Clearly, u ∈ W 1,p(x) (Ω ) and v ∈ W 1,q (Ω ). This finishes the proof. 5. Numerical implementation This section focuses on the numerical implementation of the proposed elliptic system (1)–(3) for image denoising. The proof of the existence of weak solutions also provides an iterative scheme for solving the system numerically. We shall construct a finite difference scheme for the system by using the iterative scheme and the lagged diffusivity fixed-point iteration [24]. 5.1. The iterative scheme For any given v (0) , let u(1) be the solution of subproblem (13)–(14) with w = v (0) . Then we find a v (1) by solving subproblem (20)–(21) with u = u(1) . Iteratively, we construct a sequence of functions {(u(n) , v (n) )}, which converges to the weak solution (u, v) of system (1)–(3) in the following sense.

462

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

Theorem 6. Let f ∈ L2 (Ω ), 1 < p− ≤ p(x) ≤ p+ = 2, N = 2, the initial guess v (0) ∈ L2 (Ω ) and λ > 0 be sufficiently small. Then ∥∇u − ∇u(n+1) ∥p− ≤ C∥∇u − ∇u(n) ∥p− , (38) where the constant 0 < C < 1 depends only on λ, g− , p− , Ω . Thus {u(n) } converges to u in W 1,p− (Ω ). Proof . 1. As a first step, we show the uniformly boundedness of ∥∇u(n) ∥p(·) for small λ > 0. According to the construction, u(n+1) is the weak solution of the p(x)-Laplacian equation (13)–(14) with w = v (n) , and v (n) is the weak solution of the Poisson equation (20)–(21) with u = u(n) . Recalling the results of Theorem 4 (step 3 of the proof) and Theorem 5, and noticing that (p∗− )′ < 2, 1 ≤ q < 2, we have p −1

p −1

+ − min{∥∇u(n+1) ∥p(·) , ∥∇u(n+1) ∥p(·) } ≤ λC∥v (n) ∥(p∗− )′ ≤ λC1 ∥v (n) ∥2 ,

(39)

and ∥v (n) ∥2 ≤ C∥v (n) ∥1,q ≤C(∥f ∥1 + ∥u(n) ∥1 ) ( ≤C2

∥f ∥1 +

p+

p−

)

p+ p− , ∥∇u(n) ∥p(·) } max{∥∇u(n) ∥p(·)

(40)

.

Here the constants C1 and C2 depend only on g− , p− , Ω . We now let n = 0 in (39) and let λ > 0 be small such that λC1 ∥v (0) ∥2 ≤ 1. Then we have ∥∇u(1) ∥p(·) ≤ 1. We let n = 1 in (40) to find ∥v (n) ∥2 ≤ C2 (∥f ∥1 + 1) =: C3 . Then we return to (39) with n = 1 and λ ≤ C11C3 to obtain ∥∇u(2) ∥p(·) ≤ 1. By iteration, we conclude that ∥∇u(n) ∥p(·) ≤ 1 for all integer n > 0 if λ > 0 is sufficiently small. 2. Since f ∈ L2 (Ω ), by the same method used in Section 3, there exist u ∈ W 1,p(x) (Ω ) and v ∈ H 1 (Ω ) satisfying system (1)–(3) in the weak sense. Then we have ∫ ( ) p(x)−2 p(x)−2 g(x) |∇u| ∇u − |∇u(n+1) | ∇u(n+1) ∇ϕdx Ω ∫ (41) (n) +λ (v − v )ϕdx = 0, Ω

∫ (

)

∇v − ∇v (n+1) ∇φdx −

Ω 1,p(x)

∫

(u − u(n+1) )φdx = 0,

(42)

Ω

for any integer n ≥ 0, ϕ ∈ W (Ω ) and φ ∈ H 1 (Ω ). (n+1) Let ϕ = ∇u − ∇u in (41) and utilizing the monotonicity, we have ∫ ( ) p(x)−2 2 p(x) p(x) p(x) L := g− (p− − 1)|∇u − ∇u(n+1) | |∇u| + |∇u(n+1) | dx

(43)

Ω

≤λ∥v − v

(n)

∥

(p∗ )′ −

(n+1)

∥u − u

∥

p∗ −

≤ λC∥v − v

(n)

∥2 ∥∇u − ∇u

(n+1)

∥ p− .

By the reverse H¨ older inequality, (∫ ) p2 − (n+1) p− L ≥ g− (p− − 1) |∇u − ∇u | dx Ω

(∫ ·

(|∇u|

p(x)

) p−p −2

+

p p(x) p(x)−2 p − |∇u(n+1) | ) p(x) − −2 dx

−

.

Ω p −2

p(x) − We now estimate the last term in the above. Notice that s(x) := p(x)−2 older’s p− ≥ 1, we deduce from H¨ inequality that ∫   p(x) 1 p(x) 1   p(x) p(x) (|∇u| + |∇u(n+1) | ) s(x) dx ≤ C (|∇u| + |∇u(n+1) | ) s(x)  . Ω

s(x)

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

463

Applying (9) to the right-hand side and recalling that p+ = 2, we have ∫ p(x) 1 p(x) (|∇u| + |∇u(n+1) | ) s(x) dx ≤ C max {Mn+1 , 1}, Ω

∫

p(x)

where Mn := Ω |∇u| + |∇u(n) | three inequalities, we obtain

p(x)

dx is uniformly bounded by the result of step 1. Combining the above

(∫ L≥C

(n+1) p−

|∇u − ∇u

|

) p2

−

dx

p− −2 p−

max {Mn+1 , 1}

Ω

(∫ ≥C

p−

|∇u − ∇u(n+1) |

) p2 dx

−

,

Ω

which together with (43) lead us to ∥∇u − ∇u(n+1) ∥p− ≤ λC∥v − v (n) ∥2 .

(44)

Now let φ = v − v (n+1) in (42) and utilize H¨ older’s inequality, we have ∥∇v − ∇v (n+1) ∥2 ≤ ∥u − u(n+1) ∥2 ≤ C∥u − u(n+1) ∥p∗− ≤ C∥∇u − ∇u(n+1) ∥p− . By applying Poincar´e’s inequality to the left-hand side of the above inequality, we have ∥v − v (n+1) ∥2 ≤ C∥∇u − ∇u(n+1) ∥p− . □

This and (44) finish the proof of (38). 5.2. The finite difference scheme

We now present the finite difference scheme for system (1)–(3). According to the iterative scheme, we need to discretize the p(x)-Laplacian equation (1) and the Poisson equation (2) with homogeneous Neumann boundary conditions, respectively. For the p(x)-Laplacian equation, a computational regularization is needed p(x)−2 to overcome the singularity of the nonlinear gradient term |∇u| . We adopt the common technique that p(x)−2 replacing the gradient term by |∇u + ε| , where ε is a small constant. Let M ×N be the size of the image and ∆x = ∆y = 1 be the space step. Then the finite difference scheme for the equation and the boundary condition is given by Di− 1 ,j (ui−1,j − ui,j ) + Di,j− 1 (ui,j−1 − ui,j ) + Di+ 1 ,j (ui+1,j − ui,j ) 2

2

2

+ Di,j+ 1 (ui,j+1 − ui,j ) = λvi,j , 2

u0,j = u1,j ,

uM +1,j = uM,j ,

ui,0 = ui,1 ,

(45)

ui,N +1 = ui,N ,

where 1 ≤ i ≤ M , 1 ≤ j ≤ N and Di− 1 ,j = gi− 1 ,j |∇ui− 1 ,j + ε| 2

2

p

i− 1 2 ,j

2

−2

,

gi− 1 ,j = (gi−1,j − gi,j )/2, 2

∇ui− 1 ,j = ((ux )i− 1 ,j , (uy )i− 1 ,j ) = (ui,j − ui−1,j , (uy )i− 1 ,j ), 2

2

2

2

(uy )i− 1 ,j = (ui−1,j+1 + ui,j+1 − ui−1,j−1 − ui,j−1 )/4. 2

Other functions evaluated on midpoints are defined in a similar way. While for the Poisson equation, the difference finite scheme is given as vi−1,j + vi,j−1 + vi+1,j + vi,j+1 − 4vi,j = fi,j − ui,j , v0,j = v1,j ,

vM +1,j = vM,j ,

vi,0 = vi,1 ,

vi,N +1 = vi,N ,

(46)

464

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

where 1 ≤ i ≤ M , 1 ≤ j ≤ N . Both systems of algebraic equations (45) and (46) will be solved by iterative solvers. In order to solve the nonlinear system (45) efficiently, we use the lagged diffusivity fixed-point iteration to linearize the system at each iteration and use the Gauss–Seidel iteration to solve the linearized system. System (46) can be solved by the Gauss–Seidel iteration directly [25]. We summarize the finite difference scheme for system (1)–(3) as follows. (n+1)

(n+1) ui,j

=

(n+1)

2

2

(n)

2

2

Di− 1 ,j + Di,j− 1 + Di+ 1 ,j + Di,j+ 1 2

(n+1)

(n+1) vi,j

(n)

(n)

Di− 1 ,j ui−1,j + Di,j− 1 ui,j−1 + Di+ 1 ,j ui+1,j + Di,j+ 1 ui,j+1 − λvi,j (n+1)

2

(n)

2

(n)

,

(47)

2

(n+1)

vi−1,j + vi,j−1 + vi+1,j + vi,j+1 − fi,j + ui,j = 4

.

(48) (0)

Here we use an alternating direction scheme to speed up it. That is to say, for the initial guesses ui,j = fi,j (0) (1) (1) and vi,j = 0, we find ui,j by solving (47) with n = 0. Then we find vi,j by solving (48) with n = 1. We return to (47) with n = 1 and proceed iteratively until the best result is obtained. 5.3. Numerical experiments The variable exponent p(x) in the p(x)-Laplacian equation is usually taken as p(x) = 1 +

1

2,

1 + k|(∇Gσ ∗ f )(x)|

(49)

for image processing. Consequently, the p(x)-Laplacian equation acts like the TV equation at the image’s boundaries and acts as the heat equation at homogeneous regions. It avoids the staircasing effect appearing in the TV regularization, but still suffers from two problems. First, edges are not preserved as good as in the TV regularization. Second, noise in an image f (x) will be wrongly regarded as boundaries. To fix these problems, we propose to replace the smooth variable exponent by a discontinuous one. We shall construct a piecewise constant variable exponent from p(x) defined in (49) by utilizing the local K-means algorithm [26]. More precisely, the image of p(x) in (49) is split into a series of equal size subimages. And each subimage is converted into a piecewise constant image by using the K-means algorithm (K = 2). The advantages of the piecewise constant variable exponent can be illustrated by a concrete experiment clearly. Fig. 1 shows the denoising results of system (1)–(3) with smooth variable exponent and piecewise constant variable exponent for an artificial image. It can be observed that the result with piecewise constant variable exponent has better visual quality and higher PSNR value than the result with smooth variable exponent. Here PSNR (Peak Signal-to-Noise Ratio) for two images fi,j and ui,j is defined as PSNR(f, u) = 10 log10

1 MN

2552 db. 2 i,j (ui,j − fi,j )

∑

Fig. 2 shows the images of both variable exponents. The background of Fig. 2(b) is much cleaner than the background of Fig. 2(a). Besides, the black rings in Fig. 2(b) are wider than the corresponding rings in Fig. 2(a). These explain the results in Figs. 1(c) and 1(d). At last , we compare the proposed elliptic system (1)–(3) with three related classical methods. The first one is the non-standard growth variation method, which corresponds to the p(x)-Laplacian equation [5]. The second one is the TV method [3], which is famous for its edge preserving ability. The last one is the PM method [27], which is a forward–backward diffusion equation. Figs. 3 and 4 show the denoising results of the four methods for two real images. The results of the TV method suffer from the staircasing effect. The PM method enhances images in regions where |∇u| is large. Consequently, some noise is remaining in the results. The p(x)-Laplacian equation avoids the staircasing effect. However, it wrongly detects some noise as edges. The results of the proposed elliptic system have the best visual quality and the highest PSNR values.

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

465

Fig. 1. The denoising results of system (1)–(3) with smooth variable exponent and piecewise constant variable exponent for a 300 × 300 artificial image.

Fig. 2. The images of the two variable exponents, σ = 0.5, k = 0.02.

466

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

Fig. 3. The denoising results of the four methods. The regions in read squares are enlarged for better visual comparison. (c) λ = 0.1, σ = 0.5 and k = 0.02 in p(x), 60 iterations. (d) λ = 0.1, σ = 0.5 and k = 0.02 in p(x), 20 iterations. (e) 130 iterations. (f) K = 11, 60 iterations.

Acknowledgments The authors would like to thank the reviewer for the valuable suggestions and comments that helped to improve the initial version of the manuscript significantly.

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

467

Fig. 4. The denoising results of the four methods. The regions in read squares are enlarged for better visual comparison. (c) λ = 0.1, σ = 0.5 and k = 0.02 in p(x), 120 iterations. (d) λ = 0.1, σ = 0.5 and k = 0.02 in p(x), 30 iterations. (e) 230 iterations. (f) K = 11, 60 iterations.

This work was partially supported by the National Science Foundation of China (Grant Nos. U1637208, 11871133), the National Key Research and Development Program of China (2017YFB1401801) and the Natural Science Foundation of Heilongjiang Province of China (Grant Nos. LC2018001, G2018006).

468

D. Zhang, K. Shi, Z. Guo et al. / Nonlinear Analysis: Real World Applications 50 (2019) 448–468

References [1] S. Osher, A. Sol´ e, L. Vese, Image decomposition and restoration using total variation minimization and the H −1 norm, Multiscale Model. Simul. 1 (3) (2003) 349–370. [2] R. Dautray, J.-L. Lions, Mathematical Analysis and Numerical Methods for Science and Technology, in: Physical Origins and Classical Methods, vol. 1, Springer Science & Business Media, 2012. [3] L.I. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D 60 (1–4) (1992) 259–268. [4] P. Blomgren, T.F. Chan, P. Mulet, C.K. Wong, Total variation image restoration: numerical methods and extensions, in: International Conference on Image Processing, 1997. Proceedings, Vol. 3, 1997, pp. 384–387. [5] Y. Chen, S. Levine, M. Rao, Variable exponent, linear growth functionals in image restoration, SIAM J. Appl. Math. 66 (4) (2006) 1383–1406. [6] Z. Guo, Q. Liu, J. Sun, B. Wu, Reaction–diffusion systems with p(x)-growth for image denoising, Nonlinear Anal. RWA 12 (5) (2011) 2904–2918. [7] T.F. Chan, S. Esedoglu, Aspects of total variation regularized L1 function approximation, SIAM J. Appl. Math. 65 (5) (2005) 1817–1837. [8] P. Athavale, E. Tadmor, Integro-differential equations based on (BV, L1 ) image decomposition, SIAM J. Imaging Sci. 4 (1) (2011) 300–312. [9] Z. Guo, J. Yin, Q. Liu, On a reaction–diffusion system applied to image decomposition and restoration, Math. Comput. Modelling 53 (5) (2011) 1336–1350. [10] L. Afraites, A. Atlas, F. Karami, D. Meskine, Some class of parabolic systems applied to image processing, Discrete Contin. Dyn. Syst. 21 (6) (2016) 1671–1687. [11] Q. Liu, Z. Guo, C. Wang, Renormalized solutions to a reaction–diffusion system applied to image denoising, Discrete Contin. Dyn. Syst. 21 (6) (2016). [12] P. B´ enilan, L. Boccardo, T. Gallou¨ et, R. Gariepy, M. Pierre, J.L. V´ azquez, An L1 -Theory of Existence and Uniqueness of Solutions of Nonlinear Elliptic Equations, Citeseer, 1995. [13] M. Sanch´ on, J.M. Urbano, Entropy solutions for the p(x)-Laplace equation, Trans. Amer. Math. Soc. 361 (12) (2009) 6387–6405. [14] R.J. DiPerna, P.-L. Lions, On the Cauchy problem for Boltzmann equations: global existence and weak stability, Ann. of Math. (1989) 321–366. [15] M. Bendahmane, P. Wittbold, Renormalized solutions for nonlinear elliptic equations with variable exponents and L1 data, Nonlinear Anal. TMA 70 (2) (2009) 567–583. [16] C. Zhang, S. Zhou, Renormalized and entropy solutions for nonlinear parabolic equations with variable exponents and L1 data, J. Differential Equations 248 (6) (2010) 1376–1400. [17] V.V. Zhikov, On Lavrentiev’s phenomenon, Russ. J. Math. Phys. 3 (2) (1995). [18] O. Kov´ aˇ cik, J. R´ akosn´ık, On spaces Lp(x) and W k,p(x) , Czechoslovak Math. J. 41 (4) (1991) 592–618. [19] P. Harjulehto, P. H¨ ast¨ o, An overview of variable exponent lebesgue and Sobolev spaces, in: D. Herron (Ed.), Future Trends in Geometric Function Theory, Jyv¨ askyl¨ a, 2003, pp. 85–93, RNC Work-shop. [20] S. Samko, On a progress in the theory of lebesgue spaces with variable exponent: maximal and singular operators, Integral Transforms Spec. Funct. 16 (5–6) (2005) 461–482. [21] M. Bul´ıˇ cek, A. Glitzky, M. Liero, Systems describing electrothermal effects with p(x)-Laplacian-like structure for discontinuous variable exponents, SIAM J. Math. Anal. 48 (5) (2016) 3496–3514. [22] L. Orsina, Elliptic equations with measure data, http://www.ugr.es/∼edp/Talks/Luigi-Abril 2009/Curso.pdf. [23] X.-L. Fan, Q.-H. Zhang, Existence of solutions for p(x)-Laplacian dirichlet problem, Nonlinear Anal. TMA 52 (8) (2003) 1843–1852. [24] T. Chan, J. Shen, Image Processing and Analysis: Variational, PDE, Wavelet, and Stochastic Methods, Society for Industrial and Applied Mathematics, 2005. [25] J.W. Thomas, Numerical Partial Differential Equations: Finite Difference Methods, Springer, 1995. [26] K. Alsabti, S. Ranka, V. Singh, An efficient k-means clustering algorithm, 1997. [27] P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Trans. Pattern Anal. Mach. Intell. 12 (7) (1990) 629–639.