STABILITY AND DISCRETIZATION ERROR ANALYSIS FOR THE CAHN–HILLIARD SYSTEM VIA RELATIVE ENERGY ESTIMATES

. The stability of solutions to the Cahn–Hilliard equation with concentration dependent mobility with respect to perturbations is studied by means of relative energy estimates. As a by-product of this analysis, a weak-strong uniqueness principle is derived on the continuous level under realistic regularity assumptions on strong solutions. The stability estimates are further inherited almost verbatim by appropriate Galerkin approximations in space and time. This allows to derive sharp bounds for the discretization error in terms of certain projection errors and to establish order-optimal a priori error estimates for semi-and fully discrete approximation schemes. Numerical tests are presented for illustration of the theoretical results.

The second equation (2) defines the chemical potential  =   ℰ() as the variational derivative of the associated energy functional which together with (1) encodes the gradient flow structure of the problem and implies the decay of the energy ℰ() along weak solutions.This ensures thermodynamic consistency of the model and allows to establish existence of weak solutions by Galerkin approximation, energy estimates, and compactness arguments.On this basis, existence and regularity of weak solutions for the Cahn-Hilliard equation has been established in [18] for constant mobility and polynomial potential.Logarithmic potentials and concentration dependent mobilities were treated in [5,15]; we refer to [3,7,26] for further results concerning the extension to multi-component systems and multiphysical problems.In recent work [11,34,37], we considered logarithmic potentials and concentration dependent mobility functions in the context of models for viscoelastic phase separation.
Many theoretical results about the Cahn-Hilliard equation are based on certain approximations in space and/or time and appropriate energy estimates.Finite element approximations of the fourth-order system resulting after elimination of the chemical potential were analyzed in [18].A mixed finite element approximation for constant mobilities treating  and  as separate variables was proposed in [21] and further analyzed in [17,20].For extensions to logarithmic potentials and degenerate mobilities, we again refer to [5,15].In [27,28], the analysis of finite element approximations has been extended to study the thin-interface limit  → 0. Apart from finite element methods, alternative discretization schemes, like discontinuous Galerkin methods [31,33,40] and Fourier-spectral approximations [32] have been investigated.Extensive research has further been devoted to developing stable second order approximations in time; see [38] for an overview and comparison of different approaches.In a recent paper [16], which is probably closest to our investigations, an unconditionally well-posed fully discrete two-step approximation was proposed and a full error analysis was presented yielding order optimal convergence rates.Quantitative convergence results in most works are derived for the case of constant mobility, which allows to apply arguments of linear theory and to cover the terms stemming from the nonlinearity of the chemical potential  by perturbation arguments.
In this paper, we consider problems with concentration dependent mobility, and we utilize genuinely nonlinear arguments, namely relative energy estimates, to conduct a quantitative stability and error analysis.Related techniques have been used intensively for the study of nonlinear evolution problems and, more recently, also for the convergence analysis of corresponding discretization methods; see [24,30] for an introduction to the topic and [23,25,29,35] for convergence and asymptotic analysis for fluid flow problems via relative energy estimates.For discretization, we here consider conforming finite element approximations of second order in space and time, but our arguments, in principle, also apply to higher order approximations and inexact Galerkin approximations.
We here study the case of nonlinear but non-degenerate parameters, e.g., strictly positive mobility () and polynomially bounded double well potential  (), for which one can guarantee existence of smooth solutions which is required for the convergence rate analysis.If one assumes the existence of a smooth solution bounded away from the pure states, these conditions could be further relaxed.
The first basic result of this paper is a formal relative energy estimate which allows to deduce quantitative perturbation bounds for sufficiently regular solutions of (1) and (2).As a by-product of this analysis, we also obtain a weak-strong uniqueness principle and thus a rather general argument for uniqueness.Due to the variational character, these stability estimates are inherited almost verbatim by Galerkin approximations in space and Petrov-Galerkin approximation in time, which is our basic approach towards a systematic error analysis.The structure of the relative energy estimates further provides guidelines for the choice of appropriate projection operators required in the error analysis.The discrete relative energy estimates then allow to estimate the discretization error by more or less standard projection error estimates, which finally leads to optimal convergence rates under minimal and less restrictive smoothness requirements than in previous works.This nonlinear convergence rate analysis can be seen as the main contribution of our manuscript.
Compared to previous works, the main novelty of our approach consists of (a) the development of a genuinely nonlinear stability and error analysis via relative energy estimates, which can be extended to more complex systems; (b) a concise and rather sharp error analysis of a variational discretization method with order optimal convergence rates under minimal regularity assumptions.
The remainder of the manuscript is organized as follows: In Section 2, we introduce our notation and basic assumptions and recall some results about existence and regularity of solutions.In Section 3, we introduce the relative energy functional and present a formal relative energy estimate which serves as the basis for the following considerations.Furthermore, we will also deduce the weak-strong uniqueness principle in course of the analysis.In Section 4, we study the semi-discretization in space by a mixed finite element method.We will see that the relative energy estimate translates almost verbatim to the semi-discrete setting.This allows us to estimate the difference between the semi-discrete solution and a particular projection of the continuous solution by projection errors and to derive order optimal error convergence rates.In Section 5, we then consider the time discretization by a Petrov-Galerkin approximation, which allows us to extend our arguments almost verbatim to the fully discrete setting.For illustration of our theoretical results, we present some preliminary numerical results in Section 6.Some technical details are provided in the appendix.
Remark 1.Here we consider the non-degenerate case with polynomially bounded double well potential  () and strictly positive mobility function ().These assumptions allow to prove existence and uniqueness of smooth solutions which is required to establish quantitative convergence rates below.Even in the case of more general potentials and mobilities, we would expect smooth solutions to stay away from the pure states, i.e., from the possible singularities or degeneracies of  (•) and (•).Validity of the assumptions may then be established by suitable regularization, which is used for the proof of existence in the degenerate case; see e.g., [6,19].
Under the above assumptions, the existence of periodic weak solutions can be deduced from classical results.For later reference, we make a precise statement.
Lemma 2. Let (A1)-(A3) hold.Then for any  0 ∈  1  (Ω), there exists at least one periodic weak solution (, ) of problem (1) and (2) with (0) =  0 and satisfying with   (‖ 0 ‖  ) depending only on the bounds for the coefficients, the domain Ω, the time horizon  , and the bounds for the initial value.In dimension  = 3 the estimates  > 1 are only valid for sufficiently small  .
Proof.Existence of weak solutions and the a priori bounds for  = 1 and any time  > 0 in dimension  = 2, 3 follow from standard arguments; see [4,6] for similar results under even more general assumptions on the problem data.Conservation of mass and dissipation of energy follow from the variational identities ( 7) and ( 8) by formally testing with (, ) = (1, 0) and (, ) = (,   ), respectively.Improved regularity and the bounds with  > 1, which require a restriction on the maximal time  in dimension  = 3, can be obtained by a boot-strap argument and regularity results for the Poisson problem; details are given in Appendix C [8].
Remark 3. From the estimates of Lemma 2 and the embedding theorem for Bochner spaces, see e.g., Chapter 25 of [39], one can see that weak solutions (, ) are continuous functions of time, e.g.,  0 ∈  3  (Ω), we have and hence  is uniformly bounded on Ω × (0,  ).This will be used in Section 4 below.

