# The Spectral Method for the Cahn-Hilliard Equation with Concentration-Dependent Mobility

Authors: Shimin Chai; Yongkui Zou
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

This paper is concerned with the numerical approximations of the Cahn-Hilliard-type equation with concentration-dependent mobility. Convergence analysis and error estimates are presented for the numerical solutions based on the spectral method for the space and the implicit Euler method for the time. Numerical experiments are carried out to illustrate the theoretical analysis.


---

## Body

## 1. Introduction

In this paper, we apply the spectral method to approximate the solutions of Cahn-Hilliard equation, which is a typical class of nonlinear fourth-order diffusion equations. Diffusion phenomena is widespread in the nature. Therefore, the study of the diffusion equation caught wide concern. Cahn-Hilliard equation was proposed by Cahn and Hilliard in 1958 as a mathematical model describing the diffusion phenomena of phase transition in thermodynamics. Later, such equations were suggested as mathematical models of physical problems in many fields such as competition and exclusion of biological groups [1], moving process of river basin [2], and diffusion of oil film over a solid surface [3]. Due to the important application in chemistry, material science, and other fields, there were many investigations on the Cahn-Hilliard equations, and abundant results are already brought about.The systematic study of Cahn-Hilliard equations started from the 1980s. It was  Elliott and Zheng [4] who first study the following so-called standard Cahn-Hilliard equation with constant mobility:
(1.1)∂u∂t+div[k∇Δu-∇A(u)]=0.
Basing on global energy estimates, they proved the global existence and uniqueness of classical solution of the initial boundary problem. They also discussed the blow-up property of classical solutions. Since then, there were many remarkable studies on the Cahn-Hilliard equations, for example, the asymptotic behavior of solutions [5–8], perturbation of solutions [9, 10], stability of solutions [11, 12], and the properties of the solutions for the Cahn-Hilliard equations with dynamic boundary conditions [13–16]. In the mean time, a number of the numerical techniques for Cahn-Hilliard equations were produced and developed. These techniques include the finite element method [4, 17–24], the finite difference method [25–29], the spectral method, and the pseudospectral method [30–35]. The finite element method for the Cahn-Hilliard equation is well investigated by many researchers. For example, in [23, 36], (1.1) was discreted by conforming finite element method with an implicit time discretization. In [22], semidiscrete schemes which can define a Lyapunov functional and remain mass constant were used for a mixed formulation of the governing equation. In [37], a mixed finite element formulation with an implicit time discretization was presented for the Cahn-Hilliard equation (1.1) with Dirichlet boundary conditions. The conventional strategy to obtain numerical solutions by the finite differece method is to choose appropriate mesh size based on the linear stability analysis for different schemes. However, this conventional strategy does not work well for the Cahn-Hilliard equation due to the bad numerical stability. Therefore, an alternative strategy is proposed for general problems, for example, in [26, 38] the strategy was to design such that schemes inherit the energy dissipation property and the mass conservation by Furihata. In [27, 28], a conservative multigrid method was developed by Kim.The advantage of the spectral method is the infinite order convergence; that is, if the exact solution of the Cahn-Hilliard equation isCk smooth, the approximate solution will be convergent to the exact solution with power for exponent N-k, where N is the number of the basis function. This method is superior to the finite element methods and finite difference methods, and a lot of practice and experiments convince the validity of the spectral method [39]. Many authors have studied the solution of the Cahn-Hilliard equation which has constant mobility by using spectral method. For example, in [33–35], Ye studied the solution of the Cahn-Hilliard equation by Fourier collocation spectral method and Legendre collocation spectral method under different boundary conditions. In [30], the author studied a class of the Cahn-Hilliard equation with pseudospectral method. However, the Cahn-Hilliard equation with varying mobility can depict the physical phenomena more accurately; therefore, there is practical meaning to study the numerical solution for the Cahn-Hilliard equation with varying mobility. Yin [40, 41] studied the Cahn-Hilliard equation with concentration-dependent mobility in one dimension and obtained the existence and uniqueness of the classical solution. Recently, Yin and Liu [42, 43] investigate the regularity for the solution in two dimensions. Some numerical techniques for the Cahn-Hilliard equation with concentration-dependent mobility are already studied with the finite element method [28] and with finite difference method [44].In this paper, we consider an initial-boundary value problem for Cahn-Hilliard equation of the following form:(1.2)∂u∂t+D[m(u)(D3u-DA(u))]=0,(x,t)∈(0,1)×(0,T),(1.3)Du|x=0,1=D3u|x=0,1=0,t∈(0,T),(1.4)u(x,0)=u0,x∈(0,1),
where D=∂/∂x and
(1.5)A(s)=-s+γ1s2+γ2s3,γ2>0.
Here, u(x,t) represents a relative concentration of one component in binary mixture. The function m(u) is the mobility which depends on the unknown function u, which restricts diffusion of both components to the interfacial region only. Denote QT=(0,1)×(0,T). Throughout this paper, we assume that
(1.6)0<m0≤m(s)≤M0,|m′(s)|≤M1,∀s∈ℝ,
where m0, M0, and M1 are positive constants. The existence and uniqueness of the classical solution of the problems (1.2)–(1.4) were proved by Yin [41].In this paper, we will apply the spectral method to discretize the spatial variables of (1.2) to construct a semidiscrete system. We prove the existence and boundedness of the solutions of this semidiscrete system. Then, we apply implicit midpoint Euler scheme to discretize the time variable and obtain a full-discrete scheme, which inherits the energy dissipation property. The property of the mobility depending on the solution of (1.2) causes much troubles for the numerical analysis. Furthermore, with the aid of Nirenberg inequality we investigate the boundedness and convergence of the numerical solutions of the full-discrete equations. We also obtain the error estimation for the numerical solutions to the exact ones.This paper is organized as follows. In Section2, we study the spectral method for (1.2)–(1.4) and obtain the error estimate between the exact solution u and the spectral approximate solution uN. In Section 3, we use the implicit Euler method to discretize the time variable and obtain the error estimate between the exact solution u and the full-discrete approximate solution UNj. Finally in Section 4, we present a numerical computation to illustrate the theoretical analysis.


## 2. Semidiscretization with Spectral Method

