A B-spline finite element method for solving a class of nonlinear parabolic equations modeling epitaxial thin-film growth with variable coefficient

*Correspondence: hwenzhu@gmc.edu.cn 4School of Biology and Engineering, Guizhou Medical University, Guiyang, P.R. China Full list of author information is available at the end of the article Abstract In this paper, we propose an efficient B-spline finite element method for a class of fourth order nonlinear differential equations with variable coefficient. For the temporal discretization, we choose the Crank–Nicolson scheme. Boundedness and error estimates are rigorously derived for both semi-discrete and fully discrete schemes. A numerical experiment confirms our theoretical analysis.


Introduction
The epitaxial growth of nanoscale thin films has attracted a lot of attention in recent years [1][2][3][4][5][6]. The key reason for this concern is that compositions like YBa 2 Cu 3 O 7-δ (YBCO) are expected to be high temperature superconducting materials that can be used in semiconductor design. King et al. [1] proved the existence, uniqueness, and regularity of solution in an appropriate function space for the initial boundary value problem of the epitaxial thin-film growth. Kohn et al. [3] considered a fourth order parabolic equation, which is a specific example of energy-driven coarsening in two dimension space, and proved that the time-averaged energy per unit area decays no faster than t - 1 3 . The finite element method (FEM) plays an important role in solving differential equations [7][8][9][10][11]. There are some papers which have already been published to study the FEM for fourth order nonlinear parabolic equations [12][13][14][15][16][17][18]. Liu et al. [12] considered a nonlinear model describing epitaxial thin-film growth with constant coefficient and demonstrated that the Hermite FEM has the convergence rate of O( t + h 3 ). Choo [13] constructed a finite element scheme for the viscous Cahn-Hilliard equation with a nonconstant gradient energy coefficient and obtained the error estimate using the extended Lax-Richtmyer equivalence theorem. In [18], Qiao et al. presented a mixed FEM for the molecular beam epitaxy mode and showed that the semi-discrete and fully discrete schemes satisfy the nonlinear energy stability property. Moreover, the authors gave the error analysis.
In 1946, the B-spline method was first introduced by Schoenberg [19]. In 1966, Curry and Schoenberg [20] presented one element B-spline functions. In 1976, B-splines were extended to multiple situations [21]. As a class of piecewise polynomials, B-splines are often used in finite element analysis [22][23][24][25][26][27]. Erfanian et al. [25] used the linear B-spline FEM and cubic B-spline FEM for solving linear Volterra integro-differential equation in the complex plane. Dhawan et al. [26] applied the linear and quadratic B-spline functions to the advection-diffusion equations.
The main advantages of B-splines are the freedom to choose the order and smoothness, the simple data structure with one parameter in [0, 1], and the exact representation of boundary conditions. Compared with Lagrange and Hermite type elements, B-spline shape functions involve only one type of basis function. Thus the scale of matrix from B-spline FEM is smaller than that from Lagrange and Hermite elements. Moreover, Bspline shape functions are smoother. It is known that quadratic B-splines, which are in C 1 (-∞, +∞), satisfy the weak form of fourth order differential equations. However, to deal with boundary conditions, the B-spline basis functions need to be modified. In the present work, we choose the cubic B-spline FEM for a fourth order nonlinear parabolic equation with variable coefficient. It is proved that the convergence order of the Crank-Nicolson scheme is higher than that of the backward Euler scheme in [12].
The following sections are organized as follows. In Sect. 2, we introduce the model and some basic preliminaries. In Sect. 3, we show the boundedness and error estimates for the semi-discrete scheme. In Sect. 4, a fully discrete scheme based on the Crank-Nicolson method is studied. In Sect. 5, a numerical experiment is provided to confirm theoretical results.
In this work, we denote L 2 , L k , L ∞ , H k norms in I by · , · L k , | · | ∞ , and · k , respectively.

Initial boundary value problem and some preliminaries
In this paper, we consider the following problem: where I = [0, 1] and u t = ∂u ∂t . For the variable coefficient, the following assumptions hold: where s, S, M 1 , and M 2 are positive constants. Considering the boundary value conditions, we need the following space: The weak formulation associated with problem (1) is: where Du = ∂u ∂x . According to [12], the solution of problem (1) exists.

