Finite difference schemes for the parabolic p-Laplace equation

We propose a new finite difference scheme for the degenerate parabolic equation ∂tu-div(|∇u|p-2∇u)=f,p≥2.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \partial _t u - \text{ div }(|\nabla u|^{p-2}\nabla u) =f, \quad p\ge 2. \end{aligned}$$\end{document}Under the assumption that the data is Hölder continuous, we establish the convergence of the explicit-in-time scheme for the Cauchy problem provided a suitable stability type CFL-condition. An important advantage of our approach, is that the CFL-condition makes use of the regularity provided by the scheme to reduce the computational cost. In particular, for Lipschitz data, the CFL-condition is of the same order as for the heat equation and independent of p.

The main result is the pointwise convergence of our scheme given Hölder continuous data (f and u 0 ) and a stability type CFL-condition.See Theorem 2.2 for the precise statement and (CFL) for the CFLcondition.One of the advantages of our approach is that the CFL-condition makes use of the regularity provided by the scheme.As a consequence, for Lipschitz continuous data, the CFL-condition is of the same order as the one for the heat equation.In general, the order of the CFL-condition depends on p and on the regularity of the data.
1.1.Related literature.Equation (1.1) has attracted much attention in the last decades.We refer to [11] and [12] for the theory for weak solutions of this equation and to [23] for the relation between viscosity solutions and weak solutions.To the best of our knowledge, the best regularity results known are C 1,α −regularity in space for some α > 0 (see [11, Chapter IX]) and C 0,1/2 −regularity in time (see [3,Theorem 2.3]).
The literature regarding finite difference schemes for parabolic problems involving the p-Laplacian is quite scarce.One reason for that is naturally that, since the p-Laplacian is in divergence form, it is very well suited for methods based on finite elements, see for instance [2,21,14,1,19] for related results.
In the stationary setting, there has been some development of finite difference methods the past 20 years.Section 1.1 in [28] provides an accurate overview of such results, we will only mention a few.In [5,10,18,28], finite difference schemes for the p-Laplace equation based on the mean value formula for the normalized p-Laplacian (cf.[25]) are considered.Since the corresponding parabolic equation for the normalized p-Laplacian is completely different in nature (see [22,15]), these methods do not seem very well suited to be used for the parabolic equation considered in this paper.In [9], the authors of the present paper studied a monotone finite difference discretization of the p-Laplacian based on the mean value property presented in [8,4].We also seize the opportunity to mention [27], where difference schemes for degenerate elliptic and parabolic equations (but not for equation (1.1)) are discussed.
It is noteworthy that, in dimension d = 1, the spatial derivative of a solution of (1.1) is a solution of the Porous Medium Equation (PME).See [29,30] for a general presentation of the PME, and [20] for a proof of this fact.Finite difference schemes for the PME are well known, see [13,16,26,6,7].

Assumptions and main results
In this section, we introduce a general form of finite difference discretizations of ∆ p and the associated numerical scheme for (1.1).This is followed by our assumptions, the notion of solutions for (1.1) and the formulation of our main result.
2.1.Discretization and scheme.In order to treat (1.1), we consider a general discretization of ∆ p of the form (2.1) We also need to introduce a time discretization.We will employ an explicit and uniform-in-time discretization.Let N ∈ N and consider a discretization parameter τ > 0 given by τ = T /N .Consider also the sequence of times {t j } N j=0 defined by t 0 = 0 and t j = t j−1 + τ = jτ .The time grid, T τ , is given by Then, our general form of an explicit finite difference scheme of (1.1) is given by (2.2) where f α := f (x α ), (u 0 ) α = u 0 (x α ) and D h p is given by (2.1).

2.2.
Assumptions.In order to ensure convergence of the scheme (2.2), we impose the following hypotheses on the data and the discretization parameters.This entails a regularity assumption on the data, some assumptions on the discretization and a nonlinear CFL-condition on the parameters, as is customary for explicit schemes.
Hypothesis on the data.We assume that (A u0,f ) u 0 , f : R d → R are bounded and globally Hölder continuous functions for some a ∈ (0, 1].
Hypothesis on the spatial discretization.For the discretization, we assume the following type of monotonicity and boundedness: (A ω ) ω β = ω −β ≥ 0, w β = 0 for y β ∈ B r for some r > 0, and Here M = M (p, d) > 0. In addition, we assume the following consistency for the discretization: Examples of discretizations satisfying these properties can be found in Section 5.
Hypothesis on the discretization parameters.We assume the following stability condition on the numerical parameters: with and K a constant given in (3.4), depending on p, the modulus of continuity in time of the discretized solution and some universal constants coming from a mollifier.
Remark 2.1.For Lipschitz data u 0 and f , the condition (CFL) reads τ ≤ Cr 2 for a certain constant C = C(u 0 , f, d, p, T ) > 0. We note that, regardless of the constant C, the relation between τ and r is always quadratic (as in the linear case p = 2) and independent of p.It is important to mention that this is computationally very relevant, especially if we want to deal with problems related to large p.