In this section, we apply the spectral method to discretize (1.2)–(1.4) and study the error estimate between the exact solution and the semidiscretization solution.Denote by∥·∥k and |·|k the norm and seminorm of the Sobolev spaces Hk(0,1)(k∈ℕ), respectively. Let (·,·) be the standard L2 inner product over (0,1). Define
(2.1)HE2(0,1)={v∈H2(0,1);Dv|x=0,1=0}.A functionu is said to be a weak solution of the problems (1.2)–(1.4), if u∈L∞(0,T;HE2(0,1)) and satisfies the following equations:
(2.2)(∂u∂t,v)+(D2u-A(u),D(m(u)Dv))=0,∀v∈HE2,(2.3)(u(·,0),v)-(u0,v)=0,∀v∈HE2.For any integerN>0, let SN=span{coskπx,k=0,1,2,…,N}. Define a projection operator PN:HE2↦SN by
(2.4)∫01(PN)u(x)v(x)dx=∫01u(x)v(x)dx,∀v∈SN.
We collect some properties of this projection PN in the following lemma (see [39]).Lemma 2.1.
(i)PN commutes with the second derivation on HE2(I), that is,
(2.5)PND2v=D2PNv,∀v∈HE2(I).
(ii) For any0≤μ≤σ, there exists a positive constant C such that
(2.6)∥v-PNv∥μ≤CNμ-σ|v|σ,v∈Hσ(I).The following Nirenberg inequality is a key tool for our theoretical estimates.Lemma 2.2.
Assume thatΩ⊂ℝn is a bounded domain, u∈Wm,r(Ω), then we have
(2.7)∥Dju∥Lp⩽C1∥Dmu∥Lra∥u∥Lq1-a+C2∥u∥Lq,
where
(2.8)jm⩽a<1,1p=jn+a(1r-mn)+(1-a)1q.By [41], we have the following.Lemma 2.3.
Assume thatm(s)∈C1+α(ℛ), u0∈C4+α(I-), Diu0(0)=Diu0(1)=0(i=0,1,2,3,4),m(s)>0, then there exists a unique solution of the problems (1.2)–(1.4) such that
(2.9)u∈C1+α/4,4+α(Q-T).The spectral approximation to (2.2) is to find an element uN(·,t)∈SN such that
(2.10)(∂uN∂t,vN)+(D2uN-PNA(uN),D(m(uN)DvN))=0,∀vN∈SN,(2.11)(uN(·,0),vN)-(u0,vN)=0,∀vN∈SN.Now we study theL∞ norm estimates of the function uN(·,t) and DuN(·,t) for 0≤t≤T.Theorem 2.4.
Assume (1.6) and u0∈HE2. Then there is a unique solution of (2.10) and (2.11) such that
(2.12)∥uN∥∞≤C,∥DuN∥≤C,0≤t≤T,
where C=C(u0,T,γ1,γ2) is a positive constant.Proof.
From (2.11) it follows that uN(·,0)=PNu0(·). The existence and uniqueness of the initial problem follow from the standard ODE theory. Now we study the estimate.
Define an energy function:(2.13)F(t)=(H(uN),1)+12∥DuN∥2,
where H(s)=∫0s‍A(t)dt=(1/4)γ2s4+(1/3)γ1s3-(1/2)s2. Direct computation gives
(2.14)dFdt=(A(uN),∂uN∂t)+(DuN,D∂uN∂t)=(PNA(uN),∂uN∂t)-(D2uN,∂uN∂t)=(∂uN∂t,PNA(uN)-D2uN).
Noticing that PNA(uN)-D2uN∈SN and setting vN=PNA(uN)-D2uN in (2.10), applying integration by part, we obtain
(2.15)dFdt=(m(uN)(D3uN-DPNA(uN)),D(PNA(uN)-D2uN))=-(m(uN)(D3uN-DPNA(uN)),D3uN-DPNA(uN))≤-m0∥D(PNA(uN)-D2uN)∥2≤0.
Hence,
(2.16)F(t)≤F(0)=∫01H(PNu0)ds+12∥DPNu0∥2,∀0≤t≤T.
Applying Young inequality, we obtain
(2.17)uN2≤εuN4+C1ε,uN3≤εuN4+C2ε,
where C1ε and C2ε are positive constants. Letting ε=3γ2/(8|γ1|+12), then for all 0≤t≤T,
(2.18)∫01‍H(uN)dx≥∫01(14γ2uN4-13|γ1uN3|-12uN2)dx≥γ28∫01(uN)4dx-K1,
where K1 is a positive constant depending only on γ1 and γ2. Therefore, we get
(2.19)12∥DuN(·,t)∥2+γ28∫01(uN(·,t))4dx-K1≤12∥DuN(·,t)∥2+∫01H(uN(·,t))dx=F(t)≤F(0)=∫01H(PNu0)dx+∥DPNu0∥2.
Thus,
(2.20)∥DuN(·,t)∥2≤2K1+2F(0)=2K1+2∫01H(PNu0)dx+∥DPNu0∥2≜C,∫01(uN(x,t))4dx≤8γ2(K1+F(0))=8γ2(K1+∫01H(PNu0)dx+12∥DPNu0∥2)≜C,
where C=C(u0,γ1,γ2) is a constant. By Hölder inequality, we obtain
(2.21)∥uN∥2=∫01uN2dx≤(∫01uN4dx)1/2(∫0112dx)1/2=(∫01uN4dx)1/2.
Therefore,
(2.22)∥uN∥≤C,∥DuN∥≤C.
From the embedding theorem it follows that
(2.23)∥uN∥∞≤C,∀0≤t≤T.Theorem 2.5.
Assume (1.6) and let uN(·,t) be the solution of (2.10) and (2.11). Then there is a positive constant C=C(u0,m,T,γ1,γ2)>0 such that
(2.24)∥DuN∥∞≤C,∥D2uN∥≤C.Proof.
SettingvN=D4uN in (2.10) and integrating by parts, we get
(2.25)12ddt∥D2uN∥2+(m(uN)(D4uN-D2PNA(uN)),D4uN)+(m′(uN)DuN(D3uN-DPNA(uN)),D4uN)=0.
Consequently,
(2.26)12ddt∥D2uN∥2+(m(uN)D4uN,D4uN)=I1+I2+I3,
where
(2.27)I1=(m(uN)D2PNA(uN),D4uN),I2=-(m′(uN)DuND3uN,D4uN),I3=(m′(uN)DuNDPNA(uN),D4uN).
Noticing (1.6), we have
(2.28)12ddt∥D2uN∥2+m0∥D4uN∥2≤12ddt∥D2uN∥2+(m(uN)D4uN,D4uN).
In terms of the Nirenberg inequality (2.7), there is a constant C>0 such that
(2.29)∥DuN∥∞≤C(∥D4uN∥3/8∥u∥5/8+∥u∥),∥DuN∥∞≤C(∥D2uN∥3/4∥u∥1/4+∥u∥),∥D3uN∥∞≤C(∥D4uN∥7/8∥u∥1/8+∥u∥).
Noticing the definition of the function A and the estimates in (2.12), we have
(2.30)∥A(uN)∥∞≤C,∥A′(uN)∥∞≤C,∥A′′(uN)∥∞≤C,
for some constant C=C(u0,T,γ1,γ2)>0. Applying Hölder inequality and Young inequality, for any ε>0,
(2.31)|I1|≤M0(∥A′(uN)D2uN∥+∥A′′(uN)(DuN)2∥)·∥D4uN∥≤M0∥A′(uN)∥∞·∥D2uN∥∥D4uN∥+M0∥A′′(uN)∥∞·∥DuN∥∥DuN∥∞∥D4uN∥≤ε∥D4uN∥2+C∥D2uN∥2+C(∥D4uN∥3/8+1)∥D4uN∥≤ε∥D4uN∥2+C∥D2uN∥2+C∥D4uN∥11/8+C∥D4uN∥≤ε∥D4uN∥2+C∥D2uN∥2+ε2∥D4uN∥2+C+ε2∥D4uN∥2+C≤2ε∥D4uN∥2+C∥D2uN∥2+C,
where C=C(u0,m,T,γ1,γ2,ε)>0 is a positive constant. Similarly, we obtain
(2.32)|I2|≤M1∥DuN∥∥D3uN∥∞∥D4uN∥≤C(∥D4uN∥7/8+1)∥D4uN∥≤C(∥D4uN∥15/8+∥D4uN∥)≤ε2∥D4uN∥2+C+ε2∥D4uN∥2+C≤ε∥D4uN∥2+C,|I3|≤M1∥DPNA(uN)∥∞·∥DuN∥·∥D4uN∥≤C∥DPNA(uN)∥∞∥D4uN∥≤C(∥D2PNA(uN)∥3/4∥PNA(uN)∥1/4+∥PNA(uN)∥)∥D4uN∥≤C(∥D2PNA(uN)∥3/4∥A(uN)∥1/4+∥A(uN)∥)∥D4uN∥≤C(∥D2PNA(uN)∥3/4∥A(uN)∥∞1/4+∥A(uN)∥∞)∥D4uN∥≤C(∥D2PNA(uN)∥3/4+1)∥D4uN∥≤C(∥D2PNA(uN)∥+1)∥D4uN∥≤C(∥D2A(uN)∥+1)∥D4uN∥≤C(∥A′(uN)D2uN∥+∥A′′(uN)(DuN)2∥)·∥D4uN∥≤ε∥D4uN∥2+C.
Hence,
(2.33)12ddt∥D2uN∥2+m0∥D4uN∥2≤4ε∥D4uN∥2+C∥D2uN∥2+C.
Taking ε=m0/8, we have
(2.34)12ddt∥D2uN∥2+m02∥D4uN∥2≤C∥D2uN∥2+C,
where C=C(u0,m,T,γ1,γ2)>0 is a positive constant. Therefore,
(2.35)12ddt∥D2uN∥2≤C∥D2uN∥2+C.
From Gronwall inequality it follows that
(2.36)∥D2uN∥2≤expCt(∥D2u0∥2+1)≤expCT(∥D2u0∥2+1)≤C,
where C=C(u0,m,T,γ1,γ2)>0 is a positive constant. According to the embedding theorem, we have
(2.37)∥DuN∥∞≤C,0≤t≤T,
where C=C(u0,m,T,γ1,γ2)>0 is a positive constant.Now, we study the error estimates between the exact solutionu and the semidiscrete spectral approximation solution uN. Set the following decomposition:
(2.38)u-uN=η+e,η=u-PNu,e=PNu-uN.
From the inequality (2.6) it follows that
(2.39)∥η∥≤CN-4.
Hence, it remains to obtain the approximate bounds of e.Theorem 2.6.
Assume thatu is the solution of (1.2)–(1.4), uN∈SN is the solution of (2.10) and (2.11), and m is smooth and satisfies (1.6), then there exists a constant C=C(u0,m,γ1,γ2) such that
(2.40)∥e∥≤C(N-2+∥e(0)∥).Before we prove this theorem, we study some useful approximation properties.Lemma 2.7.
For anyε>0, we have
(2.41)-(D2η+D2e,D(m(uN)De))≤-m0∥D2e∥2+3ε∥D2e∥2+C(∥e∥2+N-4),
where C=C(u0,m,T,γ1,γ2,ε)>0 is a positive constant.Proof.
Direct computation gives(2.42)-(D2η+D2e,D(m(uN)De))=-(D2e,m(uN)D2e)-(D2η,m(uN)D2e)-(D2η,m′(uN)DuNDe)-(D2e,m′(uN)DuNDe)≤-m0∥D2e∥2+M0∥D2e∥∥D2η∥2+M1∥D2η∥∥DuN∥∞∥De∥+M1∥D2e∥∥DuN∥∞∥De∥≤-m0∥D2e∥2+ε∥D2e∥2+M024ε∥D2η∥2+M1∥DuN∥∞(∥De∥2+∥D2η∥2)+ε∥D2e∥2+M12∥DuN∥∞24ε∥De∥2≤-m0∥D2e∥2+2ε∥D2e∥2+(M024ε+M1∥DuN∥∞)∥D2η∥2+(M1∥DuN∥∞+M12∥DuN∥∞24ε)∥De∥2≤-m0∥D2e∥2+2ε∥D2e∥2+(M024ε+M1∥DuN∥∞)∥D2η∥2+(M1∥DuN∥∞+M12∥DuN∥∞24ε)∥D2e∥·∥e∥≤-m0∥D2e∥2+2ε∥D2e∥2+C∥D2η∥2+ε∥D2e∥2+C∥e∥2≤-m0∥D2e∥2+3ε∥D2e∥2+C(∥e∥2+∥D2η∥2)≤-m0∥D2e∥2+3ε∥D2e∥2+C(∥e∥2+N-4),
where C=C(u0,m,T,γ1,γ2,ε)>0 is a positive constant.Lemma 2.8.
For anyε>0, we have
(2.43)(A(u)-PNA(uN),D(m(uN)De))≤3εM02∥D2e∥2+C(∥e∥2+N-4),
where C=C(u0,m,T,γ1,γ2,ε)>0 is a positive constant.Proof.
Noticing that(2.44)(A(u)-PNA(uN),D(m(uN)De))=(A(u)-PNA(u),D(m(uN)De))+(PNA(u)-PNA(uN),D(m(uN)De))≤∥A(u)-PNA(u)∥·∥D(m(uN)De)∥+∥PN(A-(u,uN)(u-uN))∥·∥D(m(uN)De)∥≤∥A(u)-PNA(u)∥·∥D(m(uN)De)∥+∥A-(u,uN)(η+e)∥·∥D(m(uN)De)∥,
where
(2.45)A-(s,τ)=γ2(s2+sτ+τ2)+γ1(s+τ)-1.
From the boundedness of uN in Theorem 2.4 and the property of u in Lemma 2.3, it follows that
(2.46)∥A-(u,uN)∥∞≤C.
Then we obtain
(2.47)∥A(u)-PNA(u)∥2≤C(∥u3-PNu3∥2+∥u2-PNu2∥2+∥u-PNu∥2)≤CN-8(|u3|4+|u2|4+|u|4)2≤CN-8,∥A-(u,uN)(η+e)∥2≤∥A-(u,uN)∥∞2(∥e∥2+∥η∥2)≤C(N-8+∥e∥2),∥D(m(uN)De)∥2≤2∥m(uN)D2e∥2+2∥m′(uN)DuNDe∥2≤2M02∥D2e∥2+2M12∥DuN∥∞2(∥D2e∥2∥e∥2)≤2M02∥D2e∥2+M02∥D2e∥2+M14∥DuN∥∞4M02∥e∥2≤3M02∥D2e∥2+C∥e∥2,
where C=C(u0,m,T,γ1,γ2)>0 is a positive constant. By Cauchy inequality, for any ε>0, we have
(2.48)(A(u)-PNA(uN),D(m(uN)De))≤∥A(u)-PNA(u)∥·∥D(m(uN)De)∥+∥A-(u,uN)(η+e)∥·∥D(m(uN)De)∥≤ε2∥D(m(u)De)∥2+12ε∥A(u)-PNA(u)∥2+ε2∥D(m(u)De)∥2+12ε∥A-(u,uN)(η+e)∥2≤3εM02∥D2e∥2+C(∥e∥2+N-8),
where C=C(u0,m,T,γ1,γ2,ε)>0 is a positive constant.Lemma 2.9.
Assume thatu is the solution of (1.2)–(1.4), there exists a positive constant C=C(u0,m,γ1,γ2) such that
(2.49)(D2u-A(u),D(m(u)-m(uN))De)≤m02∥D2e∥2+C(∥e∥2+N-4),
where C=C(u0,m,T,γ1,γ2)>0 is a positive constant.Proof.
By Lemma2.3, we have
(2.50)∥D3u-DA(u)∥∞≤∥D3u∥∞+∥A′(u)Du∥∞≤C.
In the other hand, it follows that
(2.51)∥m(u)-m(uN)∥2≤M12∥u-uN∥2≤M12(∥e∥2+∥η∥2).
By the Young inequality,
(2.52)(D2u-A(u),D(m(u)-m(uN)De))=-(D3u-DA(u),(m(u)-m(uN))De)≤∥D3u-DA(u)∥∞∥m(u)-m(uN)∥·∥De∥≤C∥m(u)-m(uN)∥∥De∥≤ε∥De∥2+Cε(∥e∥2+∥η∥2)≤ε∥D2e∥2+Cε(∥e∥2+∥η∥2).
Choosing ε=m0/2 in the previous inequality, we obtain (2.49).Proof of Theorem2.6.
Settingv=e in (2.2), we obtain
(2.53)(∂u∂t,e)+(D2u-A(u),D(m(u)De))=0.
Setting vN=e in (2.10), we get
(2.54)(∂uN∂t,e)+(D2uN-PNA(uN),D(m(uN)De))=0.
(2.53) minus (2.54) gives
(2.55)(∂e∂t,e)=-(Du2-A(u),D(m(u)-m(uN)De))-(D2u-D2uN,D(m(uN)De))+(A(u)-PNA(uN),D(m(uN)De)).
According to Lemmas 2.7, 2.8 and 2.9, we have
(2.56)12ddt∥e∥2≤-m0∥D2e∥2+3ε∥D2e∥2+C(∥e∥2+N-4)+3εM02∥D2e∥2+C(N-8+∥e∥2)≤(-m0+(4+3M02)ε)∥D2e∥2+C(N-4+∥e∥2).
Set ε=m0/(4+3M02), then there exists a positive constant C=C(u0,m,γ1,γ2) such that
(2.57)ddt∥e∥2≤C(∥e∥2+N-4).
By Gronwall inequality, we have
(2.58)∥e∥≤expCt(∥e0∥+N-2)≤expCT(∥e0∥+N-2)≤C(∥e0∥+N-2),
where C=C(u0,m,T,γ1,γ2)>0 is a constant.Summig up the properties above, we obtain the following.Theorem 2.10.
Assumem(s) is sufficiently smooth and satisfies (1.6), u is the solution of (1.2)–(1.4), and uN is the solution of (2.10) and (2.11). Then there exists a positive constant C=C(u0,m,T,γ1,γ2)>0 such that
(2.59)∥u-uN∥≤C(N-2+∥u0-u0N∥),∀t∈(0,T).


