Finite Element Method for a Nonlinear Differential Equation Describing Crystal Surface Growth

In this paper, for a nonlinear differential equation describing crystal surface growth, the finite element method is presented. A nice order error estimates is derived by means of a finite element projection approximation.


Introduction
Molecular beam epitaxy (MBE) is a widely practiced technique for depositing atoms from a vapor phase onto a surface. The atoms stick to the surface and then diffuse until they find a favorable site. This technique is very important in growing thin films (see [8]). Recently, there have been some attempts to link MBE with surface growth ideas from statistical physics and several experimental studies exhibiting MBE growth (see [9,14]).
Under conditions typical of molecular beam epitaxy (MBE), evaporation and the formation of bulk defects can be neglected. The height H(x, t) of the surface above the substrate plane then satisfies a continuity equation, where F is the incident mass flux out of the molecular beam. Since we are interested in large scale features we neglect fluctuations in F ("shot noise") and in the surface current ("diffusion noise"). As a general rule, J surface , which is the systematic current, depends on the whole surface configuration. If we keep only the most important terms in a gradient expansion, subtracting the mean height H = F t, use appropriately rescaled units of height, distance and time [13], Equation (1.1) attains the followsing form where γ and µ are two positive constants. Here, we use the common practice and disregard contributions to the current that are even in h, such as ∇(∇h) 2 , though they may well be relevant for the coarsening behavior of the surface [12,13]. In Equation (1.2), the linear term describes relaxation through adatom diffusion driven by the surface free energy, the second nonlinear term models the nonequilibrium current [13,14]. Assuming in-plane symmetry, it follows that the nonequilibrium current is (anti)parallel to the local tilt ∇h, with a magnitude f (∇h 2 ) depending only on the magnitude of the tilt. For the function f (∇h 2 ), within a Burton-Cabrera-Frank-type theory (see [10,12]), for small tilts the current is proportional to |∇h|, and in the opposite limit it is proportional to |∇h| −1 . This suggests the interpolation formula f (s 2 ) = 1 1+s 2 [13]. Using u(x, t) to replace h(x, t), considering the 1D case, we have the following form ∂u ∂t where Ω = (0, π). On the basis of physical consideration, the equation is supplemented by the following boundary conditions and the initial condition Remark 1. There are four boundary conditions (BCs) for the equation (1.3), e.g. two BCs at x = 0, π. Noticing that the essential BCs are u x (0, t) = u x (π, t) = 0. So, we define and the weak formulation of problem (1.3)-(1.5) in the following.
During the past years, many authors have paid much attention to Equation (1.3). In [11], in the limit of weak desorption, O. Pierre-Louis et al. derived Equation (1.3) for a vicinal surface growing in the step flow mode. This limit turned out to be singular, and nonlinearities of arbitrary order need to be taken into account. Recently, M. Grasselli, G. Mola and A. Yagi [7] showed that the 2D case of equation (1.3) endowed with no-flux boundary conditions generates a dissipative dynamical system under very general assumptions on ∂Ω on a phase-space of L 2 -type. In [16], based on the iteration technique for regularity estimates and the classical existence theorem of global attractors, Zhao and Liu proved that the 2D case of problem (1.3)-(1.5) possesses a global attractor on some affine space of H k (0 ≤ k < +∞).
The finite element method is an essentially discretization method for the approximate solution of partial differential equations. It has the natural advantage in keeping the physical properties of primitive problems. There are many papers which have already been published to study the finite element method for fourth order nonlinear parabolic equation. In [15], using finite element method, Zhang consider the semi-discrete approximation and Full-discrete approximation for 1D Cahn-Hilliard equation which is also a fourth-order nonlinear parabolic equation. In Elliott and French's paper [6], for 2D Cahn-Hilliard equation, a continuous in-time finite-element Galerkin approximation was considered. They use the nonconforming Morley element and derive optimal order error bounds in L 2 . Barrett, Blowey and Garcke [1] have consider a fully practical finite element approximation of the fourth order nonlinear degenerate parabolic equation where generically b(u) := |u| p for any given p ∈ (0, +∞). An iterative scheme for solving the resulting nonlinear discrete system was analysed. In addition to showing well-posedness of their approximation, they proved convergence in one space dimension. For more recent results we refer the reader to [2,4,5] and the references therein.
In this paper, we consider the finite element analysis for problem (1.3)-(1.5). Noticing that the existence of a solution locally in time is proved by the standard Picard iteration, global existence results are obtained by proving a priori estimate for the appropriate norms of u(x, t). Adjusted to our needs, the results is given in the following theorem. The outline of this paper is as follows. In the next section, we establish a semi-discrete approximation and derive its error bound. In Section 3, the fulldiscrete approximation for problem (1.3)-(1.5) is studied. In the last section, some numerical experiments which confirm our results are performed.
Throughout this paper, we denote L 2 , L p , L ∞ , H k norm in (0, 1) simply by · , · L p , | · | ∞ and · k . Define the inner product of L 2 space as (·,·), the space On the other hand, the letter C denotes generic constant independent of the finite element division size not necessarily the same at different occurrences.