Main result.
We now state our main result regarding the convergence of the scheme.Several other properties of the scheme are also obtained, but we will state them later.
2.4.Viscosity solutions.Throughout the paper, we will use the notion of viscosity solutions.For completeness, we define the concept of viscosity solutions of (1.1), adopting the definition in [23].
Remark 2.3.We remark that it is not necessary to require strict inequality in the definition above.It is We also state a necessary uniqueness result that will ensure convergence of the scheme.Without such a result, we would only be able to establish convergence up to a subsequence.The theorem below is a consequence of the fact that viscosity solutions are weak solutions (see Corollary 4.7 in [23]) and that bounded weak solutions are unique (see Theorem 6.1 in [11]).

Properties of the numerical scheme
In this section we will study properties of the numerical scheme (2.2).More precisely, we establish existence and uniqueness for the numerical solution, stability in maximum norm, as well as conservation of the modulus of continuity of the data.
3.1.Existence and uniqueness.We have the following existence and uniqueness result for the numerical scheme.Proposition 3.1.Assume (A u0,f ), (A ω ), p ≥ 2 and r, h, τ > 0. Then there exists a unique solution Proof.First we note that, for a function ψ ∈ ∞ (G h ), we have that Then, for each α ∈ Z, U j α is defined recursively using the values of U j−1 β for β ∈ Z, and we have that sup The conclusion follows since sup 3.2.Stability and preservation of the modulus of continuity in space.First we will prove that the scheme preserves the regularity of the data.
For every j = 0, . . ., N , we have Remark 3.3.In particular, if both u 0 and f are Lipschitz functions with constants L u0 and L f respectively, the above result reads, Proof of Proposition 3.2.By assumption (A u0,f ), for any given x α , x γ ∈ G h , we have that Assume by induction that Using the scheme at x α and x γ we get Now, since p ≥ 2, we have, by Taylor expansion, that for some η β ∈ R between (U j α+β − U j α ) and (U j γ+β − U j γ ).Thus, Now observe that, by the induction assumption, we have By (A ω ), we have w β = 0 for y β ∈ B r for some r > 0, and we deduce that Thus, by (CFL), we get τ (p − 1) Using the above estimate and the induction hypothesis in (3.1), we get that which concludes the proof.
We are now ready to state and prove the stability result: solutions with bounded data remain bounded (uniformly in the discretization parameters) for all times.Proposition 3.4.Under the assumptions of Proposition 3.2, we have that , for all j = 0, . . ., N.
Proof.By assumption (A u0,f ), we have that sup Assume by induction that sup Direct computations lead to By Proposition 3.2 we have that , which together with assumptions (A ω ) and (CFL) imply that Direct computations plus the induction hypothesis allow us to conclude that which concludes the proof.
3.3.Time equicontinuity for a discrete in time scheme.Now we extend the scheme from G h to R d by considering U : R d × T τ defined by x ∈ R d .
Remark 3.5.Clearly, if we restrict the solution of (3.2) to G h , we recover the solution of (2.2).
Proposition 3.6 (Continuous dependence on the data).Assume (A u0,f ), (A ω ), p ≥ 2, r, h, τ > 0 and (CFL).Let U, U be the solutions of (2.2) corresponding to u 0 , u 0 and f, f .For every j = 0, . . ., N , we have By assumption (A u0,f ), we have that Assume by induction that . Similar computations as the ones in the proof of Proposition 3.2 yield where η β ∈ R is some number between (U j (x + y β ) − U j (x)) and ( U j (x + y β ) − U j (x)).From here, the proof follows as in the proof of Proposition 3.2.Proposition 3.7 (Equicontinuity in time).Assume (A u0,f ), (A ω ), p ≥ 2, r, h, τ > 0 and (CFL).Let U be the solution of (3.2).Then where M comes from assumption (A ω ), and K 1 and K 2 are constants given in Section A.1 (depending on a certain choice of mollifiers).Remark 3.8.Actually, a close inspection of the proof reveals that for Proof.Consider a mollification of the initial data u 0,δ = u 0 * ρ δ where ρ δ (x) is a standard mollifier (as defined in Appendix A).Let (U δ ) j be the corresponding solution of (3.2) with u 0,δ as initial data.Then, Define U j δ := U j+1 δ for all j = 0, . . ., N .Clearly, U j δ is the unique solution of (3.2) with initial data U 0 δ = U 1 δ and right hand side f .By Proposition 3.6 A repeated use of the triangle inequality yields The symmetry of the weights ω β together with Lemma B.1 implies Now note that, by the a-Hölder regularity of u 0 given by assumption (A u0,f ), Lemma A.1 and Lemma A.2 imply where K 1 and K 2 depend only on the mollifier ρ.Now note that, by (A ω ), we have (3.8) Combining (3.5) and (3.8), we obtain Using the triangle inequality, the above estimate and applying Proposition 3.6 several times we obtain By choosing δ = ( K 2) in the above estimate, we get the desired result 3.4.Equiboundedness and equicontinuity estimates for a scheme in R d × [0, T ].We now need to extend the numerical scheme in time in a continuous way.This is done by continuous interpolation, i.e., (3.9) where U j is the solution of (3.2).
Remark 3.9.It is standard to check that, for all t ∈ [t j , t j+1 ], we have that the original scheme is preserved also outside the grid points, i.e., (3.10) We have the following result.Proposition 3.10 (Stability and equicontinuity).Assume (A u0,f ), (A ω ), p ≥ 2, r, h, τ > 0 and (CFL).Let U be the solution of (3.9).Then Proof.Equiboundedness follows easily from a continuous in space version of Proposition 3.4, since Equicontinuity in space follows from the translation invariance of the scheme and Proposition 3.6: To prove equicontinuity in time, we first consider t, t ∈ [t j , t j+1 ] for some j = 0, . . ., N − 1.In this case we have Then, from Proposition 3.7, we get Now consider t ∈ [t j , t j+1 ) and t ∈ [t j+k , t j+k+1 ) for k ≥ 1.By the triangle inequality, the previous step and Proposition 3.7 Since t ≤ t j+1 ≤ t and t ≤ t j+k ≤ t, the above estimate yields Finally, we conclude space-time equicontinuity combining the above estimates to get By Arzelà-Ascoli, we obtain as a corollary that, up to a subsequence, the numerical solution converges locally uniformly to a limit.Corollary 3.11.Assume the hypotheses of Proposition 3.10.Let {U h } h>0 be a sequence of solutions of (3.9).Then, there exist a subsequence