Stability via relative energy
To measure the difference between a solution (, ) of ( 7) and ( 8) and a solution ( φ, μ) of the perturbed problem (10) and (11), we utilize a regularized relative energy functional with regularization parameter where − 1 is the constant in the lower bound for  ′′ from assumption (A3).This implies that the regularized energy functional ℰ  () = ℰ() +  2 ‖‖ 2 0 becomes strictly convex such that ℰ  (| φ) amounts to the associated Bregman divergence [9].Moreover, the norm distance of two functions can be bounded by the relative energy.
Proof.For ease of presentation, we assume for the moment that (, ) and ( φ, μ) are sufficiently regular, such that all computations in the following are justified.The general case can then be deduced by a density argument; details are given in the appendix.Formal differentiation of the relative energy with respect to time yields From the definition of the relative energy ℰ  and using the variational identities ( 7), ( 8) and ( 10), (11), which are satisfied by the functions (, ) and ( φ, μ), we get We can now estimate the individual terms separately.Before we proceed, let us note that by the energy bounds for weak solutions (, ), Lemma 2, and by the assumptions on φ in the statement of Theorem 6, we know that Using Hölder's and Young's inequalities, we can then bound with  > 0 arbitrary, constant (,  2 ) =  2 /(4), and  2 denoting the upper bound for the function  in assumption (A2).By definition of the dual norm, a Poincaré inequality, and the bounds for the coefficients, the second term can be estimated by In the last step, we utilized Young's inequality to separate the factors with the same parameter  > 0 as before.
For the second term on the right-hand side, we can use the identities ( 8) and ( 11) with  = 1, which leads to From the bounds for the potential  in assumption (A3), we can further deduce that An application of Hölder's inequality, the norm estimates for the continuous embedding of  1 (Ω) into   (Ω), and the uniform bounds for , φ in (15), then lead to 2 + 2 (2) Using ‖r 2 ‖ 0,1 ≤ (Ω)‖r 2 ‖ 1 and the lower bound (13) for the relative energy, we get 2 , 3 Condition (13) further allows us to estimate From the bounds in assumption (A3), we can deduce that Using Hölders inequality, embedding estimates, and the uniform bounds in (15), we can further bound the fifth term in the above estimate by 2 + 2 , By combination of the individual estimates and choosing  = 1/6, we finally obtain with constants  0 ,  1 ,  2 ,  3 depending only on the bounds for the coefficients, the domain, and the bounds for ‖ 0 ‖ 1 and ‖ φ‖  ∞ ( 1 ) .An application of Gronwall's inequality (A.1) with , then leads to the stability estimate of the theorem with  =  0  +  1 Ĉ and  = max{ 2 ,  3 }.
Remark 7. The lower bound (13) for the relative energy, and the estimate for the relative dissipation immediately lead to the bounds for the error.With similar arguments as used for the estimate of the term (ii), we can also bound the full norm ‖ − μ‖  2 ( 1  ) .The stability estimate thus provides perturbation bounds in the natural norms for the problem.
Remark 9. Note that for regular initial values, e.g., φ(0) =  0 ∈  2  (Ω), the existence of a weak solution ( φ, μ) with the required regularity follows from Lemma 2. In that case, we therefore obtain a unique weak solution.