Semi-Discrete Approximation
Let I h : 0 = x 0 < x 1 < · · · < x N = π be a finite element division for the interval be the piecewise polynomial spline space with the degree k ≥ 3, and The weak formulation of problem (1.3)-(1.5) reads: It is easy to see that the conservation of mass for (2.2) holds as for the classical solution. Setting where C is a positive constant depends only on γ, µ and T , independent of h.
Proof. The equation of problem (2.2) is an ordinary differential equation and according to ODE theory, there exists a unique local solution to the problem (2.2) in the interval [0, t n ). If the estimate (2.3) is obtained, then according to the extension theorem, we can obtain the existence of unique global solution. Setting where β = µ 2 γ . Integrating (2.6) with respect to the time t, we get Differencing F h (t) with respect to t, using (2.8), we get Noticing that Adding the above two inequalities together gives We also have It then follows from (2.8) that Therefore, we get (2.12) By (2.7), (2.10), (2.11) and (2.12), we complete the proof of Theorem 2.
In order to consider error estimate, we first introduce a finite element approximation projection for a steady-state problem. Let u ∈ H 2 (2.13) It then follows (2.13) that where c 0 is a positive constant depends only on γ and µ. Hence, a(u, v) is a symmetrical positive determined bilinear form, and there exist a unique solution h for problem (2.13). Based on the standard finite element method (see [3]), we have Now, we consider the error estimate for the semi-discrete finite element solution. Let u be the solution of (2.1), u h be the solution of (2.2). Denote (2.16) Note that u, u h and R h u satisfy (2.1), (2.2) and (2.13). Then, θ(t) satisfies Lemma 1. Let u be the solution of (2.1), u h be the solution of (2.2), u ∈ L ∞ (0, T ; H 2 (I)). Then, ∀ε ∈ (0, 1), there exists a constant C = C(u, u h (0), ε) such that Proof. First of all, we give some estimates which will be used in this proof. It follows from Theorem 2, (2.13), (2.14) that Set g(s) = 1 1+s 2 . Then, g (s) = 1 1+s 2 − 2ss (1+s 2 ) 2 . By Theorems 1 and 2, we have D 2 u ≤ C, D 2 u h ≤ C. Therefore, there exists a positive constant ∈ (0, 1) such that Hence, the proof of Lemma 1 is completed.
Setting v h = θ in (2.17), using Cauchy's inequality, we immediately conclude that It then follows from the above inequality that Using Lemma 1, setting ε = γ 2 2 , we have By Gronwall's inequality, we deduce that Combing (2.20) and (2.15) (Noticing that (R h u) t = R h u t ), using triangle inequality, we complete the proof of Theorem 3.