Convergence of the numerical scheme
From Corollary 3.11, we have that the sequence of numerical solutions has a subsequence converging locally uniformly to some function v.We will now show that v is a viscosity solution of (1.1).Proof.For notational simplicity, we avoid the subindex j and consider First of all, by the local uniform convergence, v(x, 0) = lim h→0 U h (x, 0) = u 0 (x), locally uniformly.We will now show that v is a viscosity supersolution.The proof that v is a viscosity subsolution is similar.Now let ϕ be a suitable test function for v at (x * , t * ) ∈ R d × (0, T ).We may assume that ϕ satisfies The local uniform convergence ensures (see Section 10.1.1 in [17]) that there exists a sequence ] (note that the index j might depend on h, but this fact plays no role in the proof).By Remark 3.9, Clearly, U h (x h , t h ) = ϕ(x h , t h ) and U h ≥ ϕ, which implies that (4.1) ϕ(x h , t h ) = U h (x h , t j ) + (t h − t j ) Now consider the function g : R → R given by We will check now that g (ξ) ≥ 0 for any ξ ∈ [ϕ(x h , t j ), U (x h , t j )].Indeed, where we have used that U (x h , t h ) = ϕ(x h , t h ), Proposition 3.10 and the fact that |t h − t j | ≤ τ .By (CFL), and taking τ small enough, we have Thus, where we have used (A ω ) and where the last inequality is due to the (CFL) condition.We can use this fact in (4.1) to get Passing to the limit as h, τ → 0, we get the desired result by the regularity of ϕ and the fact that t h , t j → t * and x h → x * as h → 0.
We are now ready to prove convergence of the scheme.
Proof of Theorem 2.2.By Corollary 3.11 and Theorem 4.1, we know that, up to a subsequence, the sequence U h converges to a viscosity solution of (1.1).Moreover, since viscosity solutions are unique (cf.Theorem 2.4), the whole sequence converges to the same limit.

Discretizations
In this section, we present two examples of discretizations and verify that the assumptions (A c ) and (A ω ) are satisfied.Moreover, we also give the precise form of corresponding CFL-condition.5.1.Discretization in dimension d = 1.We consider the following finite difference discretization of ∆ p in dimension d = 1 A proof of consistency (A c ) can be found in Theorem 2.1 in [8].Assumption (A ω ) is trivially true for r = h since ffl ∂B1 |y 1 | p dσ(y).When p ∈ N, a more explicit value of this constant is given in [9].In general, the explicit value is given by .
A proof of consistency (A c ) can be found in Theorem 1.1 in [9].Assumption (A ω ) trivially holds for h = o(r α ) for some α > 0 according to (5.1) since To check (A ω ) we rely on the following estimate given in the proof of Theorem 1.1 in [9]: In particular, taking for example h ≤ r/ √ d, we have

Numerical experiments
We will perform the numerical tests comparing the numerical solution with the explicit Barenblatt solution of (1.1).For p > 2 this is given by Let us now comment on the CFL-condition (CFL).Clearly, u 0 is a Lipschitz function, and we can estimate its Lipschitz constant as follows Thus, for all p > 2, the CFL condition (CFL) can be take as τ ∼ h 2 (since f = 0 in this case).For completeness, we find the value of K in dimension d = 1.Note that In Figure 1, we show the numerical errors obtained.As it can be seen there, the errors seem to behave like O(h p/(p−1) ).

Theorem 4 . 1 .
Let the assumptions of Corollary 3.11 hold.Then v is a viscosity solution of (1.1).

. 2 .
Discretization in dimension d > 1.The following discretization was introduced in [9]: D h p φ(x) = h d D d,p ω d r p+d y β ∈Br J p (φ(x + y β ) − φ(x)), where ω d denotes the measure of the unit ball in R d , the relation between r and h is given by (if p ∈ (3, ∞), and D d,p = d 2(d+p)