## 3. Full-Discretization Spectral Scheme

In this section, we apply implicit midpoint Euler scheme to discretize time variable and get a full-discrete form. Furthermore, we investigate the boundedness of numerical solution and the convergence of the numerical solutions of the full-discrete system. We also obtain the error estimates between the numerical solution and the exact ones.Firstly, we introduce a partition of[0,T]. Let 0=t0<t1<⋯<tΛ, where tj=jh and h=T/Λ is time-step size. Then the full-discretization spectral method for (1.2)–(1.4) reads: ∀v∈SN, find UNj∈SN(j=0,1,2,…,Λ) such that
(3.1)(UNj+1-UNjh,v)+(D2U-Nj+1/2-PNA~(UNj+1,UNj),D(m(U-Nj+1/2)Dv))=0,(3.2)(UN0,v)-(u0,v)=0,
where U-Nj+1/2=(UNj+1+UNj)/2 and
(3.3)A~(ϕ,φ)=γ24(ϕ3+φ3+ϕ2φ+ϕφ2)+γ13(ϕ2+ϕφ+φ2)-12(ϕ+φ).
The solution UNj has the following property.Lemma 3.1.
Assume thatUNj∈SN(j=1,2,…,Λ) is the solution of (3.1)-(3.2). Then there exists a constant C=C(u0,m,γ1,γ2)>0 such that
(3.4)∥UNj∥∞≤C,∥DUNj∥≤C.Proof.
Define a discrete energy function at timetj by
(3.5)F(j)=12∥DUNj∥2+(H(UNj),1).
Notice that
(3.6)1h(F(j+1)-F(j))=(DU-Nj+1/2,DUNj+1-UNjh)+(PNA~(UNj+1,UNj),UNj+1-UNjh)=-(D2U-Nj+1/2-PNA~(UNj+1,UNj),UNj+1-UNjh).
Setting v=D2U-Nj+1/2-PNA~(UNj+1,UNj)∈SN in (3.1), we obtain
(3.7)1h(F(j+1)-F(j))=-(m(U-Nj+1/2)(D3U-Nj+1/2-DPNA~(UNj+1,UNj)),D3U-Nj+1/2-DPNA~(UNj+1,UNj))≤-m0∥D3U-Nj+1/2-DPNA~(UNj+1,UNj)∥2≤0,
which implies
(3.8)F(j)≤F(0)=12∥DUN0∥2+(H(UN0),1).
By (2.18), we have
(3.9)∫01‍H(UNj)dx≥γ28∫01‍(UNj)4dx-K1.
Then
(3.10)12∥DUNj∥2+γ28∫01‍(UNj)4dx-K1≤12∥DUNj∥2+∫01‍H(UNj)dx=F(j)≤F(0)=12∥DUN0∥2+(H(UN0),1).
So we obtain
(3.11)∥DUNj∥2≤2K1+∥DUN0∥2+2(H(UN0),1)≡C1,∫01‍(UNj)4dx≤8γ2(K1+12∥DUN0∥2+(H(UN0),1))≡C2,
where C1=C1(u0,m,γ1,γ2) and C2=C2(u0,m,γ1,γ2) are positive constants. By Hölder inequality, we get
(3.12)∥UNj∥2=∫01‍(UNj)2dx≤(∫01‍(UNj)4dx)1/2≤C21/2.
Therefore,
(3.13)∥UNj∥≤C21/4.
By the embedding theorem, we obtain
(3.14)∥UNj∥∞≤C,
where C=C(u0,m,γ1,γ2)>0 is a constant.Lemma 3.2.
Assume thatUNj∈SN(j=1,2,…,Λ) is the solution of the full-discretization scheme (3.1)-(3.2), then there is a constant C=C(u0,m,T,γ1,γ2)>0 such that
(3.15)∥DUNj∥∞≤C,∥D2UNj∥≤C.Proof.
Settingv=2D4U-Nj+1/2∈SN in (3.1), we have
(3.16)(UNj+1-UNjh,2D4U-Nj+1/2)+(D(m(U-Nj+1/2)D3U-Nj+1/2-DPNA~(UNj+1,UNj))2D4U-Nj+1/2)=(D2UNj+1-D2UNjh,D2UNj+1+D2UNj)+(m(U-Nj+1/2)D4U-Nj+1/2,2D4U-Nj+1/2)+(m′(U-Nj+1/2)DU-Nj+1/2D3U-Nj+1/2,2D4U-Nj+1/2)-(m(U-Nj+1/2)D2PNA~(UNj+1+UNj),2D4U-Nj+1/2)-(m′(U-Nj+1/2)DU-Nj+1/2DPNA~(UNj+1,UNj),2D4U-Nj+1/2)=0.
Therefore,
(3.17)1h(∥D2UNj+1∥2-∥D2UNj∥2)+2(m(U-Nj+1/2)D4U-Nj+1/2,D4U-Nj+1/2)=2(m′(U-Nj+1/2)D3U-Nj+1/2DU-Nj+1/2,D4U-Nj+1/2)+2(m(U-Nj+1/2)D2PNA~(UNj+1,UNj),D4U-Nj+1/2)+2(m′(U-Nj+1/2)DU-Nj+1/2DPNA~(UNj+1,UNj+1),D4U-Nj+1/2)=I1j+I2j+I3j.
By Nirenberg inequality (2.7), we have
(3.18)∥DUNj∥∞≤C(∥D4UNj∥3/8∥UNj∥5/8+∥UNj∥),∥DUNj∥∞≤C(∥D2UNj∥3/4∥UNj∥1/4+∥UNj∥),∥D3UNj∥∞≤C(∥D4UNj∥7/8∥UNj∥1/8+∥UNj∥).
According to (3.14), we obtain
(3.19)∥A~1′(UNj+1,UNj)∥∞≤C,∥A~2′(UNj+1,UNj)∥∞≤C,∥A~11′′(UNj+1,UNj)∥∞≤C,∥A~12′′(UNj+1,UNj)∥∞≤C,∥A~22′′(UNj+1,UNj)∥∞≤C,
where C=C(u0,m,γ1,γ2) is a positive constant. By Young inequality, for any positive constant ε>0, it follows that
(3.20)|I1j|≤2M1∥D3U-Nj+1/2∥∞∥DU-Nj+1/2∥∥D4U-Nj+1/2∥≤C(∥D4U-Nj+1/2∥7/8+1)∥D4U-Nj+1/2∥≤C∥D4U-Nj+1/2∥15/8+C∥D4U-Nj+1/2∥≤ε2∥D4U-Nj+1/2∥2+Cε+ε2∥D4U-Nj+1/2∥2+Cε≤ε∥D4U-Nj+1/2∥2+Cε,|I2j|≤2M0∥D2A~(UNj+1,UNj)∥∥D4U-Nj+1/2∥≤2M0∥A11′′(DUNj+1)2+2A12′′DUNj+1DUNj+A22′′(DUNj)22M02+A'1D2UNj+1+A'2D2UNj(DUNj+1)2∥∥D4U-Nj+1/2∥≤C(∥D2UNj+1∥+∥D2UNj∥+∥(DUNj+1)2∥+∥(DUNj)2∥)∥D4U-Nj+1/2∥≤C(∥D2UNj+1∥+∥D2UNj∥+∥DUNj+1∥∥DUNj+1∥∞+∥DUNj∥∥DUNj∥∞)∥D4U-Nj+1/2∥≤C(∥D2UNj+1∥+∥D2UNj∥+∥DUNj+1∥∞+∥DUNj∥∞)∥D4U-Nj+1/2∥≤C(∥D2UNj+1∥+∥D2UNj∥+∥D2UNj+1∥3/4+∥D2UNj∥3/4+1)∥D4U-Nj+1/2∥≤C(∥D2UNj+1∥+∥D2UNj∥+1)∥D4U-Nj+1/2∥≤ε∥D4U-Nj+1/2∥2+Cε(∥D2UNj+1∥2+∥D2UNj∥2+1),|I3j|≤M1∥DU-Nj+1/2∥∥DPNA~(UNj+1,UNj)∥∞∥D4U-Nj+1/2∥≤C(∥D2PNA~(UNj+1,UNj)∥3/4+1)∥D4U-Nj+1/2∥≤C(∥D2PNA~(UNj+1,UNj)∥+1)∥D4U-Nj+1/2∥≤C(∥D2A~(UNj+1,UNj)∥∥D4U-Nj+1/2∥+∥D4U-Nj+1/2∥)≤ε∥D4U-Nj+1/2∥2+Cε(∥D2UNj+1∥2+∥D2UNj∥2+1)+ε∥D4U-Nj+1/2∥2+Cε≤2ε∥D4U-Nj+1/2∥2+Cε(∥D2UNj+1∥2+∥D2UNj∥2+1),
where Cε=Cε(u0,m,γ1,γ2)>0 is a constant. Therefore,
(3.21)1h(∥D2UNj+1∥2-∥D2UNj∥2)+2m0∥D4U-Nj+1/2∥2≤1h(∥D2UNj+1∥2-∥D2UNj∥2)+2(m(U-Nj+1/2)D4U-Nj+1/2,D4U-Nj+1/2)≤4ε∥D2U-Nj+1/2∥2+Cε(∥D2UNj+1∥2+∥D2UNj∥2+1).
Setting ε=m0/2, there is a positive constant C=C(u0,m,γ1,γ2) such that
(3.22)∥D2UNj+1∥2≤(1+2Ch1-Ch)∥D2UNj∥2+C1-Chh.
Denoted by C^=C/(1-Ch), if h is sufficiently small such that Ch≤1/2, we have
(3.23)∥D2UNj+1∥2≤(1+C^h)∥D2UNj∥2+C^2h≤(1+C^h)j+1∥D2UN0∥2+C^2h∑i=1j‍(1+C^h)i≤exp(j+1)C^h(∥D2UN0∥2+C^jh)≤expC^T(∥D2UN0∥2+C^T)≡C',
where C'=C'(u0,m,T,γ1,γ2)>0 is a positive constant. By the embedding theorem, the estimate (3.15) holds.Next, we investigate the error estimates for the numerical solutionUNj to the exact solution u(tj). Our analysis is based on the error decomposition denoted by
(3.24)u(tj)-UNj=ηj+ej,ηj=u(tj)-PNu(tj),ej=PNu(tj)-UNj.
The boundedness estimate of ηj follows from the inequality (2.6), that is, for any 0≤j≤Λ, there is a positive constant C=C(u0,m,γ1,γ2) such that
(3.25)∥ηj∥≤CN-4.
Hence, it remains to obtain the approximate bounds of ej. If no confusion occurs, we denote the average of the two instant errors ej and ej+1 by e-j+1/2:
(3.26)e-j+1/2=ej+ej+12,η-j+1/2=ηj+ηj+12.
For later use, we give some estimates in the next lemmas.Lemma 3.3.
Assume that the solution of (1.2)–(1.4) is such that uttt∈L2(QT), then
(3.27)∥ej+1∥2≤∥ej∥2+2h(ut(tj+1/2)-UNj+1-UNjh,e-j+1/2)+1320h4∫tjtj+1‍∥uttt∥2dt+h∥e-j+1/2∥2.Proof.
Applying Taylor expansion abouttj+1/2, we have
(3.28)uj=uj+1/2-h2utj+1/2+h28uttj+1/2-12∫tjtj+1/2‍(t-tj)2utttdt,uj+1=uj+1/2+h2utj+1/2+h28uttj+1/2+12∫tj+1/2tj+1‍(tj+1-t)2utttdt.
Then
(3.29)utj+1/2-uj+1-ujh=-12h(∫tj+1/2tj+1‍(tj+1-t)2utttdt+∫tjtj+1/2‍(t-tj)2utttdt).
From Hölder inequality it follows that
(3.30)∥utj+1/2-uj+1-ujh∥2≤12h2(∥∫tj+1/2tj+1‍(tj+1-t)2utttdt∥2+∥∫tjtj+1/2‍(t-tj)2utttdt∥2)≤h3320∫tjtj+1‍∥uttt∥2dt.
Noticing that for any v∈SN, we have
(3.31)(utj+1/2-UNj+1-UNjh,v)=(utj+1/2-uj+1-ujh,v)+1h(ej+1-ej,v).
Taking v=2e-j+1/2 in (3.31), we obtain
(3.32)∥ej+1∥2=∥ej∥2+2h(utj+1/2-UNj+1-UNjh,e-j+1/2)-2h(utj+1/2-uj+1-ujΔt,e-j+1/2)≤∥ej∥2+2h(utj+1/2-UNj+1-UNjh,e-j+1/2)+1320h4∫tjtj+1‍∥uttt∥2dt+h∥e-j+1/2∥2.Takingv=e-j+1/2 in (2.2) and (3.1), respectively, we have
(3.33)(∂u(tj+1/2)∂t,e-j+1/2)+(D2uj+1/2-A(uj+1/2),D(m(uj+1/2)De-j+1/2))=0,(3.34)(UNj+1-UNjh,e-j+1/2)+(D2U-Nj+1/2-PNA~(UNj+1,UNj),D(m(U-Nj+1/2)De-j+1/2))=0.
Comparing (3.33) and (3.34), we have
(3.35)(utj+1/2-UNj+1-UNjh,e-j+1/2)=-(D2uj+1/2-D2UNj+1+D2UNj2,D(m(U-Nj+1/2)De-j+1/2))+(A(uj+1/2)-PNA~(UNj+1,UNj),D(m(U-Nj+1/2)De-j+1/2))+(D2uj+1/2-A(uj+1/2),D(m(U-Nj+1/2)-m(uj+1/2))De-j+1/2).
Now we investigate the error estimates of the three items in the right-hand side of the previous equation.Lemma 3.4.
Assume thatu is the solution of (1.2)–(1.4) such that Dutt∈L2(QT), then there exists a positive constant C1=C1(m,u0,T,γ1,γ2)>0 such that
(3.36)-(D2uj+1/2-D2UNj+1+D2UNj2,D(m(U-Nj+1/2)De-j+1/2))≤-m02∥D2e-j+1/2∥2+C1(N-4+∥ej∥2+∥ej+1∥2+h3∫tjtj+1‍∥D2utt∥2dt).Proof.
By Taylor expansion and Hölder inequality, we obtain(3.37)uj=uj+1/2-h2utj+1/2+∫tjtj+1/2‍(t-tj)uttdt,uj+1=uj+1/2+h2utj+1/2+∫tj+1/2tj+1‍(tj+1-t)uttdt.
Therefore,
(3.38)12(uj+uj+1)-uj+1/2=12(∫tjtj+1/2‍(t-tj)uttdt+∫tj+1/2tj+1‍(tj+1-t)uttdt).
By Hölder inequality, we have
(3.39)∥D2(uj+1/2-12(uj+uj+1))∥2=14(∥D2(∫tjtj+1/2‍(t-tj)uttdt+∫tj+1/2tj+1‍(tj+1-t)uttdt)∥2)≤14(∥(∫tjtj+1/2‍(t-tj)2dt∫tjtj+1/2‍(D2utt)2dt+∫tj+1/2tj+1‍(tj+1-t)2dt∫tj+1/2tj+1‍(D2utt)2dt)∥2)=14(∥(h324∫tjtj+1/2‍(D2utt)2dt+h324∫tj+1/2tj+1‍(D2utt)2dt)∥2)≤14·h324∫tjtj+1‍∥D2utt∥2dt.
Direct computation gives
(3.40)∥D2(uj+1/2-12(uj+uj+1))∥2≤h396∫tjtj+1‍∥D2utt∥2dt.
Then
(3.41)-(D2(uj+1/2-UNj+1+UNj2),D(m(U-Nj+1/2)De-j+1/2))=-(D2(uj+1/2-uj+1+uj2),D(m(U-Nj+1/2)De-j+1/2))-(D2(uj+1+uj2-UNj+1+UNj2),D(m(U-Nj+1/2)De-j+1/2))≤∥D2(uj+1/2-uj+1+uj2)∥·∥D(m(U-Nj+1/2)De-j+1/2)∥-(D2(η-j+1/2+e-j+1/2),D(m(U-Nj+1/2)De-j+1/2))≤{h396∫tjtj+1‍∥D2utt∥2dt}1/2·(∥m(U-Nj+1/2)D2e-j+1/2∥+∥m′(U-Nj+1/2)DU-Nj+1/2De-j+1/2∥)-(D2η-j+1/2,m(U-Nj+1/2)D2e-j+1/2)-(D2η-j+1/2,m′(U-Nj+1/2)DU-Nj+1/2De-j+1/2)-(D2e-j+1/2,m(U-Nj+1/2)D2e-j+1/2)-(D2e-j+1/2,m′(U-Nj+1/2)DU-Nj+1/2De-j+1/2)≜σ1j+σ2j+σ3j.
By Cauchy inequality, it follows that
(3.42)|σ1j|≤{h396∫tjtj+1‍∥D2utt∥2dt}1/2·∥m(U-Nj+1/2)D2e-j+1/2∥+{h396∫tjtj+1‍∥D2utt∥2dt}1/2∥m′(U-Nj+1/2)DU-Nj+1/2De-j+1/2∥≤M0{h396∫tjtj+1‍∥D2utt∥2dt}1/2∥D2e-j+1/2∥+M1∥DU-Nj+1/2∥∞{h396∫tjtj+1‍∥D2utt∥2dt}1/2∥De-j+1/2∥≤h3C1ε∫tjtj+1‍∥D2utt∥2dt+ε∥D2e-j+1/2∥2+ε∥De-j+1/2∥2≤h3C1ε∫tjtj+1‍∥D2utt∥2dt+ε∥D2e-j+1/2∥2+ε∥D2e-j+1/2∥∥e-j+1/2∥≤h3C1ε∫tjtj+1‍∥D2utt∥2dt+2ε∥D2e-j+1/2∥2+ε∥e-j+1/2∥2,|σ2j|≤M0∥D2ηj+1/2∥∥D2e-j+1/2∥+M1∥DUNj+1/2∥∞∥D2ηj+1/2∥∥De-j+1/2∥≤M022ε∥D2ηj+1/2∥2+ε2∥D2e-j+1/2∥2+M12∥DU-Nj+1/2∥∞22εε2∥De-j+1/2∥2≤C2ε∥D2ηj+1/2∥2+ε2∥D2e-j+1/2∥2+ε2∥D2e-j+1/2∥∥e-j+1/2∥≤C2ε∥D2ηj+1/2∥2+ε∥D2e-j+1/2∥2+ε2∥e-j+1/2∥2,σ3j≤-m0∥D2e-j+1/2∥2+M1∥DU-Nj+1/2∥∞∥D2e-j+1/2∥·∥De-j+1/2∥≤-m0∥D2e-j+1/2∥2+ε2∥D2e-j+1/2∥2+M12∥DU-Nj+1/2∥∞22ε∥De-j+1/2∥2≤-m0∥D2e-j+1/2∥2+ε∥D2e-j+1/2∥2+C3ε∥e-j+1/2∥2.
Then we obtain
(3.43)-(D2uj+1/2-D2UNj+1+D2UNj2,D(m(U-Nj+1/2)De-j+1/2))≤|σ1j|+|σ2j|+σ3j≤(4ε-m0)∥D2e-j+1/2∥2+h3C1ε∫tjtj+1‍∥D2utt∥2dt+C2ε∥D2ηj+1/2∥2+C3ε∥e-j+1/2∥2,
where C1ε,C2ε, and C3ε are positive constants. Choosing ε=m0/8, and terms of the properties of the projection operator PN, we complete the proof of the estimate (3.36).Lemma 3.5.
Assume thatu is the solution of (1.2)–(1.4) such that utt∈L2(QT) and ut∈L∞(QT). Then for any positive constant ε>0, there exists a constant C2=C2(m,u0,T,γ1,γ2)>0, such that
(3.44)(A(uj+1/2)-PNA~(UNj+1,UNj),D(m(U-Nj+1/2)De-j+1/2))≤m04∥D2e-j+1/2∥2+C2{∫tjtj+1‍∥ej+1∥2+∥ej∥2+N-8≤m04∥D2e-j+1/2∥2+C2C2+h3(∫tjtj+1‍∥utt∥2dt+∥utj+1/2∥∞2∫tjtj+1‍∥ut∥2dt)}.Proof.
Firstly, we consider(3.45)A(uj+1/2)-A~(uj+1,uj)=γ2(uj+1/2)3-γ24[(uj+1)3+(uj+1)2uj+uj+1(uj)2+(uj)3]+γ1(uj+1/2)2-γ13[(uj+1)2+uj+1uj+(uj)2]-(uj+1/2-12(uj+1+uj))≜γ2ρ1j+γ1ρ2j-ρ3j.
Direct computation gives
(3.46)∥ρ3j∥=12∥∫tjtj+1/2‍(t-tj)uttdt+∫tj+1/2tj+1‍(tj+1-t)uttdt∥≤(Δt396∫tjtj+1‍∥utt∥2dt)1/2,∥ρ2j∥=∥16[(2uj+1/2)2-(uj+1+uj)2]+16[2(uj+1/2)2-(uj+1)2-(uj)2]∥≤16∥[2uj+1/2-(uj+1+uj)][2uj+1/2+(uj+1+uj)]∥+16∥(uj+1/2-uj+1)(uj+1/2+uj+1)+(uj+1/2-uj)(uj+1/2+uj)∥≤16∥2uj+1/2-(uj+1+uj)∥·∥2uj+1/2+(uj+1+uj)∥∞+16∥(-Δt2utj+1/2-∫tj+1/2tj+1‍(tj+1-t)uttdt)(uj+1/2+uj+1)+16+16+(Δt2utj+1/2-∫tjtj+1/2‍(t-tj)uttdt)(uj+1/2+uj)∥≤CΔt3/2(∫tjtj+1‍∥utt∥2dt)1/2+Δt3/212∥utj+1/2∥∞·(∫tjtj+1‍∥ut∥2dt)1/2,∥ρ1j∥=112∥((2uj+1/2)3-(uj+1+uj)3)+16[2(uj+1/2)3-(uj+1)3-(uj)3]∥≤112∥[2uj+1/2-(uj+1+uj)]·[4(uj+1/2)2+2(uj+1+uj)uj+1/2+(uj+1+uj)2]∥+16∥(uj+1/2-uj+1)((uj+1/2)2+uj+1uj+1/2+(uj+1)2)+16+16+(uj+1/2-uj)((uj+1/2)2+ujuj+1/2+(uj)2)∥≤CΔt3/2{(∫tjtj+1‍∥utt∥2dt)1/2+∥utj+1/2∥∞·(∫tjtj+1‍∥ut∥2dt)1/2}.
Then
(3.47)(A(uj+1/2)-PNA~(UNj+1,UNj),D(m(U-Nj+1/2)De-j+1/2))≤34ε∥A(uj+1/2)-PNA(uj+1/2)∥2+ε3∥D(m(U-Nj+1/2)De-j+1/2)∥2+34ε∥A(uj+1/2)-A~(uj+1,uj)∥2+ε3∥D(m(U-Nj+1/2)De-j+1/2)∥2+34ε∥A~(uj+1,uj)-A~(UNj+1,UNj)∥2+ε3∥D(m(U-Nj+1/2)De-j+1/2)∥2≤ε∥D(m(U-Nj+1/2)De-j+1/2)∥2+34ε∥A(uj+1/2)-PNA(uj+1/2)∥2+34ε∥A(uj+1/2)-A~(uj+1,uj)∥2+34ε∥A~(uj+1,uj)-A~(UNj+1,UNj)∥2.
Direct computation yields
(3.48)∥D(m(U-Nj+1/2)De-j+1/2)∥2≤(M02+M12)∥D2e-j+1/2∥2+C∥e-j+1/2∥2,∥A(uj+1/2)-PNA(uj+1/2)∥2≤CN-8,∥A(uj+1/2)-A~(uj+1,uj)∥2≤C(∥ρ1j∥2+∥ρ2j∥2+∥ρ3j∥2),∥A~(uj+1,uj)-A~(UNj+1,UNj)∥2≤(∥G1j∥∞2+∥G2j∥∞2)(∥ej+1∥2+∥ej∥2+CN-8),
where
(3.49)G1j=γ28((uj+1+UNj+1)2+(UNj+1+UNj)2+(uj+1+UNj)2)+γ13(UNj+1+uj+1+UNj)-12,G2j=γ28((uj+UNj)2+(uj+uj+1)2+(uj+1+UNj)2)+γ13(uj+1+uj+UNj)-12.
Applying Lemma 2.3 and Theorem 3.7, we obtain that ∥G1∥∞≤C(u0,m,T,γ1,γ2) and ∥G2∥∞≤C(u0,m,T,γ1,γ2). Taking ε=m0/4(M02+M12) in (3.47),we have
(3.50)(A(uj+1/2)-PNA~(UNj+1,UNj),D(mj+1/2De-j+1/2))≤m04∥D2e-j+1/2∥2+C2{∫tjtj+1‍∥ej+1∥2+∥ej∥2+N-8≤m04∥D2e-j+1/2∥2+C2m04C2C2+Δt3(∫tjtj+1‍∥utt∥2dt+∥utj+1/2∥∞2∫tjtj+1‍∥ut∥2dt)},
where C2=C2(u0,m,T,γ1,γ2)>0 is a constant.Lemma 3.6.
Assume thatu is the solution of (1.2)–(1.4). Then there exists a positive constant C=C(u0,m,γ1,γ2)>0 such that
(3.51)(D2uj+1/2-A(uj+1/2),D(m(U-Nj+1/2)-m(uj+1/2))De-j+1/2)≤m04∥D2e-j+1/2∥2+C(∥e-j+1/2∥2+N-8+h3∫tjtj+1‍∥utt∥dt).Proof.
By (2.9), we have
(3.52)∥D3uj+1/2-DA(uj+1/2)∥∞≤∥D3uj+1/2∥∞+∥A′(uj+1/2)Duj+1/2∥∞≤∥D3uj+1/2∥∞+∥A′(uj+1/2)∥∞∥Duj+1/2∥∞≤C.
In the other hand,
(3.53)∥m(U-Nj+1/2)-m(uj+1/2)∥2≤M12∥(U-Nj+1/2-uj+1/2)∥2≤C∥U-Nj+1/2-uj+uj+12∥2+∥uj+uj+12-uj+1/2∥2≤C∥e-j+1/2∥2+∥η-j+1/2∥2+h396∫tjtj+1‍∥utt∥dt.
By Young inequality, we obtain
(3.54)(D2uj+1/2-A(uj+1/2),D(m(U-Nj+1/2)-m(uj+1/2))De-j+1/2)=-(D3uj+1/2-DA(uj+1/2),m(U-Nj+1/2)-m(uj+1/2)De-j+1/2)≤∥D3uj+1/2-DA(uj+1/2)∥∞∥m(U-Nj+1/2)-m(uj+1/2)∥∥De-j+1/2∥≤ε∥De-j+1/2∥2+Cε∥m(U-Nj+1/2)-m(uj+1/2)∥2≤ε∥D2e-j+1/2∥∥e-j+1/2∥+Cε∥m(U-Nj+1/2)-m(uj+1/2)∥2≤ε∥D2e-j+1/2∥2+Cε(∥e-j+1/2∥2+N-8+h3∫tjtj+1‍∥utt∥dt).
Choosing ε=m0/4 in the previous inequality leads to (3.51).Finally, we obtain the main theorem of this paper.Theorem 3.7.
Assume thatu(x,t) is the solution of (1.2)–(1.4) and satisfies that
(3.55)D4u∈L∞(QT),ut∈L∞(QT),D2utt∈L2(QT),uttt∈L2(QT).UNj∈SN(j=1,2,…,Λ) is the solution of the full-discretization (3.1) and (3.2). If h is sufficiently small, there exists a positive constant C such that, for any j=0,1,2,…,Λ,
(3.56)∥ej+1∥=∥PNu(tj+1)-UNj+1∥≤C(N-2+∥e0∥+Bh2),
where B=∫0tj+1‍(∥D2utt∥2+∥utt∥2+∥uttt∥2+max0≤l≤j{∥utl+1/2∥∞2}·∥ut∥2)dt.Proof.
By (3.27), (3.36), (3.44), and (3.51), we have (3.57)∥ej+1∥2≤∥ej∥2+h4320∫tjtj+1‍∥uttt∥2dt+h∥e-j+1/2∥2+2h{-(D2uj+1/2-D2UNj+1+D2UNj2,D(m(U-Nj+1/2)De-j+1/2))+(A(uj+1/2)-A~(UNj+1,UNj),D(m(U-Nj+1/2)De-j+1/2))+(D2uj+1/2-A(uj+1/2),D(m(U-Nj+1/2)-m(uj+1/2))De-j+1/2)D2UNj+1+D2UNj2}≤∥ej∥2+hC1(N-4+∥ej+1∥2+∥ej∥2)+C2h4∫tjtj+1‍(∥D2utt∥2+∥utt∥2+∥uttt∥2+∥utj+1/2∥∞2∥ut∥2)dt,
where C1=C1(m,u0,T,γ1,γ2)>0 and C2=C2(m,u0,T,γ1,γ2)>0 are constants. For sufficiently small h such that C1h≤1/2, denoting C~=2(C1+C2), we obtain
(3.58)∥ej+1∥2≤(1+C~h)∥ej∥2+C~(hN-4+h4Bj),
where
(3.59)Bj=∫tjtj+1‍(∥D2utt∥2+∥utt∥2+∥uttt∥2+∥utj+1/2∥∞2∥ut∥2)dt.
By the Gronwall inequality of the discrete form, we obtain
(3.60)∥ej+1∥2≤(1+C~h)j+1∥e0∥2+C~∑l=0j‍(1+C~h)l(htN-4+h4Bl)≤exp{C~(j+1)h}{∥e0∥2+C~∑l=0j‍(hN-4+h4Bl)}≤exp{C~(j+1)h}{∥e0∥2+C~(jhN-4+h4∑l=0j‍Bl)}.
Direct computation gives
(3.61)∑l=0j‍Bl≤∫0tj+1‍(∥D2utt∥2+∥utt∥2+∥uttt∥2+max0≤l≤j{∥utl+1/2∥∞2}·∥ut∥2)dt.
Then we complete the conclusion (3.56).Furthermore, we get the following theorem.Theorem 3.8.
Assumeu is the solution of (1.2)–(1.4) and satisfies that
(3.62)D4u∈L∞(QT),ut∈L∞(QT),D2utt∈L2(QT),uttt∈L2(QT).UNj∈SN(j=1,2,…,Λ) is the solution of the full-discretization scheme (3.1)-(3.2), and U0 satisfies ∥e0∥=∥PNu0-U0∥≤CN-4. If h is sufficiently small, there exists a constant C=C(u0,m,T,γ1,γ2)>0 such that
(3.63)∥u(x,tj)-UNj∥≤C(N-2+h2),j=1,2,…,Λ.