Galerkin semi-discretization
We now turn to the systematic discretization in space, for which we consider a Galerkin approximation of the variational principle (7) and (8) with second order conforming finite elements.As will become clear from our analysis, higher order and non-conforming approximations could be treated with similar arguments.
Let  ℎ be geometrically conforming partition of Ω ⊂ R  ,  = 2, 3 into triangles or tetrahedra [14,22].We denote by   and ℎ  the inner-circle radius and diameter of the element  ∈  ℎ and call ℎ = max ∈ ℎ ℎ  the global mesh size.We assume that  ℎ is quasi-uniform, i.e., there exists a constant  > 0 such that ℎ ≤   ≤ ℎ  ≤ ℎ for all  ∈  ℎ , and that  ℎ is periodic in the sense that it can be extended periodically to periodic extensions of the domain Ω.We then denote by the space of continuous periodic piecewise quadratic polynomials over the mesh  ℎ .We further introduce the approximation spaces W ℎ (0,  ) :=  1 (0,  ;  ℎ ) and Q ℎ (0,  ) :=  2 (0,  ;  ℎ ).
The semi-discrete approximation for ( 7) and ( 8) now reads as follows.

Semi-discrete stability estimate
With similar arguments as on the continuous level, we now establish stability of the semi-discrete solution with respect to perturbations.For a given pair of functions ( φℎ , μℎ ) ∈ W ℎ (0,  ) × Q ℎ (0,  ), we define semi-discrete residuals (r 1,ℎ , r2,ℎ ) ∈  2 (0,  ;  ℎ ×  ℎ ) by the variational identities for all  ℎ ,  ℎ ∈  ℎ and 0 ≤  ≤  .The functions ( φℎ , μℎ ) can again be understood as solutions of a perturbed semi-discrete problem.With almost identical arguments as used in the proof of Theorem 6, we now obtain the following stability estimate.
The assertion follows with the very same arguments as used in the proof of Theorem 6; details are left to the reader.Lemma 12 allows to investigate the stability of the semi-discrete solution ( ℎ ,  ℎ ) with respect to perturbations in the problem data.By choosing ( φℎ , μℎ ) as a particular discrete approximation for the solution (, ) of ( 1) and ( 2), we will be able to derive quantitative error estimates for the semi-discrete approximation.