Full-Discrete Approximation
For any given positive integer M , let ∆t = T /M denote the size of time discretization. Denote U n = U (x, t n ) for t n = n∆t, n = 0, 1, . . . , M . We introduce the backward Euler difference formula Here, if g(t) is a continuous function of t, let g n = g(t n ). Now, we define the full-discrete finite element form to approximate problem (2.1): Find U n ∈ S (k) h (n = 1, 2, . . . , M ) such that For the above form, if U n−1 is known and ∆t is sufficiently small, by solving a positive definite system of linear equations which is equal to equation (3.2), we can obtain U n . Let Using (3.2) and (2.13), we derive that Theorem 4. Let u be the solution of (2.1), U n be the solution of (3.2), u(0) ∈ H k+1 (I), u t ∈ L 2 (0, T ; H k+1 (I)), u tt ∈ L 2 (0, T ; L 2 (I)), ∆t/h 2 ≤ c, and U 0 ∈ S where i = 0, 1. Then if h is sufficiently small, there exists a constant C = C(u) which is independent of h, ∆t and n, such that Proof. First of all, we give a posterior hypothesis: There exits a h 0 , when 0 < h ≤ h 0 , we have We will prove the correctness of (3.6) in the end of this proof.
Setting v h = θ n in (3.4), we derive that (3.7) Using ( In addition, we have We also have 2C Dθ n 2 = −2C θ n , D 2 θ n ≤ γ 2 D 2 θ n 2 + 2C 2 γ θ n 2 . Taking above estimates into (3.7), we derive that Taking the sum of n, noticing that Dθ 0 ≤ Ch k u(0) 1 , n∆t = t n ≤ T , we obtain Let ∆t be sufficiently small, which satisfies 2C 2 γ ∆t ≤ 1 2 , we deduce that Using the discrete Gronwall inequality, (2.13) and triangle inequality, we obtain (3.5). Now, in order to complete the proof of Theorem 4, we only need to prove the posterior hypothesis (3.6). We can use inductive method. when m = 0, based on the initial approximation assumption and the finite element inverse inequality, letting h ≤ h 0 and h 0 be sufficiently small, we obtain (3.6). If we assume that (3.6) is correct for m = l − 1, based on the above proof, we can easily obtain that the estimate (3.5) is correct for n = l, where C is a constant independent of n, ∆t and h (noticing that t n ≤ T ). Using finite element inverse inequality, interpolation approximation properties and (3.5), we have where u i ∈ S (3) h is the Hermite type interpolation approximation of the function u. Hence, (3.6) is correct for m = l. Then, using the inductive method, the rationality of (3.6) is proved. Therefore, we complete the proof of Theorem 4.
We get the solution which evolves from t = 0 to t = 1 (cf. Figure 1). In addition, we consider the change of error when the time t = 0.5. Since there is no exact solution to problem (1.3), we make a comparison between the solution of (3.2) on coarse mesh and the fine mesh. We choose ∆t = 0.1, 0.05, 0.025, 0.0125, respectively to solve (3.2), set u min N (x, 0.5) as the solution for ∆t min = 0.001 and denote Then the error is shown in Table 1. In Table 1, the third column err (0.5,∆t) ∆t is monotone decreasing along with the time step's waning and the fourth column err (0.5,∆t) (∆t) 2 is not monotone de-  creasing along with the time step's waning. Then the order of convergence for time belongs O(∆t) and O(∆t 2 ). It is easy to see that the result of numerical analysis on time is better than theoretical result. The reason may be the exist of nonlinear term or the limit of theoretical proof tool. Now, we consider the error for difference h at t = 0.5. We choose h = π 40 , π 80 , π 160 , π 320 , respectively to solve (3.2), set u min N (x, 0.5) as the solution for h = π 400 , ∆t = 1 100 , denote Then the error is shown in Table 2.
In Table 2, it is easy to see that the fourth column err (h,0.5) h 2 is monotone decreasing along with space step's waning. The fifth column err (h,0.5) h 3 is not monotone increase along with space step's waning, and it tends to a positive constant when the space subdivision is small enough. Hence, we can find a positive constant C, such that err (h,0.5) h 3 ≤ C, which means the order of error estimates is O(h 3 ).
It is easy to see that the third column err (h,∆t,0.5) h 3 +∆t is tend to a positive constant. Hence, we can find a positive constant C, such that err (h,∆t,0.5) h 3 +∆t ≤ C, which means the order of error estimates is O(h 3 + ∆t).