## 4. Numerical Experiments

In this section, we apply the spectral method described in (3.1) and (3.2) to carry out numerical computations to illustrate theoretical estimations in the previous section. Consider (2.2) with settings:
(4.1)m(s)=m0+s2,A(s)=s3-s,
where m0>0 is a constant. The full-discretization spectral method of (2.2) and (2.3) reads: find UNj=∑l=0Nαjlcoslπx(j=1,2,…,Λ) such that
(4.2)(UNj+1-UNjh,v)+(D2UNj+1+D2UNj2-PNA~(UNj+1,UNj),D(m(D2UNj+1+D2UNj2)Dv))=0,(UN0,v)-(u0,v)=0.
In our computations we fix N=32 and choose five different time-step sizes hk(k=1,2,…,5). Let Λk be the integer with hkΛk=T. Since we have no exact solution of (2.2) and (2.3), we take N=32 and h0=0.15625×10-4 to compute an approximating solution UNΛ0 with h0Λ0=T and regard this as an exact solution. we also choose five different time-step sizes hk(k=1,2,…,5) with hkΛk=T to obtain five approximating solutions UNΛk(k=1,2,…,5) and compute the error estimation. Define an error function:
(4.3)err(T,hk)=(∫01‍(UNΛk-UNΛ0)2dx)1/2.
This function characterizes the estimations with respect to time-step size.