Auxiliary results
We start by introducing some projection operators and recall the corresponding error estimates.Let  0 ℎ :  1  (Ω) →  ℎ denote the  2 -orthogonal projection which is characterized by By definition,  0 ℎ is a contraction in  2 (Ω) and on quasi-uniform meshes for all −1 ≤  ≤  and 0 ≤  ≤ 3; see [10].We further make use of the  1 -elliptic projection  1 ℎ :  1  (Ω) →  ℎ , which is characterized by the variational problem By standard finite element error analysis and duality arguments, one can show that for all −1 ≤  ≤  and 1 ≤  ≤ 3; see again [10] for details.Since we assumed quasi-uniformity of the simplicial mesh  ℎ , we can further use inverse inequalities which hold for all functions  ℎ ∈  ℎ and all 1 ≤  ≤  ≤ ∞ in dimension .By combining the previous estimates, one can further see that in dimension  ≤ 3. Let us note that all estimates also hold in dimension one, i.e., for piecewise polynomial approximations in time.

Projection error estimates
Let (, ) be a periodic weak solution of the system (1) and (2).We then define φℎ () =  1 ℎ () ∈  ℎ , 0 ≤  ≤  , as the  1 -elliptic projection, and μℎ () ∈  ℎ by solving the elliptic variational problems for all  ℎ ∈  ℎ and all 0 ≤  ≤  .This problem is linear in μℎ and finite-dimensional, hence existence of a unique solution follows directly.For this choice of approximations ( φℎ , μℎ ), we have the following error estimates.
Proof.The estimates for  − φℎ and    −   φℎ follow directly from (26).We then use the triangle inequality to split the error in the chemical potential into With the help of ( 24) the last term can be estimated by ‖ 0 ℎ  − ‖ 1 ≤ ℎ 2 ‖‖ 3 .By the inverse inequalities (27), the discrete error can be bounded by and for the error in the  2 -norm, we can deduce from (23) that since  ℎ = μℎ −  0 ℎ  ∈  ℎ .We can then use (29) with this test function  ℎ , to get where we used the particular choice of φℎ =  1 ℎ  and ( 25), to replace the gradient term in the second step.Proceeding with standard arguments, we then obtain For the nonlinear term, we here used the mean value theorem and the polynomial bounds for  ′′ as well as ‖‖ 0,∞ + ‖ φℎ ‖ 0,∞ ≤ ‖‖ 2 .In summary, we thus obtain with constant  independent of the mesh size and uniform for all 0 ≤  ≤  .

