title: Efficient second-order semi-implicit finite element method for fourth-order nonlinear diffusion equations
authors: Sana Keita Abdelaziz Beljadid Yves Bourgault
text: 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  addr1]Sana Keita     sana.keita@um6p.ma
  [corres]Corresponding author.
  addr1,addr2]Abdelaziz Beljadidcorres abdelaziz.beljadid@um6p.ma
  addr2]Yves Bourgault     ybourg@uottawa.ca
  [addr1]International Water Research Institute, Mohammed VI Polytechnic University,  Morocco 
  [addr2]Department of Mathematics and Statistics, University of Ottawa, Canada                 
   
  
  
  
  We focus here on a class of fourth-order parabolic equations that can be written as a system of second-order equations by introducing an auxiliary variable. We design a novel second-order fully discrete mixed finite element method to approximate these equations. In our approach, we propose new techniques using the second-order backward differentiation formula for the time derivative and a special technique for the approximation of nonlinear terms. The use of the proposed technique for nonlinear terms makes the developed numerical scheme efficient in terms of computational cost since the proposed method only deals with a linear system at each time step and no iterative resolution is needed. A numerical convergence study is performed using the method of manufactured and analytical solutions of the system where we investigate different boundary conditions. With respect to the spatial discretization, convergence rates are found to at least match a priori error estimates available for linear problems. The convergence analysis is completed with an investigation of the temporal discretization where we numerically demonstrate the second-order time-accuracy of the proposed scheme using the method of reference solution. We present a series of numerical tests to demonstrate the efficiency and robustness of the proposed scheme. 
  
  
  
   Second-order scheme Numerical convergence analysis Mixed finite element Fourth-order equation  Diffuse-interface model Cahn-Hilliard equation
  
  
  
  
  
  
  
  § INTRODUCTION
  
  High order partial differential equations (PDEs) containing the bi-harmonic operator often arise in two-phase and multiphase problems from fluid mechanics and solid mechanics. Examples in the literature include the thin film equations <cit.>, the Cahn-Hilliard (CH) equations <cit.>, ocean and atmosphere circulation models <cit.>, tumor-growth models <cit.>, cancer modeling <cit.>, biological entities <cit.>, image processing <cit.>, the fourth-order model for the approximation of crack propagation <cit.>,  models of unstable infiltration through porous media <cit.>,  etc. 
  
  Nonlinearity and discretization of high-order spatial operators (fourth-order) are the main issues encountered when numerically solving these equations. The nonlinearity is often induced by non-constant mobility and nonlinear free energy functions. As a result, several different procedures have been employed over the years to deal with higher-order equations, for instance, mixed formulations, where the governing equations are split into a coupled system of lower-order differential equations <cit.>.
  
  There have been many studies on numerical methods for solving parabolic equations involving the bi-harmonic operator, including finite difference <cit.>, finite volume <cit.>, Fourier spectral <cit.> and finite element methods <cit.>. Most studies using finite element method focus on mixed formulation to approximate the CH equations <cit.>, see also the references therein. Second-order time-accurate methods for the CH equations have mostly been studied for variants of the Crank-Nicolson scheme with various treatment of the nonlinear energy term <cit.>. Convex-concave splitting of the energy term is often employed, but it leads to nonlinear schemes except for some particular forms of the energy function. Addition of stabilization terms, use of Hermite quadrature formula and/or explicit treatment or extrapolation of the energy function, although leading to a linear scheme, is very limited in term of stability constraint <cit.>. Recently, Yang et al. <cit.> proposed a family of second-order unconditionally energy-stable schemes for CH equations by reformulating the equations into an equivalent system employing a scalar auxiliary variable, giving rise  to four equations to solve at each time step.
  
  In the current study, we focus on the numerical approximation of the CH equations. We propose a novel linear semi-implicit time stepping scheme to approximate the  problem. The proposed scheme maintains second-order accuracy in time by using second-order backward difference formula (BDF2) for the time derivative, extrapolation formula for the non-constant mobility function and a special technique using semi-implicit Taylor approximation for the nonlinear energy function. The use of Taylor expansion on a term in solving partial differential equations was introduced by Rosenbrock <cit.>. A scheme employing the Crank-Nicolson time integration and Taylor series expansions of the nonlinear terms were proposed by Vignal et al. <cit.> for general classes of phase-field models with polynomial potentials.  Their resulting system is nonlinear and therefore, requires an iterative method for the linearization. Our scheme applies to general classes of phase-field models with any form of free energy function of class 𝒞^2 and does not require an iterative method. A linear system is solved at each time step. We perform numerical convergence analysis in space using both a manufactured and analytical solutions to demonstrate optimal convergence rates of the L^2- and H^1-error norms using standard Lagrange elements. Better, super-convergence of the L^2-error norm is observed with quadratic Lagrange elements. Numerical results from selected test cases also show that the proposed scheme is efficient in terms of computational cost and  satisfies the free energy dissipation and the mass conservation properties.
  
  The outline of the paper is as follows.  In section <ref>, the model equation is introduced and its main mathematical properties are briefly discussed. In section <ref>, the fully discrete finite element approximation scheme is presented. In Section <ref>, the spatial convergence rates, accuracy and robustness of the proposed scheme are numerically investigated by means of several numerical test cases. Second-order time-accuracy is verified with both constant and variable mobility function, and the discrete energy dissipation and the discrete mass conservation properties of the scheme are illustrated. Section <ref> provides a summary and some concluding remarks.
  
  
  
  § MODEL EQUATION
  
  
  Let Ω⊂ℝ^d, where d=2 or 3, be an open bounded set with sufficiently smooth boundary ∂Ω, and T>0 a fixed positive time. We focus on a class of fourth-order parabolic PDEs of the form
  
      ∂_tu=-γ∇·(f(u)∇Δ u) + ∇·(f(u)∇φ'(u))   in Ω×(0,T),
  
  with unknown u:Ω×[0,T]→ℝ, the given functions f:ℝ→ℝ_⩾ 0, φ:ℝ→ℝ are  smooth, and γ∈ℝ_⩾ 0 is a constant. Depending on the applications, u can describe the height of a liquid phase spreading on a solid surface <cit.>, the concentration, volume fraction or density of a phase in a mixture <cit.>, etc.
  
  By introducing an auxiliary unknown  w, the fourth-order equation (<ref>) can be written in the following system of second-order equations
  
      ∂_tu   =∇·(f(u)∇ w)   in Ω×(0,T),
      
      w   =-γΔ u +φ'(u)   in Ω×(0,T).
  
  System (<ref>) is to be supplemented by the initial condition
  
      u(x,0)=u_0(x)   in Ω,
  
  and the no-flux boundary conditions
  
      ∇ u·n=∇ w·n=0   on ∂Ω×(0,T),
  
  where n is the exterior unit normal vector to ∂Ω.
  
  There are two important properties for the  equations (<ref>) with the natural boundary conditions (<ref>). The first one is the conservation of mass. In fact,
  
      ddt∫_Ωu dx=∫_Ω∂_tu dx=∫_Ω∇·(f(u)∇ w) dx=∫_∂Ωf(u)∇ w·n dx=0,
  
  that is the total amount of the phase in the domain must always be equal to the given initial amount of this phase. The equation is associated with the total free energy functional
  
      J(u)=∫_Ω(γ/2|∇ u|^2+φ(u)) dx,
  
  which satisfies the total free energy dissipation, that is
  
      d Jdt=∫_Ω(γ∇ u ∇∂_tu+φ'(u)∂_tu) dx=∫_Ωw ∂_tu dx=-∫_Ωf(u)|∇ w|^2 dx⩽ 0.
  
  
  A weak formulation of problem (<ref>) can be obtained by multiplying (<ref>) by test function v. Integrating by parts and using the boundary conditions (<ref>) eventually lead to: Find u ∈ H^2(Ω) such that
  
      ∫_Ω∂_tu v dx+∫_Ωγ f(u)Δ uΔ v dx+∫_ΩγΔ u∇ f(u)·∇ v dx+∫_Ωf(u)∇φ'(u)·∇ v dx=0,
  
  for all v ∈ H^2(Ω). A weak formulation for problem (<ref>) can be obtained by multiplying the first and second equation in (<ref>) by test functions v and q, respectively: Find (u ,w) ∈ H^1(Ω)× H^1(Ω) such that
  
      ∫_Ωu_t v  dx+∫_Ω f(u)∇ w·∇ v  dx=0,  ∀ v ∈ H^1(Ω),
         ∫_Ω w q dx-∫_Ωγ∇ u·∇ q dx-∫_Ωφ'(u) q dx=0,  ∀ q ∈  H^1(Ω).
  
  Theoretical results on the existence of solutions to (<ref>) and (<ref>) in the sense (<ref>) and (<ref>), respectively, can be found in <cit.>. 
  
  
  
  
  § NUMERICAL METHODS
  
  
  
  
   §.§ Spatial discretization
  
  Kästner et al. <cit.> showed that the convergence rates for approximations that are less than second-order in space are suboptimal when using the direct weak formulation (<ref>). This is due to the discretization of higher-order differential operators.  The weak formulation (<ref>) naturally leads to mixed finite element methods which provide robust numerical approaches to compute solutions of high-order PDEs by avoiding the discretization of high-order derivatives <cit.>.  In the following, we will use the mixed weak formulation (<ref>).
  
  We discretize in space by continuous piecewise Lagrange finite elements. Given a polygonal domain Ω, we consider  𝒯_h being a uniform partitioning of the domain Ω into disjoint triangles κ, so that 
  
      Ω=⋃_κ∈𝒯_hκ.
  
  The  element mesh size h is defined as 
  
      h:=max_κ∈𝒯_hh_κ,
  
  where h_κ:=diam(κ). Associated with the partitioning 𝒯_h  is the finite element space
  
      𝒱_h={ v_h∈𝒞^0(Ω,ℝ): v_h_|_κ∈ℙ_k,  ∀κ∈𝒯_h }⊂ H^1(Ω),
  
  where ℙ_k denotes the space of polynomials of degree less than or equal to k on any element κ. 
  
  The resulting semi-discrete (in space) weak formulation of problem (<ref>)-(<ref>) is given as follows: Find ( u_h,w_h)∈𝒞^1(0,T;𝒱_h×𝒱_h) such that, for all t∈[0,T],
  
      ∫_Ω∂_tu_h(t)v_h  dx+∫_Ω f(u_h(t))∇ w_h(t)·∇ v_h dx=0,  ∀ v_h  ∈𝒱_h ,
         ∫_Ωw_h(t) q_h dx-∫_Ωγ∇ u_h(t)·∇ q_h dx-∫_Ωφ'(u_h(t)) q_h dx=0,  ∀ q_h ∈𝒱_h,
         u_h(0)=Π_hu_0,
  
  where Π_h is an interpolation or projection operator on 𝒱_h.
  
  
  
   §.§ Temporal discretization
  
  The time interval [0,T] is discretized as
  
      t_n=n Δ t,    n=0,1,...,N,    Δ t=TN,
   
  where Δ t is the time step used. We consider a uniform discretization in time by a semi-implicit stepping method, so that u_h^n≃ u_h(t_n) for n=0,1,..., N. For the temporal discretization of the equations (<ref>), all linear terms  are approximated implicitly. The time derivative ∂_tu_h is approximated using a second-order backward differentiation formula (BDF2). To allow for a linear scheme, a second-order extrapolation formula is employed for the mobility function f(u_h) while the energy term φ'(u_h) is evaluated using the Taylor approximation of second-order truncation error. We have noticed from numerical experiments that the approximation of the nonlinear energy term using the proposed technique is better than by an extrapolation which leads to severe stability constraint. Our technique introduces a linear implicit term which improves the stability constraint. We emphasize two main benefits resulting from the proposed scheme. First of all, it is easy to implement. No fixed point iteration for solving the system is required. Only one linear system needs to be solved at every time step. The proposed scheme is written as follows: given a suitable approximation of the initial solution u_h^-1=Π_hu_0∈𝒱_h and a proper initialization for u_h^0∈𝒱_h, find (u_h^n,w_h^n)∈𝒱_h ×𝒱_h such that
  
      ∫_Ω(3u_h^n-4u_h^n-1+u_h^n-22Δ t)v_h dx+∫_Ω(2f(u_h^n-1)-f(u_h^n-2))∇ w_h^n·∇ v_h dx=0,  ∀ v_h ∈𝒱_h,
         ∫_Ω w_h^nq_h  dx-∫_Ωγ∇ u_h^n·∇ q_h  dx-∫_Ω(φ'(u_h^n-1)+φ”(u_h^n-1)(u_h^n-u_h^n-1))q_h dx=0,   ∀ q_h ∈𝒱_h,
  
  for all 1⩽ n⩽ N. 
  
  The numerical algorithm (<ref>) is a two steps scheme and therefore, it requires the use of a starting procedure to obtain u_h^0. Here, we use a semi-implicit backward Euler method 
  
      ∫_Ω(u_h^0-u_h^-1Δ t)v_h dx+∫_Ωf(u_h^-1)∇ w_h^0·∇ v_h dx=0,  ∀ v_h ∈𝒱_h,
         ∫_Ω w_h^0q_h  dx-∫_Ωγ∇ u_h^0·∇ q_h  dx-∫_Ωφ'(u_h^-1)q_h dx=0,   ∀ q_h ∈𝒱_h,
  
  to compute u_h^0 starting from  u_h^-1.
  
  The numerical treatment of the energy term φ'(u_h) uses the following semi-implicit Taylor expansion of φ'(u_h(t_n)) about u_h=u_h(t_n-1): 
  
      φ'(u_h(t_n-1))+φ”(u_h(t_n-1))(u_h(t_n)-u_h(t_n-1))=φ'(u_h(t_n))-1/2(u(t_n)-u(t_n-1))^2φ^(3)(u^n),
  
  where u^n is between u_h(t_n-1) and u_h(t_n). To evaluate the local truncation error
  
      τ^n(q_h)=∫_Ω-1/2(u(t_n)-u(t_n-1))^2φ^(3)(u^n)q_h dx
  
  of the energy term, one expands again u_h(t_n-1) in (<ref>) around the point t=t_n and gets
  
      τ^n(q_h)=-Δ t^2/2∫_Ω(∂_tu_h(ξ^n))^2φ^(3)(u^n)q_h dx,
  
  where ξ^n=ξ^n(x)∈[t_n-1,t_n] for all x∈Ω. The backward differentiation formula used for the numerical treatment of ∂_t u_h is known to have error that is second-order in time and one can easily show that the extrapolation technique used for the approximation of the mobility term f(u_h) has an order of accuracy of 𝒪(Δ t^2). Hence, the scheme (<ref>) is second-order accurate in time. 
  
  
  
  
   §.§ Solution of the discrete system
  
  Equations in (<ref>) or (<ref>) consist in a coupled linear system in u_h^n and w_h^n of the following form
  
      [  M_1  A_1; -A_2  M_2 ][ U; W ]
      =
      [ F; G ],
  
  where M_i and A_i, i=1,2, are mass and stiffness matrices, respectively; U and W correspond to the finite element degrees of freedom for u_h^n and w_h^n, respectively; F and G are properly chosen constant vectors.
  
  This linear system shares characteristics with those for saddle point problems, such as the block-structured sparsity pattern, the lack of positive definiteness, etc. These linear systems are difficult to solve. The use of a direct solver may lead to a lot of fill-in, while iterative solvers may be impacted by the indefiniteness of the matrix. For our 2D test cases, we used the sparse linear solver UMFPACK <cit.>, which mitigates the fill-in and provides a relatively efficient solution. However, for 3D computations, one would need to use iterative methods. One may attempt to build a dedicated iterative method accounting for the fact that all the blocks M_i and A_i are at least semi-positive definite, a property usually not shared by general saddle point problems.
  
  
  
  
   §.§ Mass conservation
  
  As usual in finite element methods with no-flux boundary conditions, the scheme (<ref>) is conservative. In fact, taking v_h=1 in the first equation of (<ref>) and (<ref>), respectively, leads to
  
      ∫_Ω3u_h^n-4u_h^n-1+u_h^n-22Δ t dx=0     and    ∫_Ωu_h^0-u_h^-1Δ t dx=0.
  
  It follows by induction from (<ref>) that 
  
      ∫_Ωu_h^n dx=∫_Ωu_h^0 dx,  ∀ n=1,...,N,
  
  which means that the discrete total mass is conserved over time.
  
  
  
  
  § NUMERICAL RESULTS
  
  
  A series of numerical problems in two dimensions is given in this section to test the performance of the numerical method developed in section <ref>. All the algorithms for the simulations are implemented using FreeFem++ software <cit.>.
  
  
  
   §.§ Numerical results — Convergence
  
  In this subsection we numerically show the rates of convergence in space and time of the algorithm (<ref>)-(<ref>). We first conduct a detailed investigation using the method of manufactured solution for a constant mobility function. We then use an analytic solution for a non-constant mobility function to perform a similar error analysis.
  
  
  
    §.§.§ Manufactured solution
  
  
  Following the framework presented in <cit.>, we employ a manufactured solution to test our method. The idea is to use an arbitrary regular function as analytic solution for the computation of numerical errors.
  
  
  Substituting in (<ref>) the  function
  
      u=(t+1)sin(απ x),
  
  where α∈ℝ, by assuming a free energy function
  
      φ(u)=u^4/4-u^2/2,
  
  and a constant mobility function f(u)=1, yields 
  
      {   u_t=∇·(∇w)+S,
         w=-γΔu+φ'(u),
      .
  
  where the residual S is given by
  
      S=   sin(απ x) + 3(t+1)^3α^2π^2sin^3(απ x) + γ(t+1)α^4π^4sin(απ x) - (t+1)α^2π^2sin(απ x)
         -6(t+1)^3α^2π^2sin(απ x)cos^2(απ x)
  
  and the auxiliary variable w is given by
  
      w=γ(t+1)α^2π^2sin(απ x) + (t+1)^3sin^3(απ x) - (t+1)sin(απ x).
  
  We assumed here a constant mobility function f(u)≡ 1 to simplify the calculation of the residual function S, but the energy term (<ref>) is a nonlinear non-zero function. This will allow to test the impact of the new proposed technique for the discretization of the nonlinear energy term, based on a semi-implicit Taylor expansion of the latter.
  
  For the convergence analysis in space, we compute the L^2- and H^1-error norms  on the discrete solution u_h of (<ref>) at the final time, which are defined as 
  
      ‖u-u_h‖_L^2(Ω)=(∫_Ω|u-u_h|^2 dx)^1/2,
         ‖u-u_h‖_H^1(Ω)=(∫_Ω(|u-u_h|^2+|∇(u-u_h)|^2) dx)^1/2.
  
  A set of boundary conditions for the system has to be imposed. We will investigate both Dirichlet and non-homogeneous Neumann boundary conditions using the exact solution (<ref>) and (<ref>) at boundaries. The linear (ℙ_1) and quadratic (ℙ_2) Lagrange finite element discretizations and four uniform meshes (25× 25, 50× 50, 100× 100 and 200× 200 elements) of the unit square (centered at the origin) domain  have been considered. We use a fixed time step Δ t= 10^-5 and consider a small value for the parameter γ=10^-4. The final computational time is T=10^-3.
  
  We first impose Dirichlet conditions u_h=Π_hu and w_h=Π_hw on boundaries.  The numerical errors on u_h are computed for different values of the parameter α. Figure <ref> shows theses errors as a function of the element size h. Optimal convergence rates, i.e. 𝒪(h^k+1) and 𝒪(h^k),  are obtained for the L^2- and H^1-error, respectively. Better, we get super convergence for the L^2-norm error with the quadratic Lagrange elements. The error increases with α, which controls the gradient of the solution. However, the rates of convergence remain the same, and therefore, are independent from the choice of α. 
  
  
  We next impose non-homogeneous Neumann conditions ∂u_h/∂𝐧=Π_h∂u/∂𝐧 and ∂w_h/∂𝐧=Π_h∂w/∂𝐧. Here, also optimal convergence rates in the L^2- and H^1-error are obtained. Again, super-convergence in the L^2-norm is observed with quadratic Lagrange elements. To facilitate comparison between the results obtained using Dirichlet and Neumann conditions, the errors computed for α=2 are shown in Figure <ref>.  We can see that the L^2-error with Neumann conditions are larger than with Dirichlet conditions, but the H^1-error are almost identical.
  
  
  We have performed a similar error analysis with the L^2- and H^1-error on the auxiliary variable w_h and found the same orders of convergence as for u_h. However, the H^1-error on w_h is larger and the order of magnitude of the difference increases with α. This is natural since w_h=-γΔu_h+φ'(u_h) is more difficult to approximate and the gradient of the solution increases with α more rapidly than for u_h.
  
  To test the temporal convergence rate, we compute a reference solution (u_h^ref,w_h^ref) on a grid of 50× 50 elements using a very small time step Δ t= 10^-8 at the final time T=10^-3. We choose the quadratic Lagrange element and impose Dirichlet boundary conditions to minimize the spatial error. Moreover, with quadratic elements, fixed computational time and time step give smaller error than with ℙ_1 elements (see Table <ref> in section <ref>).  We analyze the numerical convergence in time by varying the time step Δ t systematically between 0.125× 10^-4 and 2× 10^-4. For each time step, the numerical errors in the L^2- and L^∞-norms  are computed. Figure <ref> shows the L^2- and L^∞-error on u_h and w_h, respectively, and the order of convergence of the numerical solution with respect to the reference solution.  We obtain almost a second-order convergence rate in time.
  
  
  
  
    §.§.§ Analytic solution
  
  
  In the previous subsection, optimal convergence rates have been obtained for constant mobility function. It is the main purpose of this subsection to numerically study the convergence rates for non-constant mobility function by using a particular analytic solution.
  
  Setting f(u)=u, φ=0 and γ=1 in (<ref>), we get the nonlinear fourth-order equation
  
      u_t+∇·(u∇Δ u)=0   in Ω×(0,T).
   
  Equation (<ref>) has a nonnegative closed-form compactly-supported, 2-dimensional self-similar solution <cit.> given by
  
      u(x,y,t)=
      {   t^-1/3192(L^2-r^2)^2,   r<L,
         0,   r⩾ L,
      .
             r=√(x^2+y^2)t^1/6.
  
  In the framework of the mixed formulation, we need the auxiliary variable w=-Δ u, which is calculated in a straightforward way from (<ref>) to get
  
      w(x,y,t)=
      {   t^-2/324(L^2-2r^2),   r<L,
         0,   r⩾ L.
      .
  
  The auxiliary variable (<ref>) is not continuous at point where r=L. Moreover, the source-type solution (<ref>) tends to a Dirac mass as t goes to zero. Therefore, we have to choose L large enough and  a starting time t_0>0 to ensure the regularity of the solution (<ref>)-(<ref>) for the mixed formulation. Here, there is no energy term, but the mobility function depends on the solution. This will allow to test the impact of the numerical treatment for the  mobility function, based on extrapolation. 
  
  To test the spatial convergence rates, we compute the L^2- and H^1-error norms between the exact solution (<ref>)–(<ref>) and the numerical solution computed with the proposed method (<ref>)-(<ref>), where we set φ≡0. We use the same computational domain, meshes, finite element approximations as in section <ref>. Since, we focus on the spatial convergence rates, we use a small time step Δ t=10^-8 to minimize the temporal error from the extrapolation of the non constant mobility. We choose L=3 and a starting time t_0=0.001. Tables <ref> and <ref> show  the  computational time, the L^2- and H^1-error as function of the element h, computed with Dirichlet and Neumann boundary conditions, respectively. The formula to obtain the order of convergence in space, p, reads:
  
      p=ln(‖ u-u_h‖‖ u-u_h/τ‖) / lnτ,
  
  where τ>0 is the refinement factor between consecutive element sizes. The value τ=2 was chosen in our numerical experiments. By using the formula (<ref>), we  get from Table <ref> that the convergence rates in the H^1-norm are about 1 and 2 for the linear and quadratic Lagrange elements, respectively, while the convergence rates in the L^2-norm are approximately 2 and 4 for the linear and quadratic Lagrange elements, respectively. These same convergence rates were obtained with a manufactured solution in section <ref> (see Figure <ref>). Formula (<ref>), used on Table <ref>, gives approximately the same rates of convergence as in Figure <ref> (for Neumann boundary conditions). While the quadratic Lagrange elements provide smaller error than the linear Lagrange elements, they are also the most computationally expensive. However, the quadratic Lagrange elements with mesh size 2h provide superior accuracy than the linear elements with mesh size h, and both require approximately the same computational cost. Therefore, high-order finite element approximations are beneficial in terms of efficiency. We also have computed the L^2- and H^1-error on the auxiliary w_h and have found the same orders of convergence as for u_h, except in the H^1-error for quadratic Lagrange elements, which leads to super-convergence of order 𝒪(h^3).
  
  
  
  
  For the error and convergence analysis in time, we compute a reference solution (u_h^ref,w_h^ref) on  a grid of 100× 100 using linear Lagrange elements and a time step Δ t= 10^-5. We then analyze the numerical convergence in time by varying the time step Δ t systematically between 0.001 and 0.008. For each time step, the errors on u_h and w_h, computed as 
  
      ‖ u_h^ref -u_h ‖_L^2(Ω)‖ u_h^ref‖_L^2(Ω)   and   ‖ w_h^ref -w_h ‖_L^2(Ω)‖ u_h^ref‖_L^2(Ω),
  
  respectively, are recorded at the final time T=0.201. Table <ref> shows the error, computational time and order of convergence of the numerical solution with respect to the reference solution. The rate of convergence in time is calculated using the formula (<ref>), where τ is the refinement factor between consecutive time steps. Here, also τ=2 was chosen in our numerical tests. A second-order convergence rate  in time is obtained.
  
  
  
  
  
  
   §.§ Numerical results — Applications
  
  
  
    §.§.§ Lubrication
  
  Equation (<ref>) appears in the lubrication approximation of a thin viscous film, which is driven by surface tension. In this context, u describes the height of the liquid film which spreads on a solid surface <cit.>. The numerical solutions of high-order nonlinear degenerate diffusion equations like  (<ref>) require special numerical techniques in order to preserve the non-negativity of solutions.  Constructing finite element methods that guarantee the non-negativity of the solutions of (<ref>), has been initiated by Barrett et al.  <cit.> by enforcing at each time step the non-negativity of the approximate solution as a constraint, which leads to a variational inequality. The solutions of the one-dimensional version of (<ref>) using a non-negativity preserving scheme was carried out in <cit.>. The solution of a modified problem, which consists in approximating the mobility term f(u)=u by
  
      f_ξ(u)=u^5ξ u + u^4,
  
  has been shown by <cit.> to remain positive. This regularization is used in <cit.> to compute numerical solutions of (<ref>) for the initial data
  
      u(x,y,0)=δ+Ce^-σ(x^2+y^2),
  
  where δ>0 represents the thickness of an ultra-thin precursor film under the Gaussian fluid droplet centered at the origin.   
  Here, we solve (<ref>) without regularization, i.e. with the degenerate mobility function f(u)=u. The computational domain is the rectangle [-0.5,0.5]×[-1,1]. Homogeneous Neumann conditions (<ref>) are imposed on the boundary.  We use linear Lagrange elements on a uniform mesh 70× 140, i.e. a mesh size h=0.02. We set δ=0.01, σ=80 and C=2. The same values for these parameters and mesh size were used in <cit.>. The contours of the numerical solutions and a 3D view of the droplet  at different times are shown in Figure <ref>. 
  
  The algorithm (<ref>)-(<ref>) preserves the  positivity of the solution up to a maximal time step of Δ t≃ 10^-5. We have also performed numerical simulations for C=1 as in <cit.>, for which the positivity has been preserved up to a maximal time step of Δ t≃ 0.01, while Δ t=10^-5  was found  in <cit.> as the maximal time step to preserve the positivity.
  
  To show that the numerical scheme satisfies the energy decreasing property (<ref>), we define the discrete energy functional of (<ref>) as
  
      ε_h(u_h^n)=1/2∫_Ω|∇ u_h^n|^2 dx,
  
  and consider the evolution of the discrete energy. Figure <ref>(a) shows the time-evolution of the discrete total energy. Our results demonstrate that the total discrete energy is non-increasing which agrees well with the total energy dissipation property. 
  
  The maximum of the solution of (<ref>) with the initial condition (<ref>) occurs at the origin for all times. Figure <ref>(b) shows the height of the droplet at the origin as time evolves. After an initial transient phase (t>10^-5), the maximal height u(0,t) evolves in t^-1/3 for some times, and in t^-1/5 for longer times, before approaching a uniform flat state for large times (t>10^-1). The same evolutions are observed in <cit.>. 
  
  
  
  
    §.§.§ Phase separation
  
  Spinodal decomposition is a mechanism by which a solution of two or more components can separate into distinct phases. Spinodal decomposition in a binary mixture is well-described by the Cahn-Hilliard equation <cit.>. Here,  we are interested in the phase separation of a homogeneous binary mixture governed by (<ref>) to further test our method (<ref>)-(<ref>). As in <cit.>, we choose a mobility function f(u)=(1-u^2), a free energy function φ(u)=1/4ϵ^2(u^2-1)^2 and set the parameters γ=1 and ϵ=0.03.  In this case, equation (<ref>) describes the state of a relative concentration u=(2u_a/u_a+u_b-1), where u_a and u_b are the concentrations of the two components <cit.>. With this double well form of the energy function, the system is locally in one of the two phases if u is close to one of the two minima ± 1 of φ <cit.>.  
  
  The initial condition is chosen as
  
      u_0(x)=u+ ∑_i=1^M c_ir_i(x,a_i),
  
  where u∈[-1,1], a_i∈ℝ^2, c_i∈ℝ, for i=1,...,M, and r_i(x,a_i)=e^-1000|x-a_i|^2. We have randomly generated M=1000 points a_i in the square domain [0,1]×[0,1], and 1000 scalars c_i∈[-0.01,0.01]. All numerical results in this section are computed on a 64× 64 grid with quadratic Lagrange elements. Figure <ref> shows the initial condition for two different values of u. 
  
  
  We first consider the initial solution for u=0. The dynamics of the binary fluid are shown in Figure <ref>. From its initial concentration, the system separates quickly into two phases whose composition is determined by the two minima ± 1 of the free energy function. Small regions of the same phase combine together to form larger patterns, reducing the interface boundaries between the two phases.
  
  
  A particular regime of interest is when one phase is significantly more abundant in the initial mixture. In such case, the minority phase emerges and forms droplets. Larger particles grow at the expense of the smaller ones which
  shrink and disappear. This form of competitive growth is known as the Ostwald ripening <cit.>. To illustrate such regime, we use the initial solution for u=0.4. Figure <ref> shows the time-evolution of the binary mixture computed. We can observe the formation of droplets. Smaller droplets are absorbed into larger ones through diffusion, exhibiting clearly the phenomenon of Ostwald ripening. 
  
  
  We qualitatively observe the same behavior of the solution as in <cit.>. The plots of the energy evolution in Figure <ref> illustrates that the total discrete energy satisfies the energy dissipation property. Plots of the discrete total mass are shown in Figure <ref>, which  illustrate the mass-conservative property of the linear scheme.
   We observe that the magnitude and the nonlinearity of the mobility have a significant direct impact on the stability of the scheme. We have noticed from numerical experiments that a small time step is required for the stability of numerical solutions involving non constant mobility functions. The time step could be taken larger for constant mobility functions and regular energy functions. The parameter γ has an impact on the length scale of the diffuse interface by balancing the surface energy term with the bulk free energy <cit.>. In terms of the stability of the scheme, there is a strong link between the parameter γ and the free energy term. 
  
  
  
  
  
  
    §.§.§ Electrowetting
  
  Electrowetting has become one of the most widely used tools for manipulating tiny amounts of liquids on surfaces. It consists in the modification of the wetting properties of a surface with an applied electric field. A diffuse-interface model for drop motion, due to electrowetting, in a Hele-Shaw geometry is proposed in <cit.>. The model is governed by the Cahn-Hilliard model with a spatially dependent energy
  
      {   ϵ u_t-∇·(u∇ w)=0,
         w=-ϵ^2Δ u +φ'(u) -ϵρ(x),
      .
  
  with a double well energy function φ(u)=(u-ϵ)^2(u-1)^2, where ρ(x)=λχ(x) defines a local energy (electric field) with characteristic function χ(x) and strength λ. The local energy  term introduces the electrowetting into the diffuse-interface model. The parameter  ϵ=0.0427 controls the diffuse interface thickness. For more details and the experimental setup, we refer to <cit.>.
  
  We solve this problem using a Gaussian initial data
  
      u_0(x)=δ+e^-10|x-a|^2,
  
  representing a droplet centered at the point a that will be specified for each test case. All numerical results in this section are computed on a 47× 94 grid in the rectangular domain [-0.5,0.5]×[-1,1] using linear Lagrange elements. This gives rise to a mesh size h≈0.03 as in <cit.>. We use a precursor film δ=ϵ  and a fixed time step Δ t=0.001 instead of adaptative time steps as in <cit.>.
  
  We first consider a droplet centered at point a=(0,-0.3) and then apply the electric field ρ(x) on the surface y>0. The time-evolution of the droplet is displayed in Figure <ref>. The surface of the droplet is deformed during its translation toward the region with low energy (with an electric field). The variation of the curvature of the droplet increases with the strength λ of the electric field. The same evolutions are observed using numerical simulations in <cit.> and experimentally in <cit.>.
  
  
  We next consider a droplet centered at the origin a=(0,0) and apply two electric fields of same strength λ=2 on the surfaces y<-0.3 and y>0.3, respectively. For this test case, we also want to compare the time it takes for a nonlinear scheme, namely the second-order backward differentiation formula (BDF2), and for the linear scheme introduced in our study. The BDF2 scheme is second-order accurate in time and consists in approximating the temporal derivative as in (<ref>) and treating all other terms implicitly. This fully implicit method uses an iterative procedure to solve the nonlinear resulting system. Here, we use Newton method with a tolerance of 10^-6 for the iterative solver. Numerical results for both methods are displayed in Figure <ref>, which shows the  capability of the proposed method to simulate the macroscopic dynamics of drop splitting as with the nonlinear BDF2 scheme but with a relatively smaller computational time. In fact, we can see from Table <ref> that the algorithm (<ref>)-(<ref>) is in average four times faster than the BDF2 scheme. Though the maximal time step for the proposed scheme can be smaller than for BDF methods, the required CPU time is still significantly smaller than with BDF methods.
  
  
  
  
  
  § CONCLUSION
  
  
  The present study provides a new fully discrete mixed finite element method to approximate fourth-order nonlinear diffusion equations. We develop a special numerical technique of approximating the nonlinear energy term that allows to design a linear scheme for our system. The performance of our scheme is tested using relevant numerical examples, and the results clearly demonstrated the ability of the proposed scheme to resolve fourth-order nonlinear diffusion equations. Manufactured and analytic solutions have been used for the numerical investigation of the accuracy of the scheme. We have performed numerical analysis and presented convergence results showing the good behavior of the numerical method with respect to mesh refinement, order of the  approximation and boundary conditions. With respect to the spatial discretization, convergence rates are obtained that match analytical error estimates available for linear problems. We also have presented a series of numerical experiments to confirm the second-order accuracy in time for both constant and non-constant mobility function. Some complex and real-world examples are presented including a model for electrowetting on dielectric and for coarsening dynamics of binary mixture.
  
  
  
  
  
  § ACKNOWLEDGEMENT
  
  The first and second authors are supported by the UM6P/OCP group of Morocco. The third author is supported through a Discovery Grant of the Natural Sciences and Engineering Research Council of Canada.
  
  
  
  § REFERENCES
  
  plain 
  
  