### 4.1. Example 1

Takem0=1 and T=0.1. We also take two different initial functions u01(x)=x5(1-x)5 and u02=x5(1-x)5ex to carry out numerical computations. Figure 1 shows the development of the solutions for time t from t=0 to t=0.1 with fixed step-size h0=0.15625×10-4.The development of the solution of the full-discrete scheme when initial value isu0(x)=x5(1-x)5 (a) and u0(x)=x5(1-x)5ex (b).
(a)
(b)We also choose five different time-step sizeshk to carry out numerical computations and apply the error function in (4.3) to illustrate the estimation and convergence order in time variable t, see Table 1.Table 1
The error and the convergence order.
h
k
err
1
(
0.1
,
h
k
)
Order
1
err
2
(
0.1
,
h
k
)
Order
2
0
.
1
×
1
0
-
2
0
.
2441
×
1
0
-
6
—
0
.
2114
×
1
0
-
6
—
0
.
5
×
1
0
-
3
0
.
6332
×
1
0
-
7
1.9466
0
.
5414
×
1
0
-
7
1.9651
0
.
25
×
1
0
-
3
0
.
1186
×
1
0
-
7
2.4159
0
.
1115
×
1
0
-
7
2.2799
0
.
125
×
1
0
-
3
0
.
1938
×
1
0
-
8
2.6142
0
.
1901
×
1
0
-
8
2.5515
0
.
625
×
1
0
-
4
0
.
2504
×
1
0
-
9
2.9518
0
.
2483
×
1
0
-
9
2.9368