Error estimates
Using that (, ) solves ( 7) and ( 8) and by the definition of ( φℎ , μℎ ), one can see that (19) and ( 20) is satisfied with residuals r2,ℎ = 0 and By the properties of the discrete dual norm ‖•‖ −1,ℎ and standard approximation error estimates, see Lemma 13, the residual r1,ℎ can further be bounded by with constants ,  ′ depending only on bounds on the coefficients, the domain Ω, the mesh regularity, and the constant   (‖ 0 ‖ 3 ) for the solution in Lemma 2. We can now use Lemma 12 to obtain the following bounds for the discrete error.
By combination of the previous estimates we now immediately obtain the following error bounds for the Galerkin semi-discretization with quadratic finite elements.
Remark 16.The convergence rates in the theorem are optimal with respect to the approximation properties of quadratic finite elements.Moreover, the regularity assumption on the initial value is necessary to establish the predicted convergence rates.The convergence result is both, order optimal and sharp, i.e., obtained under minimal smoothness assumptions on the problem data.The constant  ′  depends exponentially on  and  −1 .
Let us emphasize that functions in W ℎ, (0,  ) are continuous in time and piecewise linear, while functions qℎ, ∈ Q ℎ, are piecewise constant in time, which is designated by the bar symbol.The discrete approximation for (1) and ( 2) then is the following.
To show existence, we use an induction argument.Let  ℎ, ( −1 ) be given.Then in the th time step, only the function values   ℎ :=  ℎ, (  ) and  −1/2 ℎ := μℎ, (  −  /2) are to be determined.From the discrete energydissipation identity, the bounds for the coefficients, and the equivalence of norms on finite dimensional spaces, one can deduce that potential solutions are necessarily bounded.Existence of a solution for the th time step then follows from Brouwer's fixed-point theorem.The uniform bounds for the solution, finally, follow directly from the energy-dissipation identity and using ( 13) and ( 16).
Remark 20.Uniqueness of the discrete solution can be shown under a mild restriction  ≤  0 (ℎ) on the time step size.In Section 5.4 below, we will show that uniqueness holds for  ≤ ℎ  with some  ≤ 1, if the solution (, ) is sufficiently regular.The choice  = ℎ, which seems reasonable in view of the estimates of Theorem 27, therefore will lead to unique solutions for the fully discrete problem.
In the following, we first establish a discrete analogue of the stability estimate derive in Theorem 6, and then derive convergence rates for the fully-discrete scheme.
Proof.By the fundamental theorem of calculus, we obtain In the last step, we use the identities ( 33) and ( 36) with wℎ, =    ℎ, ∈ Q ℎ, .Since the function    ℎ, −   φℎ, is piecewise constant in time, we can replace At this point, we can start to estimate the individual terms in the same manner, as in the proof of Theorem 6.

Auxiliary results
Similar to the semi-discrete case, we use certain projections to define suitable approximations φℎ, and μℎ, for solutions (, ) to ( 1) and ( 2) that allow us to take advantage of the discrete stability estimate.So, let denote the piecewise linear interpolation with respect to time.Furthermore, let be the  2 -orthogonal projection to piecewise constant functions in time.For later reference, we summarize some important properties of these operators.
Lemma 23.For  ∈  1, (0,  ), 1 ≤ ,  ≤ ∞, there holds and for  ∈  , (0,  ) with 1 ≤  ≤ 2 and 1 ≤  ≤  ≤ ∞, one has Moreover, interpolation and projection commutes with differentiation, i.e., The proof for these standard results can be found, e.g., in [10].The interpolation operator naturally extends to vector valued functions, and we use the same symbol in that case.For the piecewise constant  2 -projection we can show the following estimate for the product error; see Appendix B for a proof.
Lemma 24.Let ,  ∈  2, (0,  ) and ā = π0 denote the  2 -orthogonal projection onto piecewise constants.Then with a constant  independent of  and  as well as the functions  and .
As a next step, we derive bounds for the discrete residuals in terms of interpolation and projection errors.For ease of notation, we will write  , () =  , (, ; ) for different choices of the time interval (, ), which will be clear from the context.Lemma 26.Let (, ) be a sufficiently regular periodic solution of (7) and (8).Then for all 0 <   ≤  with constants (•) independent of ℎ,  , and   .
Proof.Since vℎ, is piecewise constant in time, we can use (39), the definition of φℎ, , and the bounds for the  1 -projection error, to estimate the residual by Here and in the following, we use  = π0   to abbreviate the projection onto piecewise constant functions in time and all Bochner-norms refer to the time interval ( −1 ,   ) under consideration.The remaining term can be further estimated by Using the boundedness of , the definition of μℎ, , and the stability and error estimates for the  2 -projection (24), we immediately obtain For the second term, we use a triangle inequality, the error bounds (26) for the  1 -projection, the interpolation error estimate (38), and the lower bound (13) for the relative energy.In summary, this leads to For the third term, we observe that this is a second order approximation on the midpoint of the time interval and using the estimate (40) we obtain To estimate  2 ( 2 ) norm of ()∇, we differentiate ()∇ in time and use triangle inequalities, the uniform bounds for (•), and the Hlder inequalities in space and time.The critical terms are the higher order terms ()  ∇,  ′ ()∇  , which are estimated in  2 ( 2 ).Indeed, they can be bounded by ‖‖  2 ( 1 ) and . By combination of the previous estimates, we thus obtain ) ) independent of ℎ and  .Before turning to the bound for the second residual, let us observe that which follows from the definition of φℎ, and the variational characterization (25) of  1 ℎ .The second residual can then be expressed equivalently in strong form as r2,ℎ, = ( where  = π0   denotes the piecewise constant projection of  with respect to time.This point wise representation allows us to estimate the second residual by We again estimate the individual terms separately.For the first, we use the contraction property of the  2projection  0 ℎ and the interpolation error estimate (38) to obtain By the estimate (26) for the  1 -projection  1 ℎ to obtain for the second term For the third term, we use that  and its discrete counterpart φℎ, =  1   1 ℎ  can be uniformly bounded in  ∞ (0,  ;  1,∞  (Ω)).Therefore, all terms  () (•) appearing in the following can be bounded uniformly by a constant ( ).This leads to ) .A quick inspection of the last term shows that its evaluation involves up to cubic products of  and its derivatives, with the highest order terms given by  2   ∇, (  ) 2 ∇, and   ∇  , respectively.This allows to establish the following bounds ) ) 3 .In summary, the second residual can thus be bounded by with the two solution dependent constants given by  2 (, ) independent of ℎ and  .