Semi-discrete finite element scheme
Let L be a positive integer. Define a uniform partition To deal with boundary value conditions, we define six additional knots The cubic B-spline function with integer knots can be defined as follows: then it is easy to get the cubic B-spline in the interval [x i , x i+4 ] which is The modified basis functions are defined as follows [28]: For the sake of convenience, we denote the modified basis functions with respect to {x i } by {ϕ i (x)}, which satisfies the following properties: The modified basis functions can deal with homogeneous as well as non-homogeneous boundary conditions.
Let U h be the cubic B-spline space. One can see that the cubic B-spline space is in where δ i (t) are time-dependent quantities. The semi-discrete finite element scheme based on B-splines for problem (5) is: Find The bandwidth of stiffness matrix is 7, and the matrix order is L -1, which is only half of the Lagrange and Hermite finite element scheme.
In order to estimate the errors of the B-spline FEM, we introduce the elliptic projection R h u: then R h u is uniquely defined, and a(u, u) ≥ C 0 u 2 2 , ∀u ∈ H 2 0 (I), where C 0 is a positive constant depending on α(x, t). Hence, a(u, v) is a symmetrical positive definite bilinear form and (see [12]) We shall discuss the boundedness of the semi-discrete scheme, which is important for error analysis.  (6) such that where C is a positive constant depending on α(x, t) and T, independent of mesh size h.
Proof According to ordinary differential equation theory, there exists a unique local solution to problem (5) in the interval [0, t n ). If we have (10), then according to the extension theorem, we can also obtain the existence of unique global solution. So, we only need to prove (10).
Taking v h = u h in (6), based on (3), we get By the interpolation inequality, we have Therefore According to the method of solving separable equations, we can get the result Integrating (12) with respect to t, we have Integrating (11) with respect to t, we get By (13), we obtain Choose v h = u h,t in (6) to get A direct calculation gives

Define the energy function
By (14), we have Integrating (17) with respect to t, we have It is obvious that In view of (3) and (16), we obtain It is clear to see from Cauchy's inequality that Therefore It follows that where C is a positive constant depending on α(x, t) and u h (0). Owing to the interpolation inequality, we obtain Thus (10) holds. The proof of the theorem is completed. Now, we give the error estimates between the solution to problem (5) and the solution in L 2 and H 2 norms. Theorem 3.2 Let u be the solution to (5), u h be the solution to (6), u(0) ∈ H 4 (I), u, u t ∈ L 2 (0, T; H 4 (I)), the initial value satisfies As 0 ≤ t ≤ T, we have the following error estimate: where C is a positive constant depending on α(x, t) and T, independent of mesh size h.

Theorem 3.3 Let u be the solution to
Then we have the following error estimate: where C is a positive constant depending on α(x, t), independent of mesh size h. (22), we get Using the triangle inequality and Sobolev's embedding theorem, we get Based on (4) and ε-inequality, we have Integrating (32) with respect to t, we find By (3), we have Then Note that Combining (33) and (34), we have Finally, using (29), we obtain (30). The proof is completed.

Fully discrete finite element scheme
To construct the Crank-Nicolson scheme, we define the following function: where H(Du n h ) is a double well potential function. Obviously, H (Du h ) = |Du h | 2 Du h -Du h . The fully discrete finite element scheme for problem (1) is: Find u n h ∈ U h (n = 1, 2, . . . , N) such that where N is a given positive integer, t = T/N denotes the time step size, t n = n t and Firstly, we analyze the boundedness of the fully discrete scheme (36). It is a key step for deducing the error estimate.
where C is a positive constant depending on α(x, t) and T, independent of h and t.
Letting γ = 1 2s , we have It is easy to show If t is small enough, we conclude Then we get Define the function then G(u n h , t n-1 2 ) ≥ 0. By (44) and (45), we have G u n h , t n-1 2 ≤ G u n-1 h , t n-3 2 + 1 2 α x, t n-1 2α x, t n-3 2 D 2 u n-1 h 2 , 1 .
With the differential mean value theorem and the boundedness of variable coefficient, we obtain Taking the sum over n, we get It is obvious that Therefore we know Based on (44) and u 0 h ∈ H 2 0 (I) ∩ W 1,4 (I), we have Using discrete Gronwall's inequality, we derive Based on (48), it is easy to see We also know By (42) and (49), we obtain (37). The proof is completed.