### 4.2. Example 2

Takem0=0.05 and T=0.1. We also take two different initial functions u01(x)=x5(1-x)5 and u02=x5(1-x)5ex to carry out numerical computations. Figure 2 shows the development of the solutions for time t from t=0 to t=0.1 with fixed step-sizes h0=0.15625×10-4.The development of the solution of the full-discrete scheme with initial valueu0(x)=x5(1-x)5 (a) and u0(x)=x5(1-x)5ex (b).
(a)
(b)We also choose five different time-step sizeshk to carry out numerical computations and apply the error function in (4.3) to illustrate the estimation and convergence order in time variable t, see Table 2.Table 2
The error and the convergence order.
h
k
err
1
(
0.1
,
h
k
)
Order
1
err
2
(
0.1
,
h
k
)
Order
2
0
.
1
×
1
0
-
2
0
.
4748
×
1
0
-
8
—
0
.
8864
×
1
0
-
8
—
0
.
5
×
1
0
-
3
0
.
8348
×
1
0
-
9
2.5077
0
.
1577
×
1
0
-
8
2.4907
0
.
25
×
1
0
-
3
0
.
1370
×
1
0
-
9
2.6069
0
.
2440
×
1
0
-
9
2.6922
0
.
125
×
1
0
-
3
0
.
2694
×
1
0
-
10
2.3467
0
.
4457
×
1
0
-
10
2.4525
0
.
625
×
1
0
-
4
0
.
6398
×
1
0
-
11
2.0741
0
.
1059
×
1
0
-
10
2.0736