Error estimates
Together with the discrete stability estimate of Lemma 21 and a Gronwall-type argument, similar as already used in the proof of that result, we can now obtain the following convergence rate estimates.
Proof.We may proceed almost verbatim to the proof of Lemma 21 and insert the above estimates for the residual terms, to see that bounded uniformly in time.The proof of the assertion then follows in the same manner as that of Lemma 21.Note that it suffices to consider the case that  ≤  0 is sufficiently small, since for large  the result already follows from the a priori estimates.

Uniqueness of the fully discrete solution
As a final step of our analysis, we now establish uniqueness of the discrete solution under a mild restriction on the time step size; see Remark 20.Let us start with the observation that under the conditions of the previous theorem, μℎ, is uniformly bounded in the norm of  ∞ ( 1,3   ).To see this, note that The last two terms are uniformly bounded by assumption and standard projection error estimates.For the first term on the right-hand side, we use the second inverse inequality in (27) with  = 3,  = 2,  ≤ 3 in space and with  = ∞,  = 2,  = 1 in time, as well as the convergence estimates of the previous theorem, to see that For any choice ℎ 3 ≤  ≤ ℎ 1/3 , one can thus conclude that ‖μ ℎ, ‖  ∞ ( 1,3  ) ≤  ′ .Uniqueness can now be deduced as follows: Let the assumptions of Theorem 27 be valid.Furthermore, let ( ℎ, , μℎ, ) ∈ W ℎ, (0,  ) × Q ℎ, (0,  ) and ( φℎ, , μℎ, ) ∈ W ℎ, (0,  ) × Q ℎ, (0,  ) denote two solutions of Problem 17 with the same initial data  ℎ, (0) = φℎ, (0) and with time step size ℎ 3 ≤  ≤ ℎ 1/3 and  ≤  0 sufficiently small.Then the residuals defined by ( 35) and ( 36) are r2,ℎ, = 0 and for all vℎ, ∈ Π 0 ( −1 ,   ;  ℎ ).Using the bounds for the coefficients and ( 13), the residual term r1,ℎ, can be further estimated by The last term can be treated by a Gronwall estimate, similar as in the proof of Lemma 21 and Theorem 27.
Remark 28.A brief inspection of the arguments reveal, that the regularity assumptions on the true solution could be somewhat relaxed, which will however lead to tighter bounds ℎ  ≤  ≤ ℎ 1/ with 1 ≤  ≤ 3 for the admissible time step sizes.The choice  = ℎ, seems reasonable and leads to a uniqueness result under minimal regularity assumptions.If the mobility function () ≡  is independent of the concentration, then the above considerations become obsolete, since the relevant terms in the stability estimate vanish.