Theorem 4.2 Let u n be the solution to problem
Then we have the following error estimate: where C is a positive constant depending on α(x, t) and T, independent of mesh size h.
Proof Denote u n t = u t (x, t n ) and u n = u(x, t n ). Setting t = t n-1 and t = t n in (5), respectively, we obtain Let ρ n = u n -R h u n and θ n = R h u nu n h , then u nu n h = ρ n + θ n . It is clear to get where An easy calculation gives Using Taylor's theorem, we have α x, t n = α x, t n-1 2 + t 2 ∂α ∂t x, t n-1 2 + ( t) 2 6 With (4), we get From (7), we have Setting v h = θ n + θ n-1 in (58), we get 1 t θ n 2θ n-1 2 + s 2 D 2 θ n + D 2 θ n-1 2 ≤ r n 2 + 1 4 θ n + θ n-1 2 + F Du n , Du n-1 , Du n h , Du n-1 Based on the Newton-Leibniz formula and Hölder's inequality, we have A direct calculation gives Using Hölder's inequality, we have In view of (62)-(64) and (9) Taking the sum over n, by n t = t n ≤ T, we have If t is small enough, we have By discrete Gronwall's inequality, it gives θ n ≤ C ( t) 2 + h 3 .
In the following theorem, we introduce the error estimate in H 2 norm.

Theorem 4.3
Let u n be the solution to (5), u n h be the solution to the fully discrete problem (36), u(0) ∈ H 4 (I), u t ∈ L 2 (0, T; H 4 (I)) ∩ L 2 (0, T; W 2,4 (I)), u ttt ∈ L 2 (0, T; L 2 (I)), and u 0 h ∈ U h satisfying Then we have the following error estimate: Taking the sum over n and using (4), we can obtain It follows from (3) that Then one has First, applying the triangle inequality to I 1 , we get

Numerical approximation
In this section, a numerical example is provided to illustrate the proposed B-spline FEM for solving the nonlinear parabolic equation. The efficiency of the cubic B-spline finite element scheme is tested. We consider where α(x, t) = 1 + xt is selected to satisfy the primary assumptions. We take the analytical solution u(x, t) = t 2 (1 -cos 2πx). Then the concrete functional form of f (x, t) is f (x, t) = 2t + 4π 2 t 2 -2 1 + 8π 4 x t -16π 4 cos 2πx -16π 3 sin 2πx -48π 4 sin 2 πxcos2πx. Figure 1 illustrates the behavior of the exact solution to problem (81), and Fig. 2 demonstrates the profile of the solution to the fully discrete scheme.
In this example, the numerical solution is in good accordance with the exact solution, indicating that the numerical scheme is valid and efficient.
Then choosing t = 1, the corresponding errors and convergence rates of the cubic Bspline FEM are shown in Tables 1-3. In Table 1, to analyze the spatial convergence order, we take the time step t = 1 8000 . The values in Table 1 indicate that with the decreasing of the space size, the error is monotone decreasing. We also find that the numerical solution to the scheme is fourth order convergent in L 2 norm and is second order convergent in H 2 norm.
In Table 2, we consider the error estimates and convergence orders in time direction when the space step is fixed to h = 1 1000 . It is easy to see that the orders of error estimate both are second order in L 2 and H 2 norms.
In Table 3, we analyze the convergence rate when the space and time step change at the same time. We choose ( t, h) = ( 1 100 , 1 10 ), ( 1 400 , 1 20 ), ( 1 1600 , 1 40 ), ( 1 6400 , 1 80 ), respectively. One can see that the relation between the space step and the time step is t/h 2 = 1. This shows that the B-spline finite element scheme is very stable.    The numerical experiment indicates that the cubic B-spline FEM is an efficient approximation tool for solving the fourth order nonlinear parabolic equation.

Conclusion
In this paper, we propose the B-spline FEM for a class of fourth order nonlinear parabolic equations. On the one hand, B-splines have better smoothness than the Lagrange and Hermite type elements. On the other hand, B-spline finite element only has one type of basis functions, so the scale of matrix from B-spline FEM is lower.
The coefficient α(x, t) is variable, which broadens the application fields and also increases the difficulty of analysis. By defining the biharmonic projection operator and the energy function, we prove the boundedness of the semi-discrete and fully discrete schemes based on B-splines. Further, the error estimates in L 2 norm and H 2 norm are deduced by using the boundedness, Sobolev's embedding theorem, and so on. The results of numerical example confirm our theoretical analysis.
In general, the B-spline FEM is an efficient method for solving higher order nonlinear parabolic equations. By using central difference, the convergence rate in time direction, which can be improved, is of second order.