### 4.3. Example 3

Takem0=0.005 and T=0.1. We also take two different initial functions u01(x)=x5(1-x)5 and u02=x5(1-x)5ex to carry out numerical computations. Figure 3 shows the development of the solutions for time t from t=0 to t=0.1 with fixed step-sizes h0=0.15625×10-4.The development of the solution of the full-discrete scheme with initial valueu0(x)=x5(1-x)5 (a) and u0(x)=x5(1-x)5ex (b).
(a)
(b)We also choose five different time-step sizeshk to carry out numerical computations and apply the error function in (4.3) to illustrate the estimation and convergence order in time variable t, see Table 3.Table 3
The error and the convergence order.
h
k
err
1
(
0.1
,
h
k
)
Order1
err
2
(
0.1
,
h
k
)
Order2
0.1
×
1
0
-
2
0.1150
×
1
0
-
5
—
0.5224
×
1
0
-
5
—
0.5
×
1
0
-
3
0.2872
×
1
0
-
6
2.0013
0.1304
×
1
0
-
5
2.0019
0.25
×
1
0
-
3
0.7158
×
1
0
-
7
2.0043
0.3250
×
1
0
-
6
2.0044
0.125
×
1
0
-
3
0.1769
×
1
0
-
7
2.0170
0.8031
×
1
0
-
7
2.0171
0.625
×
1
0
-
4
0.4211
×
1
0
-
8
2.0704
0.1912
×
1
0
-
7
2.0704