Numerical validation
The aim of this section is to illustrate our theoretical results, in particular, the convergence rate estimates obtained in Theorems 15 and 27.We solve a typical test problem, which is specified as follows: The interface parameter  is set to  = 0.003, the potential and mobility functions are defined as  () = 0.3( − 0.99) 2 ( + 0.99) 2 and () = (1 − ) 2 (1 + ) 2 + 10 −3 .The computational domain Ω is a unit square, Ω = (0, 1) 2 .The system (1) and ( 2) is complemented by periodic boundary conditions and following initial data for the phase fraction   0 (, ) = 0.1 sin(4) sin(2) + 0.6.For all our computations, we use the fully discrete approximation of Problem 17 on a sequence of uniformly refined triangulations  ℎ in space and equidistant grids ℐ  in time.The arising nonlinear systems at each discrete time step are solved by the Newton method with absolute residual tolerance 10 −12 .Since in our test problem the nonlinearities are polynomial, all integrals appearing in the discrete scheme can be evaluated exactly.For more general nonlinearities numerical integration can be used.Integration errors and iteration errors of the Newton method can be analyzed with slight modifications of the above-mentioned proofs.
Time snapshot of the phase fraction  ℎ, computed by our Petrov-Galerkin method are depicted in Figure 1.One can clearly observe the expected evolution from a rather uniform distribution to an almost completely separated configuration.As predicted by our theoretical results, the solution remains smooth over the whole time interval.
We now turn to the validation of the convergence rates.The parameter functions and initial conditions are chose as described above and the final time is set to  = 0.76.Since no analytical solution is available, the discretization error is estimated by comparing the computed solutions ( ℎ, , μℎ, ) with those computed on uniformly refined grids in space and time.In accordance with our convergence analysis, the error quantities for the fully-discrete scheme are thus defined by In order to evaluate the convergence rates of the semi-discrete scheme, we choose a very small step size  * , and define the corresponding error quantities as Table 1.Errors and convergence rates for the semi-discrete and fully-discrete approximations.In Table 1, we summarize the results of our computations obtained on a sequence of uniformly refined meshes with mesh size ℎ  = 2 −(3+) ,  = 0, . . ., 5 and time steps   = 0.16 • ℎ  .For the results concerning the semidiscretization, the time step is chosen  * = 0.16 • 2 −9 .Since nested grids are used in all our computations, the error quantities defined above can be computed exactly.
The experimental order of convergence (eoc) is computed by comparing the errors of two consecutive refinements.In perfect agreement with the theoretical predictions of Theorems 15 and 27, we observe second order convergence for the total errors in space and time.The same convergence rates are obtained for each error term separately that arises in the definition of the discrete error measures  ℎ respectively  ℎ, .

Discussion
In this paper, we studied the stability, regularity, and uniqueness of solutions to the Cahn-Hilliard equation with concentration-dependent mobility.The variational characterization of weak solutions and relative energy estimates were used as the main ingredients of our analysis, and the latter greatly simplified the handling of nonlinear terms in the problem.The basic tools of our analysis are applicable almost verbatim to discretization schemes based on variational principles, i.e., Galerkin finite-element approximations in space and Petrov-Galerkin approximation in time.The variational time discretization, which is tightly related to the average vector field methods, leads to fully-implicit schemes which, however, can be solved efficiently by Newton-iterations, and which allows for a structured and transparent error analysis.The convergence results obtained in the paper are of optimal order and the result for the semi-discretization is sharp concerning regularity requirements of the solution.Some additional regularity is required for the fully-discrete scheme, which can be explained by the lack of strong stability of the Petrov-Galerkin time discretization; see [2] for details.In principle, the proposed schemes can be extended immediately to higher order in space and time.Further investigations in this direction and the extension to more complex multiphase problems, e.g., the Cahn-Hilliard Navier-Stokes equations, will be topics of future research.
with constants ,  for all , , only depending on the uniform bounds Hence, the constants ,  in the above estimate can be chosen independent of , .Using the strong convergence of the residuals in the corresponding norms, we may pass to the limit in the integral on the right-hand side.Furthermore, application of Egorov's Theorem yields almost everywhere convergence of (  ).With this and standard weak convergence results using Fatou's Lemma yields which allows us to pass to the limit in the relative dissipation term.A similar process is used to derive the energy inequality for the standard weak solution of ( 1) and (2).By the continuous embedding of  (0,  ) into ([0,  ];  1  (Ω), we obtain convergence of the energy ℰ  (  ()| φ ()) → ℰ(| φ) with ,  → ∞, and in summary, we thus obtain (14).

Figure 1 .
Figure 1.Snapshots of the phase fraction  ℎ, and time evolution of the energy.