{'fraction_non_alphanumeric': 0.05956339841235786, 'fraction_numerical': 0.019014160051491095, 'mean_word_length': 3.8433562800363683, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Preconditioned Steepest Descent Methods for some Nonlinear Elliptic Equations Involving p-Laplacian Terms

We describe and analyze preconditioned steepest descent (PSD) solvers for fourth and sixth-order nonlinear elliptic equations that include p-Laplacian terms on periodic domains in 2 and 3 dimensions. The highest and lowest order terms of the equations are constant-coefficient, positive linear operators, which suggests a natural preconditioning strategy. Such nonlinear elliptic equations often arise from time discretization of parabolic equations that model various biological and physical phenomena, in particular, liquid crystals, thin film epitaxial growth and phase transformations. The analyses of the schemes involve the characterization of the strictly convex energies associated with the equations. We first give a general framework for PSD in generic Hilbert spaces. Based on certain reasonable assumptions of the linear pre-conditioner, a geometric convergence rate is shown for the nonlinear PSD iteration. We then apply the general the theory to the fourth and sixth-order problems of interest, making use of Sobolev embedding and regularity results to confirm the appropriateness of our pre-conditioners for the regularized p-Lapacian problems. Our results include a sharper theoretical convergence result for p-Laplacian systems compared to what may be found in existing works. We demonstrate rigorously how to apply the theory in the finite dimensional setting using finite difference discretization methods. Numerical simulations for some important physical application problems -- including thin film epitaxy with slope selection and the square phase field crystal model -- are carried out to verify the efficiency of the scheme.

For these time-dependent problems, convex splitting schemes have been proposed and analyzed in [33,36] to obtain unconditional unique solvability and unconditional energy stability. The convex splittings scheme for the thin film model is [36] u m − s∇ · (|∇u m | 2 ∇u m ) + sε 2 ∆ 2 u m = u m−1 − s∆u m−1 , where s > 0 is the time step size, and the superscripts indicate the time discretizations. The convex splitting scheme for the SPFC model -which can be inferred from the general principles in [36,38] is precisely u m − ∆w m = u m−1 , sγ 0 u m − s∇ · |∇u m | 2 ∇u m + sε 2 ∆ 2 u m − w m = − sγ 1 ∆u m−1 , assuming γ 0 , γ 1 ≥ 0. These schemes are nonlinear and require one to deal with the p-Laplacian term at the implicit time level. We remark that there are also second-order-in-time convex splitting schemes for such nonlinear parabolic equations, as described in [33], which have similar nonlinear structures. In any case, solving nonlinear elliptic equations with the p-Laplacian term is challenging, because of its highly nonlinear nature. In [33,36], the authors used a nonlinear conjugate gradient algorithm to solve the nonlinear system at each implicit time step. Such naive gradient methods are guaranteed to converge due to the global convexity of the equations, but are not necessarily efficient.
Several works develop and analyze numerical schemes for nonlinear elliptic equations involving the p-Laplacian operator. The works [3,4,25,29,35,40,41] are based on finite element approximations in space. Recently, the vanishing moment method for the p-Laplacian was proposed in [18]. In that method, the highest order term is purely artificial, whereas, for the models above, the surface diffusion term is small, but non-vanishing. A hybridizable discontinuous Galerkin method for the p-Laplacian was proposed in [11]. Of these works, [25,40,41] are primarily focused on efficient solvers for the elliptic equations with p-Laplacian terms, rather than, say, error estimates.
The main goal of this paper is to design a general framework of preconditioned steepest descent (PSD) methods for certain nonlinear elliptic equations with p-Laplacian terms. The main idea is to use a linearized version of the nonlinear operator as a pre-conditioner, or in other words, as a metric for choosing the search direction. We propose and analyze the preconditioned steepest descent methods for both the fourth-and sixth-order p-Laplacian problems mentioned above. Herein we present numerical simulations for the 6-Laplacian thin film epitaxy and the H −1 gradient flow SPFC model by using the proposed method. While we restrict our focus to the p-Laplacian problems herein, the search direction framework is general and can be applied to other nonlinear equations, such as the Cahn-Hilliard (CH) equation [6,8,31,34], functionalized Cahn-Hilliard (FCH) Equation [9,14,17], for example.
The convergence analyses of the nonlinear iteration algorithms we propose for the p-Laplacian equations are quite challenging, due to the highly nonlinear nature of the problems. However, we are able to recast the equations as equivalent minimization problems involving strictly convex functionals in generic Hilbert spaces. Once this is done, we are able to characterize the properties of general pre-conditioners that will result in geometric convergence rates. This general approach is applicable to both the 4th and 6th order equations at the space-continuous level, as well as the approximation of these problems in finite dimensions using finite differences as we show. Though we do not explore it here, we remark that the theory is extendible to the pseudo-spectral, spectral-Galerkin, and mixed finite element settings as well, using the appropriate discrete Gagliardo-Nirenburg inequalities. To our knowledge, the only related theoretical results available in the existing literature are to be found in [25], in which finite element PSD solvers were designed and analyzed. Specifically, it was proved in [25] that their method converges with the rate O(k −β ), where k is the iteration index and β = p p−2 > 0. In this article, we provide a theoretical analysis with a geometric convergence rate O(α k ), with 0 < α < 1, for the finite difference PSD solver applied to the regularized p-Laplacian problems.
For such nonlinear analyses, the essential difficulty has always been associated with the subtle fact that the numerical solution has to be bounded uniformly in certain functional norms, so that a bound for the iteration error could be established. For the p-Laplacian problems, typically a uniform W 1,p bound of the numerical solution is available at each iteration stage, and such a bound may be used to derive an O(k −β ) convergence rate for the PSD iteration. However, for the regularized p-Laplacian problems, one observes that a linear operator with higher-order diffusion may be utilized so that a uniform H 2 bound of the numerical solution may be obtained. Specifically, the existence of the surface diffusion term ε 2 ∆ 2 u enables us to derive a geometric convergence rate O(α k ) for the PSD iteration, which gives a sharper theoretical result than the existing one in [25]. Our strategy comes at a cost that we point out at the offset: a linear, positive, constant-coefficient operator of order 4 or 6 must be inverted to obtain the search direction. But, since we are interested in applications involving coarsening processes over periodic domains, the FFT can be utilized to make this process efficient.
The remainder of the paper is organized as follows. In section 2, we present a general precondi-tioned steepest descent (PSD) method for nonlinear equations in generic Hilbert spaces, and provide the convergence rate estimates for the PSD method. The application of the general theory to the fourth-order regularized p-Laplacian problem is formulated in section 3. The PSD scheme for the sixth-order regularized p-Laplacian problem is outlined in section 4. Subsequently, in section 5, we introduce a two-dimensional finite difference discretization and provide the fully discrete convergence analysis. Applications to thin film epitaxy with slope selection and the SPFC model and the numerical results are presented in section 6. The concluding remarks are offered in section 7. In the Appendix, we give the proof of a few discrete Sobolev inequalities.

The Classical Setting: Linear SPD Systems in Finite Dimensions
Before we get to the general case, let us quickly review the convergence theory for preconditioned steepest decent methods for solving the linear system Au = f, where A ∈ R m×m sym is positive definite. This is closely related to the preconditioned conjugate gradient (PCG) method, though may be less familiar to the reader. Solving Au = f is, of course, equivalent to minimizing the quadratic energy Here A is the stiffness matrix and L is the pre-conditioner. The idea is that L ≈ A, but the former is "easier to invert." The preconditioned steepest decent algorithm for approximating the solution to Au = f is given in Algorithm 1 [2,26].
if r k+1 < tol or k = k max − 1 then u := u k+1 ; exit for loop; else d k+1 := L −1 r k+1 ; end end Result: u Algorithm 1: Preconditioned Steepest Descent Here d k ∈ R m is called the search direction and r k ∈ R m is called the residual. We observe that We have the classical convergence result: where κ := λm λ 1 , and λ m is the largest eigenvalue of L −1 A, and λ 1 is the smallest [2,26,32].

Non-Quadratic Energy Functionals in Generic Hilbert Spaces
Here we review the general theory for preconditioned steepest descent in a generic Hilbert space [1,10,15,26]. Suppose that H is a (real) Hilbert space with the inner product ( · , · ) H and induced norm · H . We consider an energy functional E[ · ] : H → R with the following properties: and This implies the strict convexity of E.
(E3) E is coercive with respect to the norm on H, i.e., there exist constants If E satisfies (E1) -(E3), it follows [10] that there is is a unique element u ∈ H with the property that and this minimizer further satisfies We wish to construct, via preconditioned steepest descent (PSD), a sequence that converges to the unique minimizer. By H we denote the continuous dual of H. When it is convenient, we use the symbol · , · H : H × H → R to denote the dual pairing between H and H. Consider a linear operator L : H → H . This operator L, which we call the pre-conditioner induces a bilinear form on H: We assume that L satisfies the following properties: (L2) ( · , · ) L is continuous with respect to the standard topology of H, i.e., there is some (L3) ( · , · ) L is coercive with respect to H, i.e., there is some C 4 > 0 such that It follows that ( · , · ) L : H × H → R is an inner product on H, equivalent to the primary inner product ( · , · ) H . The induced norm, ν L := (ν, ν) L , is equivalent to the primary norm. By the Riesz Representation Theorem, if f ∈ H , then there exists a unique u f ∈ H such that where the second norm is the L-induced operator norm. Suppose that u k ∈ H is given. We define the following search direction problem: We call d k the k th search direction. In operator form, we write L[d k ] = −δE[u k ] in H . The functional −δE u k is called the residual of u k . By the Riesz Representation Theorem, we discover that We then define the next iterate u k+1 as where α k ∈ R is the unique solution to Therefore, we have the fundamental orthogonality relation It follows that the sequence u k ∞ k=0 ⊂ H generated by the preconditioned steepest descent algorithm converges to the unique minimizer u ∈ H. We now wish to estimate the convergence rate.

Estimates of the Convergence Rate for the PSD Method
We summarize some standard results.
Proposition 2.1. Suppose that E satisfies (E1) -(E3). It follows that, for any ν, ξ ∈ H, Proof. Using the orthogonality relation (2.7) and the convexity inequality (2.8), we find Now, suppose d k = 0. Then, by Taylor's theorem, (2.4), and (2.2), Equivalently, we get Suppose that E satisfies (E1) -(E3) and u ∈ H is the unique minimizer of E. Then, for any ξ ∈ H, and, consequently, Proof. This follows immediately from (2.8), because δE[u](ξ) = 0, for all ξ ∈ H. Now, we make the following further assumptions about the pre-conditioner L with respect to the derivatives of the energy E: (L4) There is a constant C 5 > 0 such that for all ν, ξ ∈ H.
(This is the the case if, for example, one .) There is a constant C 6 = C 6 (E 0 ) > 0 such that, for all ν ∈ B, and any ξ ∈ H, Remark 2.4. We note that, practically speaking, (L5) is harder of the last two conditions to enforce. In some sense, if the norm induced by L is not "strong" enough, then there does not exist C 6 > 0 so that (L5) is satisfied.
If the following property holds, we get a simple corollary of the last theorem.
(L6) There is a constant C 8 > 0, such that, for every ν, ξ ∈ H, This implies the strong convexity of E and is, therefore, stronger that (E2).
Corollary 2.6. Suppose that assumptions (E1) -(E3) and (L1) -(L6) are valid. Let u k ∞ k=0 ⊂ H be the sequence generated by (2.5), and define e k := u − u k . Then Proof. By Taylor's theorem and estimate (2.15), we have where θ k is in the line segment from u k to u. The result follows from (2.12).

Notation for Periodic Sobolev Spaces
For the remainder of paper Ω ⊂ R d with d = 2, 3 is a rectangular domain. In what follows, if d = 2 we assume p ∈ [2, ∞); whereas if d = 3 we suppose p ∈ [2,6]. Most of the physically relevant cases correspond to p being an even integer, however, all of our arguments hold for any value of p in the indicated ranges. The Sobolev spaces of periodic functions are defined as follows: where k ∈ N is the differentiability index. Observe that W 0,q per (Ω) =: L q per (Ω) = L q (Ω). We denote the norm of W k,q per (Ω) by · W k,q , or just · L q when k = 0. In the case q = 2 and k = 0, we denote by (·, ·) and · the inner product and norm, respectively. We set H k per (Ω) = W k,2 per (Ω) and immediately remark that, given the range of p, we have H 2 per (Ω) → W 1,p per (Ω). For k ∈ N, the continuous dual of H k per (Ω) is denoted by H −k per (Ω) and their pairing by ·, · k . We set ·, · := ·, · 1 . If L 2 0 (Ω) denotes the set of functions in L 2 (Ω) with zero mean, we define We define a linear operator T :H −1 per (Ω) →H 1 per (Ω) via the following variational problem: given From the Riesz representation theorem it immediately follows that T is well-defined. We define the inner product The induced norm is denoted · H −1 per . The following facts can be easily established [13]: Lemma 3.1. OnH −1 per (Ω) the norm · H −1 per equals the operator norm: for all ζ ∈H −1 per (Ω),

A Fourth-Order Regularized p-Laplacian Problem
We consider the following weak formulation of (1.1): given f ∈ L 2 per (Ω), find u ∈ H 2 per (Ω) such that where 0 < ε ≤ 1 and s > 0 are parameters. Equation (3.1) is mass conservative in the following sense: (u − f , 1) = 0. One can show that the solution of the weak formulation is a minimizer of the following energy: for any ν ∈ H 2 per (Ω), It is not difficult to show that E satisfies (E1) -(E3). The first derivative of E at a point ν may be calculated as follows: for any ξ ∈ H 2 per (Ω), Thus, our original problem is equivalent to the following: find u ∈ H 2 per (Ω), such that, for all ξ ∈ H 2 per (Ω), δE[u](ξ) = 0, which is equivalent to (3.1). This problem has a unique solution, which is, in turn, the unique minimizer of the energy (3.2): The following estimate holds: for all ν, ξ ∈ H 2 per (Ω), The second variation is a continuous bilinear operator. Given a fixed ν ∈ H 2 per (Ω), the action of the second variation on the arbitrary pair (ξ, η) ∈ H 2 per (Ω) × H 2 per (Ω) is given by and we have the bound For this problem we define the pre-conditioner L : Clearly, this is a positive, symmetric operator, and it satisfies assumptions (L1) -(L3), and one can see the similarities with the nonlinear operator in (3.1). We now proceed to find the positive constants for which C 5 , C 6 , C 8 assumptions (L4) -(L6) are satisfied.
Remark 3.2. We could also consider the possibility of changing the metric in the descent direction calculation by, for example, defining the linear operator L k : This is similar to the idea in [26]. The search direction is then found as follows: Our theory does not cover this case, and we will not consider it further here. We plan to examine this in a future work.
Proof. This follows from the Gagliardo-Nirenberg interpolation inequality and elliptic regularity.
Proof. Clearly In addition, the following estimate is available: The simple interpolation inequality in conjunction with Young's inequality yields As a consequence, we get and we conclude that estimate (3.5) is valid by choosing C 5 = min( 1 2 , εs − 1 2 ). Next we derive (3.6). Suppose ν ∈ B. From (3.3) we have Suppose that d = 2. An application of the Sobolev inequality (3.4) in Lemma 3.3 indicates that where Young's inequality is applied in the second step. It follows that, Substituting (3.10) in (3.9) yields We conclude that estimate (3.6) is valid by choosing Note that both C 9 and C 10 are ε and s independent. Following the similar arguments, for p ∈ [2, 6), d = 3, we get For the case p = 6, d = 3, the Sobolev inequality (3.4) degenerates to ∇ξ L p ≤ C 9 ∆ξ , for any ξ ∈ H 2 per (Ω). Hence, we have Therefore, estimate (3.6) is valid by choosing That we can take C 8 = C 5 is the result of a simple calculation that we omit for the sake of brevity. The proof is complete.

A Sixth-Order Regularized p-Laplacian Problem
We now study problem (1.2a) -(1.2b). A weak formulation is given as follows: for f , g ∈ L 2 per (Ω), find u ∈ H 2 per (Ω) and w ∈ H 1 per (Ω) such that where λ ≥ 0, and ε ∈ (0, 1]. This problem is mass-conservative, in the sense that (u − g, 1) = 0, and (w − sλg + f , 1) = 0, and it can be recast as a minimization problem with an energy that involves the H −1 per norm. In particular, for any ν ∈H 2 per (Ω) we define where byḡ we denote the average of g over Ω. Observe that ν − g +ḡ ∈H −1 per , which is required for this energy to make sense. It is straightforward to show that E satisfies (E1) -(E3), with respect to the Hilbert space H =H 2 per (Ω). The first variation of E is given as follows: for any ξ ∈H 2 per (Ω), The unique minimizer of E -let us call it u ∈H 2 per (Ω) for the moment -satisfies δE[u ](ξ) = 0, for all ξ ∈H 2 per (Ω). By the definition of the T operator, there is a unique element w ∈H 1 per (Ω) such that Therefore, we have, for all ξ ∈H 2 per (Ω), Setting u := u +ḡ and w := w + sλḡ −f and using the fact that ξ is of zero mean, we have sλ (u, ξ) + s |∇u| p−2 ∇u, ∇ξ + sε 2 (∆u, ∆ξ) − (w , ξ) = (f , ξ) , ∀ ξ ∈H 2 per .
Similar to the estimate in the fourth-order case (3.3), we have the bound which implies that  Proof. Using integration by parts we get The definition of theH −1 per (Ω) norm implies that Therefore, a combination of (4.7) and (4.8) leads to ∇ξ ≤ ξ which yields the second inequality. The first may be proved in a similar way.
Similar to before, the Gagliardo-Nirenberg inequality, together with elliptic regularity, yield the following interpolation result.   for some C 9 = C 9 (d, p) > 0.
We can now find the coefficients C 5 , C 6 , and C 8 , which establish properties (L4) -(L6) and therefore guarantee the geometric convergence of the PSD method for the sixth-order case.
Proof. The proof is similar to that of Lemma 3.4. Using (3.8) again, we have where the last step is a consequence of the interpolation inequality (4.6): We conclude that estimate (4.10) holds by choosing C 5 = min( 1 3 , ε 4 3 s − 1 3 ). Next we derive (4.11). Inequality (4.3) yields Since ν ∈ B, ∇ν L p ≤ (E(u 0 )) 1 p =: C 10 . Suppose that d = 2. An application of the Sobolev inequality (4.9) from Lemma 4.2 indicates that, for every ξ ∈H 2 per (Ω), where, in the last step, we applied Young's inequality. It follows that, As a result, estimate (4.11) is valid by choosing Similarly, For p ∈ [2, 6), d = 3, we have As a result, estimate (4.11) is valid by choosing For the case p = 6, d = 3, the Sobolev inequality (4.9) degenerates, as before. But it is straightforward to show that estimate (4.11) is valid upon choosing C 6 = 1 + (p − 1)ε −2 C 2 9 C p−2 10 . As before, we omit the simple argument that one may take C 8 = C 5 to satisfy (L6). The proof is complete.

Notation
In this subsection we define the discrete spatial difference operators, function space, inner products and norms, following the notation used in [24,33,36,37,38]. Let Ω = (0, L x ) × (0, L y ), where, for simplicity, we assume L x = L y =: L > 0. We write L = m · h, where m is a positive integer. The parameter h = L m is called the mesh or grid spacing. We define the following two uniform, infinite grids with grid spacing h > 0: Consider the following 2D discrete periodic function spaces: The functions of V per are called vertex centered functions; those of C per are called cell centered functions. The functions of E ew per are called east-west edge-centered functions, and the functions of E ns per are called north-south edge-centered functions. We also define the mean zero space We now define the important difference and average operators on the spaces: with a x , d x : E ew per → C per if is an integer, and a x , d x : V per → E ns per if is a half-integer; and with a y , d y : E ns per → C per if is an integer, and a y , d y : V per → E ew per if is a half-integer. Define the 2D center-to-vertex derivatives D x , D y : C per → V per component-wise as The utility of these definitions is that the differences D x and D y are collocated on the grid, unlike the case for D x , D y . Define the 2D vertex-to-center derivatives d x , d y : V per → C per component-wise as Now the discrete gradient operator, ∇ v h : C per → V per × V per , is defined as The standard 2D discrete Laplacian, ∆ h : C per → C per , is given by The 2D vertex-to-center average, A : V per → C per , is defined to be The 2D skew Laplacian, ∆ v h : C per → C per , is defined as The 2D discrete p-Laplacian operator is defined as . Now we are ready to define the following grid inner products: where summation-by-parts [33,38] guarantees the symmetry and the second equality. We now define the following norms for cell-centered functions. If ν ∈C per , then ν 2 −1 = (ν, ν) −1 . If ν ∈ C per , then ν 2 2 := (ν, ν) 2 ; ν p p := (|ν| p , 1) 2 (1 ≤ p < ∞), and ν ∞ := max 1≤i≤m Similarly, we define the gradient norms: for ν ∈ C per ,
The proof for p = 4, d = 2 can be found in the appendix. Following the similar arguments, the other cases can be proved.
Remark 5.2. Though we have focused on the case d = 2 in this section, we can also define our operators and norms, in particular ∇ v h ξ and ∇ v h ξ p , in three space dimensions. Then for p ∈ [2, 6], we expect for some C 9 = C 9 (d = 3, p) > 0.

Convergence for the Discretized Fourth-Order Problem
The discrete version of (1.1) can be expressed as follows: given f ∈ C per , find u ∈ C per such that This represents a second-order approximation of the solution of (1.1). As in the space continuous case, we formulate an equivalent minimization problem. Using the definitions from subsection 5.1, we have the following discrete energy: given f ∈ C per , for any ν ∈ C per , define This (discrete) energy satisfies (E1) -(E3). The discrete variational derivative at ν ∈ C per is for all ξ ∈ C per , where we have used summation-by-parts [33,38] to obtain the last equality. Given a fixed ν ∈ C per , the action of the second variation on the arbitrary pair (ξ, η) ∈ C per × C per is given by We have the bound: For this problem, we define the pre-conditioner via for all ν, ξ ∈ C per , where we have used summation-by-parts to establish the second equality. In other words, L h [ν] = ν − s∆ h ν + sε 2 ∆ 2 h ν. One will notice the similarity of the pre-conditioner with the nonlinear operator in (5.1). The induced norm is Mimicking the proofs in the continuous case, using summation-by-parts in place of integration-byparts, and Lemma 5.1, we get the following result, whose proof is omitted: For any ν, ξ ∈ C per , For any ν ∈ B and any ξ ∈ C per , we have and C 10 = (pE 0 ) 1 p . We can take C 8 = C 5 to satisfy estimate (2.15) of assumption (L6).

Convergence for the Discretized Sixth-Order Problem
The (second-order accurate) discrete version of (1.2a) -(1.2b) can be expressed as follows: given f , g ∈ C per , find u, w ∈ C per such that As before, it is convenient to switch to the mean-zero version: find u , w ∈C per such that Similar to fourth-order regularized p-Laplacian problem, we define the following discrete energy: for every ν ∈C per For the discrete sixth order problem, we define a linear operator L h :C per →C per via where the second equality may be seen using summation-by-parts [33,38]. This operator satisfies (L1) -(L3), and the next result, which we give without proof for the sake of brevity, shows that (L4) -(L6) are satisfied as well.

Numerical Experiments
In this section we perform some numerical experiments to support the theoretical results. The finite difference search direction equations and Poisson equations are solved efficiently using the Fast Fourier Transform (FFT). We would like to point out that the Fourier pseudo-spectral method can be used to discretize space, and, once again, one can utilize the FFT for the inversion of the linear systems. For descriptions of the pseudo-spectral methods, see, for example, [5,7,21].

Thin Film Epitaxy Model with Slope Selection
In this section we recall the convex splitting numerical scheme in [36] for the thin film epitaxy model with slope selection. Suppose that Ω ⊂ R 2 is a rectangular domain. The energy of an epitaxial thin film is given by where, p ≥ 4 is even, u : Ω → R is the height film, and ε is a constant. The L 2 gradient flow is ∂ t u = −w , w := δE = −∇ · |∇u| p−2 ∇u + ∆u + ε 2 ∆ 2 u, (6.1) and w is called the chemical potential. The model predicts the emergence of a faceted thin film, whose facets have slopes of magnitude approximately one, that coarsens over time. The fully-implicit convex splitting scheme in 2D [36] can be written in operator format as and s > 0 is the time step. Hence, the scheme can be reformulated as the fourth-order problem (5.1) with f = u n − s∆ v h u n and p ≥ 4 and even. In way of summary, to solve N h [u] = f , suppose that iterate u k ∈ C per is given. (Note that k is the PSD solver iteration index, not the time step index, the latter of which we usually denote by n.) We first compute the search direction d k ∈ C per via (2.3): where E h is as defined in (5.2). This equation is efficiently solved using FFT. Once d k is found, we perform a line-search according to (2.6): find α k ∈ R such that q(α k ) = 0, where The approximation sequence is then updated via u k+1 = u k + α k d k . When p = 4 (p = 6), a short calculation shows that q is a cubic (quintic) polynomial whose coefficients can be easily obtained. Moreover, the theory predicts that there is a unique global root for q.

Convergence and complexity of the PSD solver
In this subsection we demonstrate the accuracy and efficiency of the PSD solver by using the epitaxial thin film model with slope selection. We present the results of some convergence tests and perform some sample computations to demonstrate the convergence and near optimal complexity with respect to the grid size h.  To simultaneously demonstrate the spatial accuracy and the efficiency of the solver, we perform a typical time-space convergence test for the fully discrete scheme (6.2) for the slope selection model. As in [33,36], we perform the Cauchy-type convergence test using the following periodic initial data [33]: where Ω = (0, 3.2) 2 . In this test, we compute the Cauchy difference, δ u := u h f (T ) − I f c (u hc (T )), where h c = 2h f , and I f c is a bilinear interpolation operator that maps the coarse grid approximation u hc onto the fine grid. We take a quadratic refinement path, i.e., s = h 2 /10, to equalize the spatial and temporal error contributions. At the final time, T = 0.32, we expect the global error to be O(s) + O(h 2 ) = O(h 2 ) in the discrete · 2 and · ∞ norms, as h, s → 0. The other parameter is given by ε = 0.1. The norms of Cauchy difference, the convergence rates, average iteration number and average CPU time can be found on Table 1. Second-order convergence is observed. At the same time, the average iteration count for the solver remains at around 2. Since we are using a quadratic refinement path, increasing the grid size by a factor of two (decreasing the grid spacing by 2) means increasing the number of time-space degrees of freedom by a factor of 16. But the CPU time increases at a much slower rate. The complexity can be offset, of course, by the fact the starting guesses for the solver at each independent time level are better for smaller time step sizes.  To more directly investigate the complexity of the PSD solver we perform another series of tests to determine the dependences of the convergence rates on ε, h, s, and p, in particular. Consider the following spatially periodic function parametrized by s: u(x, y , s) = 1 2π sin 2πx cos 2πy cos(s). (6.4) First we calculate f := N h [I h (ũ( · , · , s))] ∈ C per , where I h : C 0 per (Ω) → C per is the canonical grid projection operator. Then we compute the sequence u k ∞ k=0 via the PSD algorithm, with the initialization u 0 i,j =ũ(p i , p j , 0) + s 2 sin 4πp i sin 6πp j , hence u k → I h (ũ( · , · , s)), as k → ∞. Define γ k := u k − I h (ũ( · , · , s)) ∞ . We stop the PSD algorithm when γ k ≤ τ := 1 × 10 −8 .
In Figure 2 we plot γ k versus k, on a semi-log scale, for various choices of h, ε, s and p. In  Figure 2(a), the convergence rate (as gleaned from the error reduction) is nearly uniform and nearly independent of h. Figures 2 (b) and (c) indicate that more PSD iterations are required for smaller values of ε and larger values of s, respectively. Figure 2(d) shows that the number of PSD iterations increases with the value of p. These general trends are expected form the theory.  solver. See, for example, [33,36]. The initial data for the simulations are taken as essentially random: where the r i,j are uniformly distributed random numbers in [0, 1]. Time snapshots of the evolution for the epitaxial thin film growth model with p = 4 can be found in Figure 3. The coarsening rates for the p = 4 case are given in Figure 4. These simulation results are consistent with earlier work on this topic in [33,36,39], showing the surface roughness, W , grows like t 1 3 and the energy, E, decays like t − 1 3 . We also present the numerical simulations for the epitaxial thin film growth model with p = 6 in Figure 5. Notice in Figure 5 that the evolution process is significantly different from the process depicted in Figure 3.

Square Phase Field Crystal Model
Suppose that Ω ⊂ R d , d = 2, 3 is a rectangular domain. The energy of square phase field crystal (SPFC) model is given by [16,19,20,30]: where u : Ω → R corresponds to the number density field of the atoms, and ε > 0, γ 0 , γ 1 ≥ 0 are parameters. The SPFC model is the H −1 gradient flow of this energy and is given by We propose the following fully-implicit, nonlinear convex-splitting scheme where g = u n and f = −γ 1 ∆ h u n . Using the techniques of [36,38], we can prove that this scheme is unconditionally energy stable. The fully discrete scheme can also be rewritten in operator format as We can shift the scheme from the affine space of solutions -whose elements ν satisfy (ν − g, 1) 2 = 0to the mean zero space, but this is not necessary for practical implementation. Otherwise, this scheme is in the scope of our theory, and, according to the prescription in Section 5.4, the pre-conditioner should be . Given u k ∈ C per , with u k − g, 1 2 = 0, we compute the search direction d k ∈C per by solving the sixth order linear problem L h [d k ] = f − N h [u k ] using FFT. Once d k is found, we perform the line-search: find α k ∈ R such that q(α k ) = 0, where After this, we update the approximation via u k+1 = u k + α k d k . As before, q is a cubic polynomial (since p = 4) whose coefficients can be precomputed. But this time, two of the coefficients involve the T h = −∆ −1 h operator. Specifically, for q(α) we need to compute where we have use the linearity and symmetry properties of the T h operator. These terms have only to be calculated once per line search, and can be efficiently computed using FFT. In fact, observe that we only need to compute T h d k , at the cost of a single FFT, per line search! The 4-Laplacian term in (6.6) gives preference to rotationally invariant patterns with square symmetry. We perform a simple test showing the emergence of these patterns in this subsection. The initial data for those simulations are similar to (6.5), but we add nucleation sites at specific locations in the domain. The rest of the parameters are given by ε = 1.0; λ = γ 0 = 0.5; γ 1 = 2.0; Ω = (0, 100) 2 ; and s = 0.01. The time snapshots of the evolution by using the given parameters are presented in Figures 6 (one nucleation site) and 7 (four nucleation sites). These tests confirm the emergence of the rotationally invariant square-symmetry patterns in the density field u.

Summary
A preconditioned steepest descent (PSD) solver is proposed and analyzed for fourth and sixth-order regularized p-Laplacian equations. Solution of the highly nonlinear equations are equivalent to the minimizations of the associated convex energies. The energy dissipation property of the PSD solver leads to a bound for the numerical solution at each iteration stage. This fact, coupled with an upperbound for the second derivative of the energy with respect to the metric induced by the pre-conditioner, leads to a geometric convergence rate for our (PSD) solver, which is proved rigorously for both the continuous and discrete space cases. Various numerical results are presented in this article, including a convergence test and a complexity analysis for the PSD solver, as well as long-time simulation results for the thin film epitaxy model with slope selection (both p = 4 and p = 6) and the square phase field crystal model.

A Proof of discrete Sobolev inequality
Herein we only present the proof of (5.1) in Lemma 5.1 with d = 2 and p = 4. The other cases can be handled in the same way. Without loss of generality, we assume that m = N = 2K + 1 is odd and L x = L y = L, so that h = L N = L 2K+1 . We use N, rather than m, for the mesh size, as it is more standard.
For simplicity of presentation, we are focused on the estimate of D x u 4 , and we aim to establish the following estimate: 0 > 0 depends upon L, but is independent of h and u. Due to the periodic boundary conditions for u and its cell-centered representation, it has a corresponding discrete Fourier transformation: Meanwhile, we observe thatû N 0,0 = 0 andf N 0,0 = 0. The first identity comes from the fact that u = 0, while the second one is due to the fact that f = D x u = 0, for any periodic grid function u.
The following preliminary estimates will play a very important role in the later analysis.
Lemma A.1. We have Proof. Parseval's identity (at both the discrete and continuous levels) implies that Based on the fact that hN = L, this in turn results in so that (A.1) is proven.
For the comparison between f = D x u and ∂ x u F , we look at the following Fourier expansions: A comparison of Fourier eigenvalues between |w | and |ν | shows that For the estimate (A.2), we look at similar Fourier expansions: Meanwhile, the derivation of the following inequality is straight forward: in which the eigenvalue estimate (A.6) was used again in the second step. In comparison with (A.5), we arrive at (A.4). The proof of Lemma A.1 is complete.
The following lemma gives a bound of the discrete · 4 norm of the grid function f , in terms of the continuous L 4 norm of its continuous version f F . Proof. We denote the following grid function g i+1/2,j = f i+1/2,j 2 .
A direct calculation shows that f 4 = ( g 2 ) 1 2 . (A.14) Note that both norms are discrete in the above identity. Moreover, we assume the grid function g has a discrete Fourier expansion as With an application of the Parseval equality at both the discrete and continuous levels, we have On the other hand, we also denote H(x, y ) = (f F (x, y )) 2 =
The reason for H ∈ P 2K is because f F ∈ P K . We note that H = G, since H ∈ P 2K , while G ∈ P K , although H and G have the same interpolation values on at the numerical grid points (x i , y j+1/2 ). In other words, g is the interpolation of H onto the numerical grid point and G is the continuous version of g in P K . As a result, collocation coefficientsĝ N c for G are not equal toĥ N for H, due to the aliasing error. In more detail, for −K ≤ , m ≤ K, we have the following representations: Therefore, a combination of (A.14), (A.17) and (A. 18) results in This finishes the proof of (A.13).
Now we proceed into the proof of Proposition 5.1.
Proof. We begin with an application of (A.13) in Lemma A.2: Meanwhile, using the fact that f F = 0, we apply the 2-D Sobolev inequality and get ) where the following elliptic regularity estimate is applied:   As a consequence, the first case of (5.1) (with d = 2, p = 4) is valid, by setting C 0 = √ 2C (1) 0 . The other cases could be analyzed in the same way. This finishes the proof of Proposition 5.1.