## 4.1. Example 1

Takem0=1 and T=0.1. We also take two different initial functions u01(x)=x5(1-x)5 and u02=x5(1-x)5ex to carry out numerical computations. Figure 1 shows the development of the solutions for time t from t=0 to t=0.1 with fixed step-size h0=0.15625×10-4.The development of the solution of the full-discrete scheme when initial value isu0(x)=x5(1-x)5 (a) and u0(x)=x5(1-x)5ex (b).
(a)
(b)We also choose five different time-step sizeshk to carry out numerical computations and apply the error function in (4.3) to illustrate the estimation and convergence order in time variable t, see Table 1.Table 1
The error and the convergence order.
h
k
err
1
(
0.1
,
h
k
)
Order
1
err
2
(
0.1
,
h
k
)
Order
2
0
.
1
×
1
0
-
2
0
.
2441
×
1
0
-
6
—
0
.
2114
×
1
0
-
6
—
0
.
5
×
1
0
-
3
0
.
6332
×
1
0
-
7
1.9466
0
.
5414
×
1
0
-
7
1.9651
0
.
25
×
1
0
-
3
0
.
1186
×
1
0
-
7
2.4159
0
.
1115
×
1
0
-
7
2.2799
0
.
125
×
1
0
-
3
0
.
1938
×
1
0
-
8
2.6142
0
.
1901
×
1
0
-
8
2.5515
0
.
625
×
1
0
-
4
0
.
2504
×
1
0
-
9
2.9518
0
.
2483
×
1
0
-
9
2.9368


## 4.2. Example 2

Takem0=0.05 and T=0.1. We also take two different initial functions u01(x)=x5(1-x)5 and u02=x5(1-x)5ex to carry out numerical computations. Figure 2 shows the development of the solutions for time t from t=0 to t=0.1 with fixed step-sizes h0=0.15625×10-4.The development of the solution of the full-discrete scheme with initial valueu0(x)=x5(1-x)5 (a) and u0(x)=x5(1-x)5ex (b).
(a)
(b)We also choose five different time-step sizeshk to carry out numerical computations and apply the error function in (4.3) to illustrate the estimation and convergence order in time variable t, see Table 2.Table 2
The error and the convergence order.
h
k
err
1
(
0.1
,
h
k
)
Order
1
err
2
(
0.1
,
h
k
)
Order
2
0
.
1
×
1
0
-
2
0
.
4748
×
1
0
-
8
—
0
.
8864
×
1
0
-
8
—
0
.
5
×
1
0
-
3
0
.
8348
×
1
0
-
9
2.5077
0
.
1577
×
1
0
-
8
2.4907
0
.
25
×
1
0
-
3
0
.
1370
×
1
0
-
9
2.6069
0
.
2440
×
1
0
-
9
2.6922
0
.
125
×
1
0
-
3
0
.
2694
×
1
0
-
10
2.3467
0
.
4457
×
1
0
-
10
2.4525
0
.
625
×
1
0
-
4
0
.
6398
×
1
0
-
11
2.0741
0
.
1059
×
1
0
-
10
2.0736


## 4.3. Example 3

Takem0=0.005 and T=0.1. We also take two different initial functions u01(x)=x5(1-x)5 and u02=x5(1-x)5ex to carry out numerical computations. Figure 3 shows the development of the solutions for time t from t=0 to t=0.1 with fixed step-sizes h0=0.15625×10-4.The development of the solution of the full-discrete scheme with initial valueu0(x)=x5(1-x)5 (a) and u0(x)=x5(1-x)5ex (b).
(a)
(b)We also choose five different time-step sizeshk to carry out numerical computations and apply the error function in (4.3) to illustrate the estimation and convergence order in time variable t, see Table 3.Table 3
The error and the convergence order.
h
k
err
1
(
0.1
,
h
k
)
Order1
err
2
(
0.1
,
h
k
)
Order2
0.1
×
1
0
-
2
0.1150
×
1
0
-
5
—
0.5224
×
1
0
-
5
—
0.5
×
1
0
-
3
0.2872
×
1
0
-
6
2.0013
0.1304
×
1
0
-
5
2.0019
0.25
×
1
0
-
3
0.7158
×
1
0
-
7
2.0043
0.3250
×
1
0
-
6
2.0044
0.125
×
1
0
-
3
0.1769
×
1
0
-
7
2.0170
0.8031
×
1
0
-
7
2.0171
0.625
×
1
0
-
4
0.4211
×
1
0
-
8
2.0704
0.1912
×
1
0
-
7
2.0704



---

*Source: 808216-2012-10-02.xml*