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

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.

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: 𝜕𝑢[(𝑢𝜕𝑡+div𝑘∇Δ𝑢−∇𝐴)]=0. (1.1) 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 is 𝐶^𝑘 smooth, the approximate solution will be convergent to the exact solution with power for exponent 𝑁^−𝑘, where 𝑁 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: 𝜕𝑢𝑚(𝑢)𝐷𝜕𝑡+𝐷³(𝑢)()∈()×(),|𝑢−𝐷𝐴=0,𝑥,𝑡0,10,𝑇𝐷𝑢_𝑥=0,1=𝐷³𝑢_𝑥=0,1(),𝑢()=0,𝑡∈0,𝑇𝑥,0=𝑢₀,𝑥∈(0,1), (1.2)(1.3)(1.4) where 𝐷=𝜕/𝜕𝑥 and 𝐴(𝑠)=−𝑠+𝛾₁𝑠²+𝛾₂𝑠³,𝛾₂>0. (1.5) Here, 𝑢(𝑥,𝑡) represents a relative concentration of one component in binary mixture. The function 𝑚(𝑢) is the mobility which depends on the unknown function 𝑢, which restricts diffusion of both components to the interfacial region only. Denote 𝑄_𝑇=(0,1)×(0,𝑇). Throughout this paper, we assume that 0<𝑚₀≤𝑚(𝑠)≤𝑀₀,𝑚^′(𝑠)≤𝑀₁,∀𝑠∈ℝ, (1.6) where 𝑚₀, 𝑀₀, and 𝑀₁ 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 Section 2, we study the spectral method for (1.2)–(1.4) and obtain the error estimate between the exact solution 𝑢 and the spectral approximate solution 𝑢_𝑁. In Section 3, we use the implicit Euler method to discretize the time variable and obtain the error estimate between the exact solution 𝑢 and the full-discrete approximate solution 𝑈^𝑗_𝑁. 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 ‖⋅‖_𝑘 and |⋅|_𝑘 the norm and seminorm of the Sobolev spaces 𝐻^𝑘(0,1)(𝑘∈ℕ), respectively. Let (⋅,⋅) be the standard 𝐿² inner product over (0,1). Define 𝐻²_𝐸(0,1)=𝑣∈𝐻²(0,1);𝐷𝑣|_𝑥=0,1=0. (2.1)

A function 𝑢 is said to be a weak solution of the problems (1.2)–(1.4), if 𝑢∈𝐿^∞(0,𝑇;𝐻²_𝐸(0,1)) and satisfies the following equations: 𝜕𝑢+𝐷𝜕𝑡,𝑣²𝑢−𝐴(𝑢),𝐷(𝑚(𝑢)𝐷𝑣)=0,∀𝑣∈𝐻²_𝐸(𝑢())−𝑢,⋅,0,𝑣₀,𝑣=0,∀𝑣∈𝐻²_𝐸. (2.2)(2.3)

For any integer 𝑁>0, let 𝑆_𝑁=span{cos𝑘𝜋𝑥,𝑘=0,1,2,…,𝑁}. Define a projection operator 𝑃_𝑁∶𝐻²_𝐸↦𝑆_𝑁 by ¹₀𝑃_𝑁𝑢(𝑥)𝑣(𝑥)𝑑𝑥=¹₀𝑢(𝑥)𝑣(𝑥)𝑑𝑥,∀𝑣∈𝑆_𝑁. (2.4) We collect some properties of this projection 𝑃_𝑁 in the following lemma (see [39]).

Lemma 2.1. (i)  𝑃_𝑁 commutes with the second derivation on 𝐻²_𝐸(𝐼), that is, 𝑃_𝑁𝐷²𝑣=𝐷²𝑃_𝑁𝑣,∀𝑣∈𝐻²_𝐸(𝐼). (2.5)
(ii) For any 0≤𝜇≤𝜎, there exists a positive constant 𝐶 such that 𝑣−𝑃_𝑁𝑣_𝜇≤𝐶𝑁^𝜇−𝜎|𝑣|_𝜎,𝑣∈𝐻^𝜎(𝐼). (2.6)

The following Nirenberg inequality is a key tool for our theoretical estimates.

Lemma 2.2. Assume that Ω⊂ℝ^𝑛 is a bounded domain, 𝑢∈𝑊^𝑚,𝑟(Ω), then we have 𝐷^𝑗𝑢_𝐿𝑝⩽𝐶₁𝐷^𝑚𝑢^𝑎_𝐿𝑟‖𝑢‖_𝐿^1−𝑎𝑞+𝐶₂‖𝑢‖_𝐿^𝑞, (2.7) where 𝑗𝑚1⩽𝑎<1,𝑝=𝑗𝑛1+𝑎𝑟−𝑚𝑛+()11−𝑎𝑞. (2.8)

By [41], we have the following.

Lemma 2.3. Assume that 𝑚(𝑠)∈𝐶^1+𝛼(ℛ), 𝑢₀∈𝐶^4+𝛼(𝐼), 𝐷^𝑖𝑢₀(0)=𝐷^𝑖𝑢₀(1)=0(𝑖=0,1,2,3,4),𝑚(𝑠)>0, then there exists a unique solution of the problems (1.2)–(1.4) such that 𝑢∈𝐶^{1+𝛼/4,4+𝛼}𝑄_𝑇. (2.9)

The spectral approximation to (2.2) is to find an element 𝑢_𝑁(⋅,𝑡)∈𝑆_𝑁 such that 𝜕𝑢_𝑁𝜕𝑡,𝑣_𝑁+𝐷²𝑢_𝑁−𝑃_𝑁𝐴𝑢_𝑁,𝐷𝑚𝑢_𝑁𝐷𝑣_𝑁=0,∀𝑣_𝑁∈𝑆_𝑁,𝑢_𝑁(⋅,0),𝑣_𝑁−𝑢₀,𝑣_𝑁=0,∀𝑣_𝑁∈𝑆_𝑁. (2.10)(2.11)

Now we study the 𝐿^∞ norm estimates of the function 𝑢_𝑁(⋅,𝑡) and 𝐷𝑢_𝑁(⋅,𝑡) for 0≤𝑡≤𝑇.

Theorem 2.4. Assume (1.6) and 𝑢₀∈𝐻²_𝐸. Then there is a unique solution of (2.10) and (2.11) such that 𝑢_𝑁∞≤𝐶,𝐷𝑢_𝑁≤𝐶,0≤𝑡≤𝑇, (2.12) where 𝐶=𝐶(𝑢₀,𝑇,𝛾₁,𝛾₂) is a positive constant.

Proof. From (2.11) it follows that 𝑢_𝑁(⋅,0)=𝑃_𝑁𝑢₀(⋅). The existence and uniqueness of the initial problem follow from the standard ODE theory. Now we study the estimate.
Define an energy function: 𝐹(𝑡)=𝐻𝑢_𝑁+1,12𝐷𝑢_𝑁², (2.13) where 𝐻(𝑠)=∫^𝑠₀𝐴(𝑡)𝑑𝑡=(1/4)𝛾₂𝑠⁴+(1/3)𝛾₁𝑠³−(1/2)𝑠². Direct computation gives 𝑑𝐹=𝐴𝑢𝑑𝑡^𝑁,𝜕𝑢_𝑁+𝜕𝑡𝐷𝑢_𝑁,𝐷𝜕𝑢_𝑁=𝑃𝜕𝑡^𝑁𝐴𝑢_𝑁,𝜕𝑢_𝑁−𝐷𝜕𝑡²𝑢_𝑁,𝜕𝑢_𝑁=𝜕𝑡𝜕𝑢_𝑁𝜕𝑡,𝑃_𝑁𝐴𝑢_𝑁−𝐷²𝑢_𝑁. (2.14) Noticing that 𝑃_𝑁𝐴(𝑢_𝑁)−𝐷²𝑢_𝑁∈𝑆_𝑁 and setting 𝑣_𝑁=𝑃_𝑁𝐴(𝑢_𝑁)−𝐷²𝑢_𝑁 in (2.10), applying integration by part, we obtain 𝑑𝐹=𝑚𝑢𝑑𝑡^𝑁𝐷³𝑢_𝑁−𝐷𝑃_𝑁𝐴𝑢_𝑁,𝐷𝑃_𝑁𝐴𝑢_𝑁−𝐷²𝑢_𝑁=−𝑚𝑢_𝑁𝐷³𝑢_𝑁−𝐷𝑃_𝑁𝐴𝑢_𝑁,𝐷³𝑢_𝑁−𝐷𝑃_𝑁𝐴𝑢_𝑁≤−𝑚₀𝐷𝑃_𝑁𝐴𝑢_𝑁−𝐷²𝑢_𝑁²≤0. (2.15) Hence, 𝐹(𝑡)≤𝐹(0)=¹₀𝐻𝑃_𝑁𝑢₀1𝑑𝑠+2𝐷𝑃_𝑁𝑢₀²,∀0≤𝑡≤𝑇. (2.16) Applying Young inequality, we obtain 𝑢²_𝑁≤𝜀𝑢⁴_𝑁+𝐶_1𝜀,𝑢³_𝑁≤𝜀𝑢⁴_𝑁+𝐶_2𝜀, (2.17) where 𝐶_1𝜀 and 𝐶_2𝜀 are positive constants. Letting 𝜀=3𝛾₂/(8|𝛾₁|+12), then for all 0≤𝑡≤𝑇, ¹₀𝐻𝑢_𝑁𝑑𝑥≥¹₀14𝛾₂𝑢⁴_𝑁−13𝛾₁𝑢³_𝑁−12𝑢²_𝑁𝛾𝑑𝑥≥²8¹₀𝑢_𝑁⁴𝑑𝑥−𝐾₁, (2.18) where 𝐾₁ is a positive constant depending only on 𝛾₁ and 𝛾₂. Therefore, we get 12𝐷𝑢_𝑁(⋅,𝑡)²+𝛾₂8¹₀𝑢_𝑁(⋅,𝑡)⁴𝑑𝑥−𝐾₁≤12𝐷𝑢_𝑁(⋅,𝑡)²+¹₀𝐻𝑢_𝑁()(𝑡)(0)=⋅,𝑡𝑑𝑥=𝐹≤𝐹¹₀𝐻𝑃_𝑁𝑢₀𝑑𝑥+𝐷𝑃_𝑁𝑢₀². (2.19) Thus, 𝐷𝑢_𝑁(⋅,𝑡)²≤2𝐾₁+2𝐹(0)=2𝐾₁+2¹₀𝐻𝑃_𝑁𝑢₀𝑑𝑥+𝐷𝑃_𝑁𝑢₀²≜𝐶,₁₀𝑢_𝑁(𝑥,𝑡)⁴8𝑑𝑥≤𝛾₂𝐾₁(0=8+𝐹)𝛾₂𝐾₁+¹₀𝐻𝑃_𝑁𝑢₀1𝑑𝑥+2𝐷𝑃_𝑁𝑢₀²≜𝐶, (2.20) where 𝐶=𝐶(𝑢₀,𝛾₁,𝛾₂) is a constant. By Hölder inequality, we obtain 𝑢_𝑁²=¹₀𝑢²_𝑁𝑑𝑥≤¹₀𝑢⁴_𝑁𝑑𝑥^1/21₀1²𝑑𝑥^1/2=¹₀𝑢⁴_𝑁𝑑𝑥^1/2. (2.21) Therefore, 𝑢_𝑁≤𝐶,𝐷𝑢_𝑁≤𝐶. (2.22) From the embedding theorem it follows that 𝑢_𝑁∞≤𝐶,∀0≤𝑡≤𝑇. (2.23)

Theorem 2.5. Assume (1.6) and let 𝑢_𝑁(⋅,𝑡) be the solution of (2.10) and (2.11). Then there is a positive constant 𝐶=𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂)>0 such that 𝐷𝑢_𝑁∞≤𝐶,𝐷²𝑢_𝑁≤𝐶. (2.24)

Proof. Setting 𝑣_𝑁=𝐷⁴𝑢_𝑁 in (2.10) and integrating by parts, we get 12𝑑𝐷𝑑𝑡²𝑢_𝑁²+𝑚𝑢_𝑁𝐷⁴𝑢_𝑁−𝐷²𝑃_𝑁𝐴𝑢_𝑁,𝐷⁴𝑢_𝑁+𝑚^′𝑢_𝑁𝐷𝑢_𝑁𝐷³𝑢_𝑁−𝐷𝑃_𝑁𝐴𝑢_𝑁,𝐷⁴𝑢_𝑁=0. (2.25) Consequently, 12𝑑𝐷𝑑𝑡²𝑢_𝑁²+𝑚𝑢_𝑁𝐷⁴𝑢_𝑁,𝐷⁴𝑢_𝑁=𝐼₁+𝐼₂+𝐼₃, (2.26) where 𝐼₁=𝑚𝑢_𝑁𝐷²𝑃_𝑁𝐴𝑢_𝑁,𝐷⁴𝑢_𝑁,𝐼₂=−𝑚^′𝑢_𝑁𝐷𝑢_𝑁𝐷³𝑢_𝑁,𝐷⁴𝑢_𝑁,𝐼₃=𝑚^′𝑢_𝑁𝐷𝑢_𝑁𝐷𝑃_𝑁𝐴𝑢_𝑁,𝐷⁴𝑢_𝑁. (2.27) Noticing (1.6), we have 12𝑑𝐷𝑑𝑡²𝑢_𝑁²+𝑚₀𝐷⁴𝑢_𝑁²≤12𝑑𝐷𝑑𝑡²𝑢_𝑁²+𝑚𝑢_𝑁𝐷⁴𝑢_𝑁,𝐷⁴𝑢_𝑁. (2.28) In terms of the Nirenberg inequality (2.7), there is a constant 𝐶>0 such that 𝐷𝑢_𝑁∞≤𝐶𝐷⁴𝑢_𝑁^3/8‖𝑢‖^5/8+‖𝑢‖,𝐷𝑢_𝑁∞≤𝐶𝐷²𝑢_𝑁^3/4‖𝑢‖^1/4+‖𝑢,𝐷₃‖𝑢_𝑁∞≤𝐶𝐷⁴𝑢_𝑁^7/8‖𝑢‖^1/8+‖𝑢‖. (2.29) Noticing the definition of the function 𝐴 and the estimates in (2.12), we have 𝐴𝑢_𝑁∞≤𝐶,𝐴^′𝑢_𝑁∞≤𝐶,𝐴^′′𝑢_𝑁∞≤𝐶, (2.30) for some constant 𝐶=𝐶(𝑢₀,𝑇,𝛾₁,𝛾₂)>0. Applying Hölder inequality and Young inequality, for any 𝜀>0, 𝐼₁≤𝑀₀𝐴^′𝑢_𝑁𝐷²𝑢_𝑁+𝐴^′′𝑢_𝑁𝐷𝑢_𝑁²⋅𝐷⁴𝑢_𝑁≤𝑀₀𝐴^′𝑢_𝑁∞⋅𝐷²𝑢_𝑁𝐷⁴𝑢_𝑁+𝑀₀𝐴^′′𝑢_𝑁∞⋅𝐷𝑢_𝑁𝐷𝑢_𝑁∞𝐷⁴𝑢_𝑁≤𝜀𝐷⁴𝑢_𝑁²+𝐶𝐷²𝑢_𝑁²+𝐶𝐷⁴𝑢_𝑁^3/8+1𝐷⁴𝑢_𝑁≤𝜀𝐷⁴𝑢_𝑁²+𝐶𝐷²𝑢_𝑁²+𝐶𝐷⁴𝑢_𝑁^11/8+𝐶𝐷⁴𝑢_𝑁≤𝜀𝐷⁴𝑢_𝑁²+𝐶𝐷²𝑢_𝑁²+𝜀2𝐷⁴𝑢_𝑁²𝜀+𝐶+2𝐷⁴𝑢_𝑁²𝐷+𝐶≤2𝜀⁴𝑢_𝑁²+𝐶𝐷²𝑢_𝑁²+𝐶, (2.31) where 𝐶=𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂,𝜀)>0 is a positive constant. Similarly, we obtain 𝐼₂≤𝑀₁𝐷𝑢_𝑁𝐷³𝑢_𝑁∞𝐷⁴𝑢_𝑁≤𝐶𝐷⁴𝑢_𝑁^7/8+1𝐷⁴𝑢_𝑁≤𝐶𝐷⁴𝑢_𝑁^15/8+𝐷⁴𝑢_𝑁≤𝜀2𝐷⁴𝑢_𝑁²𝜀+𝐶+2𝐷⁴𝑢_𝑁²𝐷+𝐶≤𝜀⁴𝑢_𝑁²𝐼+𝐶,₃≤𝑀₁𝐷𝑃_𝑁𝐴𝑢_𝑁∞⋅𝐷𝑢_𝑁⋅𝐷⁴𝑢_𝑁≤𝐶𝐷𝑃_𝑁𝐴𝑢_𝑁∞𝐷⁴𝑢_𝑁≤𝐶𝐷²𝑃_𝑁𝐴𝑢_𝑁^3/4𝑃_𝑁𝐴𝑢_𝑁^1/4+𝑃_𝑁𝐴𝑢_𝑁𝐷⁴𝑢_𝑁≤𝐶𝐷²𝑃_𝑁𝐴𝑢_𝑁^3/4𝐴𝑢_𝑁^1/4+𝐴𝑢_𝑁𝐷⁴𝑢_𝑁≤𝐶𝐷²𝑃_𝑁𝐴𝑢_𝑁^3/4𝐴𝑢_𝑁∞^1/4+𝐴𝑢_𝑁∞𝐷⁴𝑢_𝑁≤𝐶𝐷²𝑃_𝑁𝐴𝑢_𝑁^3/4+1𝐷⁴𝑢_𝑁≤𝐶𝐷²𝑃_𝑁𝐴𝑢_𝑁+1𝐷⁴𝑢_𝑁≤𝐶𝐷²𝐴𝑢_𝑁+1𝐷⁴𝑢_𝑁≤𝐶𝐴^′𝑢_𝑁𝐷²𝑢_𝑁+𝐴^′′𝑢_𝑁𝐷𝑢_𝑁²⋅𝐷⁴𝑢_𝑁≤𝜀𝐷⁴𝑢_𝑁²+𝐶. (2.32) Hence, 12𝑑𝐷𝑑𝑡²𝑢_𝑁²+𝑚₀𝐷⁴𝑢_𝑁²≤4𝜀𝐷⁴𝑢_𝑁²+𝐶𝐷²𝑢_𝑁²+𝐶. (2.33) Taking 𝜀=𝑚₀/8, we have 12𝑑𝐷𝑑𝑡²𝑢_𝑁²+𝑚₀2𝐷⁴𝑢_𝑁²≤𝐶𝐷²𝑢_𝑁²+𝐶, (2.34) where 𝐶=𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂)>0 is a positive constant. Therefore, 12𝑑𝐷𝑑𝑡²𝑢_𝑁²≤𝐶𝐷²𝑢_𝑁²+𝐶. (2.35) From Gronwall inequality it follows that 𝐷²𝑢_𝑁²≤exp^𝐶𝑡𝐷²𝑢₀²+1≤exp^𝐶𝑇𝐷²𝑢₀²+1≤𝐶, (2.36) where 𝐶=𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂)>0 is a positive constant. According to the embedding theorem, we have 𝐷𝑢_𝑁∞≤𝐶,0≤𝑡≤𝑇, (2.37) where 𝐶=𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂)>0 is a positive constant.

Now, we study the error estimates between the exact solution 𝑢 and the semidiscrete spectral approximation solution 𝑢_𝑁. Set the following decomposition: 𝑢−𝑢_𝑁=𝜂+𝑒,𝜂=𝑢−𝑃_𝑁𝑢,𝑒=𝑃_𝑁𝑢−𝑢_𝑁. (2.38) From the inequality (2.6) it follows that 𝜂≤𝐶𝑁⁻⁴. (2.39) Hence, it remains to obtain the approximate bounds of 𝑒.

Theorem 2.6. Assume that 𝑢 is the solution of (1.2)–(1.4), 𝑢_𝑁∈𝑆_𝑁 is the solution of (2.10) and (2.11), and 𝑚 is smooth and satisfies (1.6), then there exists a constant 𝐶=𝐶(𝑢₀,𝑚,𝛾₁,𝛾₂) such that ‖𝑒‖≤𝐶𝑁⁻²+‖𝑒(0)‖. (2.40)

Before we prove this theorem, we study some useful approximation properties.

Lemma 2.7. For any 𝜀>0, we have −𝐷²𝜂+𝐷²𝑒,𝐷𝑚𝑢_𝑁𝐷𝑒≤−𝑚₀𝐷²𝑒²+3𝜀𝐷²𝑒²+𝐶‖𝑒‖²+𝑁⁻⁴, (2.41) where 𝐶=𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂,𝜀)>0 is a positive constant.

Proof. Direct computation gives −𝐷²𝜂+𝐷²𝑒,𝐷𝑚𝑢_𝑁𝐷𝐷𝑒=−²𝑒,𝑚𝑢_𝑁𝐷²𝑒−𝐷²𝜂,𝑚𝑢_𝑁𝐷²𝑒−𝐷²𝜂,𝑚^′𝑢_𝑁𝐷𝑢_𝑁𝐷𝑒−𝐷²𝑒,𝑚^′𝑢_𝑁𝐷𝑢_𝑁𝐷𝑒≤−𝑚₀𝐷²𝑒²+𝑀₀𝐷²𝑒𝐷²𝜂²+𝑀₁𝐷²𝜂𝐷𝑢_𝑁∞‖𝐷𝑒‖+𝑀₁𝐷²𝑒𝐷𝑢_𝑁∞‖𝐷𝑒‖≤−𝑚₀𝐷²𝑒²+𝜀𝐷²𝑒²+𝑀²₀𝐷4𝜀²𝜂²+𝑀₁𝐷𝑢_𝑁∞‖𝐷𝑒‖²+𝐷²𝜂²+𝜀𝐷²𝑒²+𝑀²₁𝐷𝑢_𝑁²_∞‖‖4𝜀𝐷𝑒²≤−𝑚₀𝐷²𝑒²+2𝜀𝐷²𝑒²+𝑀²₀4𝜀+𝑀₁𝐷𝑢_𝑁∞𝐷²𝜂²+𝑀₁𝐷𝑢_𝑁∞+𝑀²₁𝐷𝑢_𝑁²_∞‖‖4𝜀𝐷𝑒²≤−𝑚₀𝐷²𝑒²+2𝜀𝐷²𝑒²+𝑀²₀4𝜀+𝑀₁𝐷𝑢_𝑁∞𝐷²𝜂²+𝑀₁𝐷𝑢_𝑁∞+𝑀²₁𝐷𝑢_𝑁²_∞𝐷4𝜀²𝑒⋅‖𝑒‖≤−𝑚₀𝐷²𝑒²+2𝜀𝐷²𝑒²+𝐶𝐷²𝜂²+𝜀𝐷²𝑒²+𝐶‖𝑒‖²≤−𝑚₀𝐷²𝑒²+3𝜀𝐷²𝑒²+𝐶‖𝑒‖²+𝐷²𝜂²≤−𝑚₀𝐷²𝑒²+3𝜀𝐷²𝑒²+𝐶‖𝑒‖²+𝑁⁻⁴, (2.42) where 𝐶=𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂,𝜀)>0 is a positive constant.

Lemma 2.8. For any 𝜀>0, we have 𝐴(𝑢)−𝑃_𝑁𝐴𝑢_𝑁,𝐷𝑚𝑢_𝑁𝐷𝑒≤3𝜀𝑀²₀𝐷²𝑒²+𝐶‖𝑒‖²+𝑁⁻⁴, (2.43) where 𝐶=𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂,𝜀)>0 is a positive constant.

Proof. Noticing that 𝐴(𝑢)−𝑃_𝑁𝐴𝑢_𝑁,𝐷𝑚𝑢_𝑁=𝐴(𝑢)𝐷𝑒−𝑃_𝑁𝐴(𝑢),𝐷𝑚𝑢_𝑁𝐷𝑒+𝑃_𝑁𝐴(𝑢)−𝑃_𝑁𝐴𝑢_𝑁,𝐷𝑚𝑢_𝑁≤𝐴(𝑢)𝐷𝑒−𝑃_𝑁𝐴(𝑢)⋅𝐷𝑚𝑢_𝑁𝐷𝑒+𝑃_𝑁𝐴𝑢,𝑢_𝑁𝑢−𝑢_𝑁⋅𝐷𝑚𝑢_𝑁≤𝐴(𝑢)𝐷𝑒−𝑃_𝑁𝐴(𝑢)⋅𝐷𝑚𝑢_𝑁𝐷𝑒+𝐴𝑢,𝑢_𝑁𝜂+𝑒⋅𝐷𝑚𝑢_𝑁𝐷𝑒, (2.44) where 𝐴(𝑠,𝜏)=𝛾₂𝑠²+𝑠𝜏+𝜏²+𝛾₁(𝑠+𝜏)−1. (2.45) From the boundedness of 𝑢_𝑁 in Theorem 2.4 and the property of 𝑢 in Lemma 2.3, it follows that 𝐴𝑢,𝑢_𝑁∞≤𝐶. (2.46) Then we obtain 𝐴(𝑢)−𝑃_𝑁𝐴(𝑢)²≤𝐶𝑢³−𝑃_𝑁𝑢³²+𝑢²−𝑃_𝑁𝑢²²+𝑢−𝑃_𝑁𝑢²≤𝐶𝑁⁻⁸𝑢³₄+𝑢²₄+|𝑢|₄²≤𝐶𝑁⁻⁸,𝐴𝑢,𝑢_𝑁𝜂+𝑒²≤𝐴𝑢,𝑢_𝑁²_∞‖𝑒‖²+𝜂²≤𝐶𝑁⁻⁸+‖𝑒‖²,𝐷𝑚𝑢_𝑁𝐷𝑒²≤2𝑚𝑢_𝑁𝐷²𝑒²+2𝑚^′𝑢_𝑁𝐷𝑢_𝑁𝐷𝑒²≤2𝑀²₀𝐷²𝑒²+2𝑀²₁𝐷𝑢_𝑁²_∞𝐷²𝑒²‖𝑒‖²≤2𝑀²₀𝐷²𝑒²+𝑀²₀𝐷²𝑒²+𝑀⁴₁𝐷𝑢_𝑁⁴_∞𝑀²₀‖𝑒‖²≤3𝑀²₀𝐷²𝑒²+𝐶‖𝑒‖², (2.47) where 𝐶=𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂)>0 is a positive constant. By Cauchy inequality, for any 𝜀>0, we have 𝐴(𝑢)−𝑃_𝑁𝐴𝑢_𝑁,𝐷𝑚𝑢_𝑁≤𝐴(𝑢)𝐷𝑒−𝑃_𝑁𝐴(𝑢)⋅𝐷𝑚𝑢_𝑁𝐷𝑒+𝐴𝑢,𝑢_𝑁𝜂+𝑒⋅𝐷𝑚𝑢_𝑁≤𝜀𝐷𝑒2‖𝐷(𝑚(𝑢)𝐷𝑒)‖²+1𝐴(𝑢)2𝜀−𝑃_𝑁𝐴(𝑢)²+𝜀2‖𝐷(𝑚(𝑢)𝐷𝑒)‖²+12𝜀𝐴𝑢,𝑢_𝑁𝜂+𝑒²≤3𝜀𝑀²₀𝐷²𝑒²+𝐶‖𝑒‖²+𝑁⁻⁸, (2.48) where 𝐶=𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂,𝜀)>0 is a positive constant.

Lemma 2.9. Assume that 𝑢 is the solution of (1.2)–(1.4), there exists a positive constant 𝐶=𝐶(𝑢₀,𝑚,𝛾₁,𝛾₂) such that 𝐷²𝑢−𝐴(𝑢),𝐷𝑚(𝑢)−𝑚𝑢_𝑁≤𝑚𝐷𝑒⁰2𝐷²𝑒²+𝐶‖𝑒‖²+𝑁⁻⁴, (2.49) where 𝐶=𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂)>0 is a positive constant.

Proof. By Lemma 2.3, we have 𝐷³𝑢−𝐷𝐴(𝑢)_∞≤𝐷³𝑢_∞+𝐴^′(𝑢)𝐷𝑢_∞≤𝐶. (2.50) In the other hand, it follows that 𝑚(𝑢)−𝑚𝑢_𝑁²≤𝑀²₁𝑢−𝑢_𝑁²≤𝑀²₁‖𝑒‖²+𝜂². (2.51) By the Young inequality, 𝐷²𝑢−𝐴(𝑢),𝐷𝑚(𝑢)−𝑚𝑢_𝑁𝐷𝑒=−𝐷³𝑢−𝐷𝐴(𝑢),𝑚(𝑢)−𝑚𝑢_𝑁≤𝐷𝐷𝑒₃𝑢−𝐷𝐴(𝑢)_∞𝑚(𝑢)−𝑚𝑢_𝑁⋅‖‖𝑚(𝑢)𝑢𝐷𝑒≤𝐶−𝑚_𝑁‖‖‖‖𝐷𝑒≤𝜀𝐷𝑒²+𝐶_𝜀‖𝑒‖²+𝜂²≤𝜀𝐷²𝑒²+𝐶_𝜀‖𝑒‖²+𝜂². (2.52) Choosing 𝜀=𝑚₀/2 in the previous inequality, we obtain (2.49).

Proof of Theorem 2.6. Setting 𝑣=𝑒 in (2.2), we obtain 𝜕𝑢+𝐷𝜕𝑡,𝑒²𝑢−𝐴(𝑢),𝐷(𝑚(𝑢)𝐷𝑒)=0. (2.53) Setting 𝑣_𝑁=𝑒 in (2.10), we get 𝜕𝑢_𝑁+𝐷𝜕𝑡,𝑒²𝑢_𝑁−𝑃_𝑁𝐴𝑢_𝑁,𝐷𝑚𝑢_𝑁𝐷𝑒=0. (2.54) (2.53) minus (2.54) gives 𝜕𝑒𝜕𝑡,𝑒=−𝐷𝑢²−𝐴(𝑢),𝐷𝑚(𝑢)−𝑚𝑢_𝑁𝐷𝑒−𝐷²𝑢−𝐷²𝑢_𝑁,𝐷𝑚𝑢_𝑁+𝐴(𝑢)𝐷𝑒−𝑃_𝑁𝐴𝑢_𝑁,𝐷𝑚𝑢_𝑁𝐷𝑒. (2.55) According to Lemmas 2.7, 2.8 and 2.9, we have 12𝑑‖𝑒‖𝑑𝑡²≤−𝑚₀𝐷²𝑒²+3𝜀𝐷²𝑒²+𝐶‖𝑒‖²+𝑁⁻⁴+3𝜀𝑀²₀𝐷²𝑒²+𝐶𝑁⁻⁸+‖𝑒‖²≤−𝑚₀+4+3𝑀²₀𝜀𝐷²𝑒²+𝐶𝑁⁻⁴+‖𝑒‖². (2.56) Set 𝜀=𝑚₀/(4+3𝑀²₀), then there exists a positive constant 𝐶=𝐶(𝑢₀,𝑚,𝛾₁,𝛾₂) such that 𝑑‖𝑒‖𝑑𝑡²≤𝐶‖𝑒‖²+𝑁⁻⁴. (2.57) By Gronwall inequality, we have ‖𝑒‖≤exp^𝐶𝑡𝑒₀+𝑁⁻²≤exp^𝐶𝑇𝑒₀+𝑁⁻²≤𝐶𝑒₀+𝑁⁻², (2.58) where 𝐶=𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂)>0 is a constant.

Summig up the properties above, we obtain the following.

Theorem 2.10. Assume 𝑚(𝑠) is sufficiently smooth and satisfies (1.6), 𝑢 is the solution of (1.2)–(1.4), and 𝑢_𝑁 is the solution of (2.10) and (2.11). Then there exists a positive constant 𝐶=𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂)>0 such that 𝑢−𝑢_𝑁≤𝐶𝑁⁻²+𝑢₀−𝑢_0𝑁,∀𝑡∈(0,𝑇). (2.59)

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,𝑇]. Let 0=𝑡₀<𝑡₁<⋯<𝑡_Λ, where 𝑡_𝑗=𝑗ℎ and ℎ=𝑇/Λ is time-step size. Then the full-discretization spectral method for (1.2)–(1.4) reads: ∀𝑣∈𝑆_𝑁, find 𝑈^𝑗_𝑁∈𝑆_𝑁(𝑗=0,1,2,…,Λ) such that 𝑈_𝑁^𝑗+1−𝑈^𝑗_𝑁ℎ,𝑣+𝐷²𝑈_𝑁^𝑗+1/2−𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁,𝐷𝑚𝑈_𝑁^𝑗+1/2𝑈𝐷𝑣=0,_0𝑁,𝑣−𝑢₀,𝑣=0, (3.1)(3.2) where 𝑈_𝑁^𝑗+1/2=(𝑈_𝑁^𝑗+1+𝑈^𝑗_𝑁)/2 and 𝐴=𝛾𝜙,𝜑²4𝜙³+𝜑³+𝜙²𝜑+𝜙𝜑²+𝛾₁3𝜙²+𝜙𝜑+𝜑²−12𝜙+𝜑. (3.3) The solution 𝑈^𝑗_𝑁 has the following property.

Lemma 3.1. Assume that 𝑈^𝑗_𝑁∈𝑆_𝑁(𝑗=1,2,…,Λ) is the solution of (3.1)-(3.2). Then there exists a constant 𝐶=𝐶(𝑢₀,𝑚,𝛾₁,𝛾₂)>0 such that 𝑈^𝑗_𝑁∞≤𝐶,𝐷𝑈^𝑗_𝑁≤𝐶. (3.4)

Proof. Define a discrete energy function at time 𝑡_𝑗 by 𝐹𝑗=12𝐷𝑈^𝑗_𝑁²+𝐻𝑈^𝑗_𝑁,1. (3.5) Notice that 1ℎ𝐹𝑗+1−𝐹𝑗=𝐷𝑈_𝑁^𝑗+1/2𝑈,𝐷^𝑁𝑗+1−𝑈^𝑗_𝑁ℎ+𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁,𝑈_𝑁^𝑗+1−𝑈^𝑗_𝑁ℎ=−𝐷²𝑈_𝑁^𝑗+1/2−𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁,𝑈_𝑁^𝑗+1−𝑈^𝑗_𝑁ℎ. (3.6) Setting 𝑣=𝐷²𝑈_𝑁^𝑗+1/2−𝑃_𝑁𝐴(𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁)∈𝑆_𝑁 in (3.1), we obtain 1ℎ𝐹𝑗𝑚𝑗+1−𝐹=−𝑈_𝑁^𝑗+1/2𝐷³𝑈_𝑁^𝑗+1/2−𝐷𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁,𝐷³𝑈_𝑁^𝑗+1/2−𝐷𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁≤−𝑚₀𝐷³𝑈_𝑁^𝑗+1/2−𝐷𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁²≤0, (3.7) which implies 𝐹𝑗(0)=1≤𝐹2𝐷𝑈⁰_𝑁²+𝐻𝑈⁰_𝑁,1. (3.8) By (2.18), we have ¹₀𝐻𝑈^𝑗_𝑁𝛾𝑑𝑥≥²8¹₀𝑈^𝑗_𝑁⁴𝑑𝑥−𝐾₁. (3.9) Then 12𝐷𝑈^𝑗_𝑁²+𝛾₂8¹₀𝑈^𝑗_𝑁⁴𝑑𝑥−𝐾₁≤12𝐷𝑈^𝑗_𝑁²+¹₀𝐻𝑈^𝑗_𝑁𝑗(0)=1𝑑𝑥=𝐹≤𝐹2𝐷𝑈⁰_𝑁²+𝐻𝑈⁰_𝑁,1. (3.10) So we obtain 𝐷𝑈^𝑗_𝑁²≤2𝐾₁+𝐷𝑈⁰_𝑁²+2𝐻𝑈⁰_𝑁,1≡𝐶₁,₁₀𝑈^𝑗_𝑁⁴8𝑑𝑥≤𝛾₂𝐾₁+12𝐷𝑈⁰_𝑁²+𝐻𝑈⁰_𝑁,1≡𝐶₂, (3.11) where 𝐶₁=𝐶₁(𝑢₀,𝑚,𝛾₁,𝛾₂) and 𝐶₂=𝐶₂(𝑢₀,𝑚,𝛾₁,𝛾₂) are positive constants. By Hölder inequality, we get 𝑈^𝑗_𝑁²=¹₀𝑈^𝑗_𝑁²𝑑𝑥≤¹₀𝑈^𝑗_𝑁⁴𝑑𝑥^1/2≤𝐶₂^1/2. (3.12) Therefore, 𝑈^𝑗_𝑁≤𝐶₂^1/4. (3.13) By the embedding theorem, we obtain 𝑈^𝑗_𝑁∞≤𝐶, (3.14) where 𝐶=𝐶(𝑢₀,𝑚,𝛾₁,𝛾₂)>0 is a constant.

Lemma 3.2. Assume that 𝑈^𝑗_𝑁∈𝑆_𝑁(𝑗=1,2,…,Λ) is the solution of the full-discretization scheme (3.1)-(3.2), then there is a constant 𝐶=𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂)>0 such that 𝐷𝑈^𝑗_𝑁∞≤𝐶,𝐷²𝑈^𝑗_𝑁≤𝐶. (3.15)

Proof. Setting 𝑣=2𝐷⁴𝑈_𝑁^𝑗+1/2∈𝑆_𝑁 in (3.1), we have 𝑈_𝑁^𝑗+1−𝑈^𝑗_𝑁ℎ,2𝐷⁴𝑈_𝑁^𝑗+1/2+𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷³𝑈_𝑁^𝑗+1/2−𝐷𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁2𝐷⁴𝑈_𝑁^𝑗+1/2=𝐷²𝑈_𝑁^𝑗+1−𝐷²𝑈^𝑗_𝑁ℎ,𝐷²𝑈_𝑁^𝑗+1+𝐷²𝑈^𝑗_𝑁+𝑚𝑈_𝑁^𝑗+1/2𝐷⁴𝑈_𝑁^𝑗+1/2,2𝐷⁴𝑈_𝑁^𝑗+1/2+𝑚^′𝑈_𝑁^𝑗+1/2𝐷𝑈_𝑁^𝑗+1/2𝐷³𝑈_𝑁^𝑗+1/2,2𝐷⁴𝑈_𝑁^𝑗+1/2−𝑚𝑈_𝑁^𝑗+1/2𝐷²𝑃_𝑁𝐴𝑈_𝑁^𝑗+1+𝑈^𝑗_𝑁,2𝐷⁴𝑈_𝑁^𝑗+1/2−𝑚^′𝑈_𝑁^𝑗+1/2𝐷𝑈_𝑁^𝑗+1/2𝐷𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁,2𝐷⁴𝑈_𝑁^𝑗+1/2=0. (3.16) Therefore, 1ℎ𝐷²𝑈_𝑁^𝑗+12−𝐷²𝑈^𝑗_𝑁²+2𝑚𝑈_𝑁^𝑗+1/2𝐷⁴𝑈_𝑁^𝑗+1/2,𝐷⁴𝑈_𝑁^𝑗+1/2=2𝑚^′𝑈_𝑁^𝑗+1/2𝐷³𝑈_𝑁^𝑗+1/2𝐷𝑈_𝑁^𝑗+1/2,𝐷⁴𝑈_𝑁^𝑗+1/2+2𝑚𝑈_𝑁^𝑗+1/2𝐷²𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁,𝐷⁴𝑈_𝑁^𝑗+1/2+2𝑚^′𝑈_𝑁^𝑗+1/2𝐷𝑈_𝑁^𝑗+1/2𝐷𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈_𝑁^𝑗+1,𝐷⁴𝑈_𝑁^𝑗+1/2=𝐼^𝑗₁+𝐼^𝑗₂+𝐼^𝑗₃. (3.17) By Nirenberg inequality (2.7), we have 𝐷𝑈^𝑗_𝑁∞≤𝐶𝐷⁴𝑈^𝑗_𝑁^3/8𝑈^𝑗_𝑁^5/8+𝑈^𝑗_𝑁,𝐷𝑈^𝑗_𝑁∞≤𝐶𝐷²𝑈^𝑗_𝑁^3/4𝑈^𝑗_𝑁^1/4+𝑈^𝑗_𝑁,𝐷³𝑈^𝑗_𝑁∞≤𝐶𝐷⁴𝑈^𝑗_𝑁^7/8𝑈^𝑗_𝑁^1/8+𝑈^𝑗_𝑁. (3.18) According to (3.14), we obtain 𝐴^′₁𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁∞≤𝐶,𝐴^′₂𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁∞𝐴≤𝐶,_′′11𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁∞≤𝐶,𝐴^′′₁₂𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁∞≤𝐶,𝐴^′′₂₂𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁∞≤𝐶, (3.19) where 𝐶=𝐶(𝑢₀,𝑚,𝛾₁,𝛾₂) is a positive constant. By Young inequality, for any positive constant 𝜀>0, it follows that 𝐼^𝑗₁≤2𝑀₁𝐷³𝑈_𝑁^𝑗+1/2_∞𝐷𝑈_𝑁^𝑗+1/2𝐷⁴𝑈_𝑁^𝑗+1/2≤𝐶𝐷⁴𝑈_𝑁^𝑗+1/27/8+1𝐷⁴𝑈_𝑁^𝑗+1/2≤𝐶𝐷⁴𝑈_𝑁^𝑗+1/215/8+𝐶𝐷⁴𝑈_𝑁^𝑗+1/2≤𝜀2𝐷⁴𝑈_𝑁^𝑗+1/22+𝐶_𝜀+𝜀2𝐷⁴𝑈_𝑁^𝑗+1/22+𝐶_𝜀≤𝜀𝐷⁴𝑈_𝑁^𝑗+1/22+𝐶_𝜀,𝐼^𝑗₂≤2𝑀₀𝐷²𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁𝐷⁴𝑈_𝑁^𝑗+1/2≤2𝑀₀𝐴^′′₁₁𝐷𝑈_𝑁^𝑗+12+2𝐴^′′₁₂𝐷𝑈_𝑁^𝑗+1𝐷𝑈^𝑗_𝑁+𝐴^′′₂₂𝐷𝑈^𝑗_𝑁²+𝐴^′₁𝐷²𝑈_𝑁^𝑗+1+𝐴^′₂𝐷²𝑈^𝑗_𝑁𝐷⁴𝑈_𝑁^𝑗+1/2≤𝐶𝐷²𝑈_𝑁^𝑗+1+𝐷²𝑈^𝑗_𝑁+𝐷𝑈_𝑁^𝑗+12+𝐷𝑈^𝑗_𝑁²𝐷⁴𝑈_𝑁^𝑗+1/2≤𝐶𝐷²𝑈_𝑁^𝑗+1+𝐷²𝑈^𝑗_𝑁+𝐷𝑈_𝑁^𝑗+1𝐷𝑈_𝑁^𝑗+1_∞+𝐷𝑈^𝑗_𝑁𝐷𝑈^𝑗_𝑁∞𝐷⁴𝑈_𝑁^𝑗+1/2≤𝐶𝐷²𝑈_𝑁^𝑗+1+𝐷²𝑈^𝑗_𝑁+𝐷𝑈_𝑁^𝑗+1_∞+𝐷𝑈^𝑗_𝑁∞𝐷⁴𝑈_𝑁^𝑗+1/2≤𝐶𝐷²𝑈_𝑁^𝑗+1+𝐷²𝑈^𝑗_𝑁+𝐷²𝑈_𝑁^𝑗+13/4+𝐷²𝑈^𝑗_𝑁^3/4+1𝐷⁴𝑈_𝑁^𝑗+1/2≤𝐶𝐷²𝑈_𝑁^𝑗+1+𝐷²𝑈^𝑗_𝑁+1𝐷⁴𝑈_𝑁^𝑗+1/2≤𝜀𝐷⁴𝑈_𝑁^𝑗+1/22+𝐶_𝜀𝐷²𝑈_𝑁^𝑗+12+𝐷²𝑈^𝑗_𝑁²,𝐼+1_𝑗3≤𝑀₁𝐷𝑈_𝑁^𝑗+1/2𝐷𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁∞𝐷⁴𝑈_𝑁^𝑗+1/2≤𝐶𝐷²𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁^3/4+1𝐷⁴𝑈_𝑁^𝑗+1/2≤𝐶𝐷²𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁+1𝐷⁴𝑈_𝑁^𝑗+1/2≤𝐶𝐷²𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁𝐷⁴𝑈_𝑁^𝑗+1/2+𝐷⁴𝑈_𝑁^𝑗+1/2≤𝜀𝐷⁴𝑈_𝑁^𝑗+1/22+𝐶_𝜀𝐷²𝑈_𝑁^𝑗+12+𝐷²𝑈^𝑗_𝑁²+1+𝜀𝐷⁴𝑈_𝑁^𝑗+1/22+𝐶_𝜀≤2𝜀𝐷⁴𝑈_𝑁^𝑗+1/22+𝐶_𝜀𝐷²𝑈_𝑁^𝑗+12+𝐷²𝑈^𝑗_𝑁²+1, (3.20) where 𝐶_𝜀=𝐶_𝜀(𝑢₀,𝑚,𝛾₁,𝛾₂)>0 is a constant. Therefore, 1ℎ𝐷²𝑈_𝑁^𝑗+12−𝐷²𝑈^𝑗_𝑁²+2𝑚₀𝐷⁴𝑈_𝑁^𝑗+1/22≤1ℎ𝐷²𝑈_𝑁^𝑗+12−𝐷²𝑈^𝑗_𝑁²+2𝑚𝑈_𝑁^𝑗+1/2𝐷⁴𝑈_𝑁^𝑗+1/2,𝐷⁴𝑈_𝑁^𝑗+1/2≤4𝜀𝐷²𝑈_𝑁^𝑗+1/22+𝐶_𝜀𝐷²𝑈_𝑁^𝑗+12+𝐷²𝑈^𝑗_𝑁²+1. (3.21) Setting 𝜀=𝑚₀/2, there is a positive constant 𝐶=𝐶(𝑢₀,𝑚,𝛾₁,𝛾₂) such that 𝐷²𝑈_𝑁^𝑗+12≤1+2𝐶ℎ𝐷1−𝐶ℎ²𝑈^𝑗_𝑁²+𝐶1−𝐶ℎℎ. (3.22) Denoted by 𝐶=𝐶/(1−𝐶ℎ), if ℎ is sufficiently small such that 𝐶ℎ≤1/2, we have 𝐷²𝑈_𝑁^𝑗+12≤1+𝐶ℎ𝐷²𝑈^𝑗_𝑁²+𝐶2ℎ≤1+𝐶ℎ^𝑗+1𝐷²𝑈⁰_𝑁²+𝐶2ℎ^𝑗_𝑖=11+𝐶ℎ^𝑖≤exp^(𝑗+1)𝐶_ℎ𝐷²𝑈⁰_𝑁²+𝐶𝐶𝑗ℎ≤exp^𝑇𝐷²𝑈⁰_𝑁²+𝐶𝑇≡𝐶^′, (3.23) where 𝐶^′=𝐶^′(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂)>0 is a positive constant. By the embedding theorem, the estimate (3.15) holds.

Next, we investigate the error estimates for the numerical solution 𝑈^𝑗_𝑁 to the exact solution 𝑢(𝑡_𝑗). Our analysis is based on the error decomposition denoted by 𝑢𝑡_𝑗−𝑈^𝑗_𝑁=𝜂^𝑗+𝑒^𝑗,𝜂^𝑗=𝑢𝑡_𝑗−𝑃_𝑁𝑢𝑡_𝑗,𝑒^𝑗=𝑃_𝑁𝑢𝑡_𝑗−𝑈^𝑗_𝑁. (3.24) The boundedness estimate of 𝜂^𝑗 follows from the inequality (2.6), that is, for any 0≤𝑗≤Λ, there is a positive constant 𝐶=𝐶(𝑢₀,𝑚,𝛾₁,𝛾₂) such that 𝜂^𝑗≤𝐶𝑁⁻⁴. (3.25) Hence, it remains to obtain the approximate bounds of 𝑒^𝑗. If no confusion occurs, we denote the average of the two instant errors 𝑒^𝑗 and 𝑒^𝑗+1 by 𝑒^𝑗+1/2: 𝑒^𝑗+1/2=𝑒^𝑗+𝑒^𝑗+12,𝜂^𝑗+1/2=𝜂^𝑗+𝜂^𝑗+12. (3.26) 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 𝑢_𝑡𝑡𝑡∈𝐿²(𝑄_𝑇), then 𝑒^𝑗+12≤𝑒^𝑗2+2ℎ𝑢_𝑡𝑡_𝑗+1/2−𝑈_𝑁^𝑗+1−𝑈^𝑗_𝑁ℎ,𝑒^𝑗+1/2+1ℎ320^4𝑡𝑗+1_𝑡𝑗𝑢_𝑡𝑡𝑡²𝑑𝑡+ℎ𝑒^𝑗+1/22. (3.27)

Proof. Applying Taylor expansion about 𝑡_𝑗+1/2, we have 𝑢^𝑗=𝑢^𝑗+1/2−ℎ2𝑢_𝑡^𝑗+1/2+ℎ²8𝑢^𝑗+1/2_𝑡𝑡−12^𝑡𝑗+1/2_𝑡𝑗𝑡−𝑡_𝑗²𝑢_𝑡𝑡𝑡𝑢𝑑𝑡,_𝑗+1=𝑢^𝑗+1/2+ℎ2𝑢_𝑡^𝑗+1/2+ℎ²8𝑢^𝑗+1/2_𝑡𝑡+12^𝑡𝑗+1_𝑡𝑗+1/2𝑡_𝑗+1−𝑡²𝑢_𝑡𝑡𝑡𝑑𝑡. (3.28) Then 𝑢_𝑡^𝑗+1/2−𝑢^𝑗+1−𝑢^𝑗ℎ1=−2ℎ^𝑡𝑗+1_𝑡𝑗+1/2𝑡_𝑗+1−𝑡²𝑢_𝑡𝑡𝑡𝑑𝑡+^𝑡𝑗+1/2_𝑡𝑗𝑡−𝑡_𝑗²𝑢_𝑡𝑡𝑡𝑑𝑡. (3.29) From Hölder inequality it follows that 𝑢_𝑡^𝑗+1/2−𝑢^𝑗+1−𝑢^𝑗ℎ²≤12ℎ^{2𝑡𝑗+1𝑡}_𝑗^+1/2𝑡_𝑗+1−𝑡²𝑢_𝑡𝑡𝑡𝑑𝑡²+^𝑡𝑗+1/2_𝑡𝑗𝑡−𝑡_𝑗²𝑢_𝑡𝑡𝑡𝑑𝑡²≤ℎ³320^𝑡𝑗+1_𝑡𝑗𝑢_𝑡𝑡𝑡²𝑑𝑡. (3.30) Noticing that for any 𝑣∈𝑆_𝑁, we have 𝑢_𝑡^𝑗+1/2−𝑈_𝑁^𝑗+1−𝑈^𝑗_𝑁ℎ,𝑣=𝑢_𝑡^𝑗+1/2−𝑢^𝑗+1−𝑢^𝑗ℎ+1,𝑣ℎ𝑒^𝑗+1−𝑒^𝑗,𝑣. (3.31) Taking 𝑣=2𝑒^𝑗+1/2 in (3.31), we obtain 𝑒^𝑗+12=𝑒^𝑗2+2ℎ𝑢_𝑡^𝑗+1/2−𝑈_𝑁^𝑗+1−𝑈^𝑗_𝑁ℎ,𝑒^𝑗+1/2−2ℎ𝑢_𝑡^𝑗+1/2−𝑢^𝑗+1−𝑢^𝑗,Δ𝑡𝑒^𝑗+1/2≤𝑒^𝑗2+2ℎ𝑢_𝑡^𝑗+1/2−𝑈_𝑁^𝑗+1−𝑈^𝑗_𝑁ℎ,𝑒^𝑗+1/2+1ℎ320^4𝑡𝑗+1_𝑡𝑗𝑢_𝑡𝑡𝑡²𝑑𝑡+ℎ𝑒^𝑗+1/22. (3.32)

Taking 𝑣=𝑒^𝑗+1/2 in (2.2) and (3.1), respectively, we have 𝜕𝑢𝑡_𝑗+1/2,𝜕𝑡𝑒^𝑗+1/2+𝐷²𝑢^𝑗+1/2−𝐴𝑢^𝑗+1/2,𝐷𝑚𝑢^𝑗+1/2𝐷𝑒^𝑗+1/2𝑈=0,_𝑁𝑗+1−𝑈^𝑗_𝑁ℎ,𝑒^𝑗+1/2+𝐷²𝑈_𝑁^𝑗+1/2−𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁,𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2=0. (3.33)(3.34) Comparing (3.33) and (3.34), we have 𝑢_𝑡^𝑗+1/2−𝑈_𝑁^𝑗+1−𝑈^𝑗_𝑁ℎ,𝑒^𝑗+1/2=−𝐷²𝑢^𝑗+1/2−𝐷²𝑈_𝑁^𝑗+1+𝐷²𝑈^𝑗_𝑁2,𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2+𝐴𝑢^𝑗+1/2−𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁,𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2+𝐷²𝑢^𝑗+1/2−𝐴𝑢^𝑗+1/2,𝐷𝑚𝑈_𝑁^𝑗+1/2−𝑚𝑢^𝑗+1/2𝐷𝑒^𝑗+1/2. (3.35) Now we investigate the error estimates of the three items in the right-hand side of the previous equation.

Lemma 3.4. Assume that 𝑢 is the solution of (1.2)–(1.4) such that 𝐷𝑢_𝑡𝑡∈𝐿²(𝑄_𝑇), then there exists a positive constant 𝐶₁=𝐶₁(𝑚,𝑢₀,𝑇,𝛾₁,𝛾₂)>0 such that −𝐷²𝑢^𝑗+1/2−𝐷²𝑈_𝑁^𝑗+1+𝐷²𝑈^𝑗_𝑁2,𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2𝑚≤−⁰2𝐷²𝑒^𝑗+1/22+𝐶₁𝑁⁻⁴+𝑒^𝑗2+𝑒^𝑗+12+ℎ^3𝑡𝑗+1_𝑡𝑗𝐷²𝑢_𝑡𝑡²𝑑𝑡. (3.36)

Proof. By Taylor expansion and Hölder inequality, we obtain 𝑢^𝑗=𝑢^𝑗+1/2−ℎ2𝑢_𝑡^𝑗+1/2+^𝑡𝑗+1/2_𝑡𝑗𝑡−𝑡_𝑗𝑢_𝑡𝑡𝑢𝑑𝑡,_𝑗+1=𝑢^𝑗+1/2+ℎ2𝑢_𝑡^𝑗+1/2+^𝑡𝑗+1_𝑡𝑗+1/2𝑡_𝑗+1−𝑡𝑢_𝑡𝑡𝑑𝑡. (3.37) Therefore, 12𝑢^𝑗+𝑢^𝑗+1−𝑢^𝑗+1/2=12^𝑡𝑗+1/2_𝑡𝑗𝑡−𝑡_𝑗𝑢_𝑡𝑡𝑑𝑡+^𝑡𝑗+1_𝑡𝑗+1/2𝑡_𝑗+1−𝑡𝑢_𝑡𝑡𝑑𝑡. (3.38) By Hölder inequality, we have 𝐷²𝑢^𝑗+1/2−12𝑢^𝑗+𝑢^𝑗+12=14𝐷^{2𝑡𝑗+1/2}_𝑡𝑗𝑡−𝑡_𝑗𝑢_𝑡𝑡𝑑𝑡+^𝑡𝑗+1_𝑡𝑗+1/2𝑡_𝑗+1−𝑡𝑢_𝑡𝑡𝑑𝑡²≤14^𝑡𝑗+1/2_𝑡𝑗𝑡−𝑡_𝑗²𝑑𝑡^𝑡𝑗+1/2_𝑡𝑗𝐷²𝑢_𝑡𝑡²𝑑𝑡+^𝑡𝑗+1_𝑡𝑗+1/2𝑡_𝑗+1−𝑡²𝑑𝑡^𝑡𝑗+1_𝑡𝑗+1/2𝐷²𝑢_𝑡𝑡²𝑑𝑡²=14ℎ³24^𝑡𝑗+1/2_𝑡𝑗𝐷²𝑢_𝑡𝑡²ℎ𝑑𝑡+³24^t𝑗+1_𝑡𝑗+1/2𝐷²𝑢_𝑡𝑡²𝑑𝑡²≤14⋅ℎ³24^𝑡𝑗+1_𝑡𝑗𝐷²𝑢_𝑡𝑡²𝑑𝑡. (3.39) Direct computation gives 𝐷²𝑢^𝑗+1/2−12𝑢^𝑗+𝑢^𝑗+12≤ℎ³96^𝑡𝑗+1_𝑡𝑗𝐷²𝑢_𝑡𝑡²𝑑𝑡. (3.40) Then −𝐷²𝑢^𝑗+1/2−𝑈_𝑁^𝑗+1+𝑈^𝑗_𝑁2,𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2=−𝐷²𝑢^𝑗+1/2−𝑢^𝑗+1+𝑢^𝑗2,𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2−𝐷²𝑢^𝑗+1+𝑢^𝑗2−𝑈_𝑁^𝑗+1+𝑈^𝑗_𝑁2,𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2≤𝐷²𝑢^𝑗+1/2−𝑢^𝑗+1+𝑢^𝑗2⋅𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2−𝐷²𝜂^𝑗+1/2+𝑒^𝑗+1/2,𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2≤ℎ³96^𝑡𝑗+1_𝑡𝑗𝐷²𝑢_𝑡𝑡²𝑑𝑡^1/2⋅𝑚𝑈_𝑁^𝑗+1/2𝐷²𝑒^𝑗+1/2+𝑚^′𝑈_𝑁^𝑗+1/2𝐷𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2−𝐷²𝜂^𝑗+1/2,𝑚𝑈_𝑁^𝑗+1/2𝐷²𝑒^𝑗+1/2−𝐷²𝜂^𝑗+1/2,𝑚^′𝑈_𝑁^𝑗+1/2𝐷𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2−𝐷²𝑒^𝑗+1/2,𝑚𝑈_𝑁^𝑗+1/2𝐷²𝑒^𝑗+1/2−𝐷²𝑒^𝑗+1/2,𝑚^′𝑈_𝑁^𝑗+1/2𝐷𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2≜𝜎^𝑗₁+𝜎^𝑗₂+𝜎^𝑗₃. (3.41) By Cauchy inequality, it follows that 𝜎^𝑗₁≤ℎ³96^𝑡𝑗+1_𝑡𝑗𝐷²𝑢_𝑡𝑡²𝑑𝑡^1/2⋅𝑚𝑈_𝑁^𝑗+1/2𝐷²𝑒^𝑗+1/2+ℎ³96^𝑡𝑗+1_𝑡𝑗𝐷²𝑢_𝑡𝑡²𝑑𝑡^1/2𝑚^′𝑈_𝑁^𝑗+1/2𝐷𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2≤𝑀₀ℎ³96^𝑡𝑗+1_𝑡𝑗𝐷²𝑢_𝑡𝑡²𝑑𝑡^1/2𝐷²𝑒^𝑗+1/2+𝑀₁𝐷𝑈_𝑁^𝑗+1/2_∞ℎ³96^𝑡𝑗+1_𝑡𝑗𝐷²𝑢_𝑡𝑡²𝑑𝑡^1/2𝐷𝑒^𝑗+1/2≤ℎ³𝐶_1𝜀^𝑡𝑗+1_𝑡𝑗𝐷²𝑢_𝑡𝑡²𝑑𝑡+𝜀𝐷²𝑒^𝑗+1/22+𝜀𝐷𝑒^𝑗+1/22≤ℎ³𝐶_1𝜀^𝑡𝑗+1_𝑡𝑗𝐷²𝑢_𝑡𝑡²𝑑𝑡+𝜀𝐷²𝑒^𝑗+1/22+𝜀𝐷²𝑒^𝑗+1/2𝑒^𝑗+1/2≤ℎ³𝐶_1𝜀^𝑡𝑗+1_𝑡𝑗𝐷²𝑢_𝑡𝑡²𝑑𝑡+2𝜀𝐷²𝑒^𝑗+1/22+𝜀𝑒^𝑗+1/22,𝜎^𝑗₂≤𝑀₀𝐷²𝜂^𝑗+1/2𝐷²𝑒^𝑗+1/2+𝑀₁𝐷𝑈_𝑁^𝑗+1/2_∞𝐷²𝜂^𝑗+1/2𝐷𝑒^𝑗+1/2≤𝑀²₀𝐷2𝜀²𝜂^𝑗+1/22+𝜀2𝐷²𝑒^𝑗+1/22+𝑀²₁𝐷𝑈_𝑁^𝑗+1/22_∞𝜀2𝜀2𝐷𝑒^𝑗+1/22≤𝐶_2𝜀‖𝐷²𝜂^𝑗+1/2‖²+𝜀2𝐷²𝑒^𝑗+1/22+𝜀2𝐷²𝑒^𝑗+1/2𝑒^𝑗+1/2≤𝐶_2𝜀‖𝐷²𝜂^𝑗+1/2‖²+𝜀‖𝐷²𝑒^𝑗+1/2‖²+𝜀2‖𝑒^𝑗+1/2‖²,𝜎^𝑗₃≤−𝑚₀𝐷²𝑒^𝑗+1/22+𝑀₁𝐷𝑈_𝑁^𝑗+1/2_∞𝐷²𝑒^𝑗+1/2⋅𝐷𝑒^𝑗+1/2≤−𝑚₀𝐷²𝑒^𝑗+1/22+𝜀2𝐷²𝑒^𝑗+1/22+𝑀²₁𝐷𝑈_𝑁^𝑗+1/22_∞𝐷2𝜀𝑒^𝑗+1/22≤−𝑚₀𝐷²𝑒^𝑗+1/22+𝜀𝐷²𝑒^𝑗+1/22+𝐶_3𝜀𝑒^𝑗+1/22. (3.42) Then we obtain −𝐷²𝑢^𝑗+1/2−𝐷²𝑈_𝑁^𝑗+1+𝐷²𝑈^𝑗_𝑁2,𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2≤𝜎^𝑗₁+𝜎^𝑗₂+𝜎^𝑗₃≤4𝜀−𝑚₀𝐷²𝑒^𝑗+1/22+ℎ³𝐶_1𝜀^𝑡𝑗+1_𝑡𝑗𝐷²𝑢_𝑡𝑡²𝑑𝑡+𝐶_2𝜀𝐷²𝜂^𝑗+1/22+𝐶_3𝜀𝑒^𝑗+1/22, (3.43) where 𝐶_1𝜀,𝐶_2𝜀, and 𝐶_3𝜀 are positive constants. Choosing 𝜀=𝑚₀/8, and terms of the properties of the projection operator 𝑃_𝑁, we complete the proof of the estimate (3.36).

Lemma 3.5. Assume that 𝑢 is the solution of (1.2)–(1.4) such that 𝑢_𝑡𝑡∈𝐿²(𝑄_𝑇) and 𝑢_𝑡∈𝐿^∞(𝑄_𝑇). Then for any positive constant 𝜀>0, there exists a constant 𝐶₂=𝐶₂(𝑚,𝑢₀,𝑇,𝛾₁,𝛾₂)>0, such that 𝐴𝑢^𝑗+1/2−𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁,𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2≤𝑚₀4𝐷²𝑒^𝑗+1/22+𝐶₂𝑒^𝑗+12+𝑒^𝑗2+𝑁⁻⁸+ℎ^3𝑡𝑗+1_𝑡𝑗𝑢_𝑡𝑡²𝑑𝑡+𝑢_𝑡^𝑗+1/22_∞^𝑡𝑗+1_𝑡𝑗𝑢_𝑡²𝑑𝑡. (3.44)

Proof. Firstly, we consider 𝐴𝑢^𝑗+1/2−𝐴𝑢^𝑗+1,𝑢^𝑗=𝛾₂𝑢^𝑗+1/23−𝛾₂4𝑢^𝑗+13+𝑢^𝑗+12𝑢^𝑗+𝑢^𝑗+1𝑢^𝑗2+𝑢^𝑗3+𝛾₁𝑢^𝑗+1/22−𝛾₁3𝑢^𝑗+12+𝑢^𝑗+1𝑢^𝑗+𝑢^𝑗2−𝑢^𝑗+1/2−12𝑢^𝑗+1+𝑢^𝑗≜𝛾₂𝜌^𝑗₁+𝛾₁𝜌^𝑗₂−𝜌^𝑗₃. (3.45) Direct computation gives 𝜌^𝑗₃=12^𝑡𝑗+1/2_𝑡𝑗𝑡−𝑡_𝑗𝑢_𝑡𝑡𝑑𝑡+^𝑡𝑗+1_𝑡𝑗+1/2𝑡_𝑗+1−𝑡𝑢_𝑡𝑡𝑑𝑡≤Δ𝑡³96^𝑡𝑗+1_𝑡𝑗𝑢_𝑡𝑡²𝑑𝑡^1/2,𝜌^𝑗₂=162𝑢^𝑗+1/22−𝑢^𝑗+1+𝑢^𝑗2+162𝑢^𝑗+1/22−𝑢^𝑗+12−𝑢^𝑗2≤162𝑢^𝑗+1/2−𝑢^𝑗+1+𝑢^𝑗2𝑢^𝑗+1/2+𝑢^𝑗+1+𝑢^𝑗+16𝑢^𝑗+1/2−𝑢^𝑗+1𝑢^𝑗+1/2+𝑢^𝑗+1+𝑢^𝑗+1/2−𝑢^𝑗𝑢^𝑗+1/2+𝑢^𝑗≤162𝑢^𝑗+1/2−𝑢^𝑗+1+𝑢^𝑗⋅2𝑢^𝑗+1/2+𝑢^𝑗+1+𝑢^𝑗_∞+16−Δ𝑡2𝑢_𝑡^𝑗+1/2−^𝑡𝑗+1_𝑡𝑗+1/2𝑡_𝑗+1−𝑡𝑢_𝑡𝑡𝑑𝑡𝑢^𝑗+1/2+𝑢^𝑗+1+Δ𝑡2𝑢_𝑡^𝑗+1/2−^𝑡𝑗+1/2_𝑡𝑗𝑡−𝑡_𝑗𝑢_𝑡𝑡𝑑𝑡𝑢^𝑗+1/2+𝑢^𝑗≤𝐶Δ𝑡^{3/2𝑡𝑗+1}_𝑡𝑗𝑢_𝑡𝑡²𝑑𝑡^1/2+Δ𝑡^3/2𝑢12^{𝑡𝑗+1/2∞}⋅^𝑡𝑗+1_𝑡𝑗𝑢_𝑡²𝑑𝑡^1/2,𝜌^𝑗₁=1122𝑢^𝑗+1/23−𝑢^𝑗+1+𝑢^𝑗3+162𝑢^𝑗+1/23−𝑢^𝑗+13−𝑢^𝑗3≤1122𝑢^𝑗+1/2−𝑢^𝑗+1+𝑢^𝑗⋅4𝑢^𝑗+1/22+2𝑢^𝑗+1+𝑢^𝑗𝑢^𝑗+1/2+𝑢^𝑗+1+𝑢^𝑗2+16𝑢^𝑗+1/2−𝑢^𝑗+1𝑢^𝑗+1/22+𝑢^𝑗+1𝑢^𝑗+1/2+𝑢^𝑗+12+𝑢^𝑗+1/2−𝑢^𝑗𝑢^𝑗+1/22+𝑢^𝑗𝑢^𝑗+1/2+𝑢^𝑗2≤𝐶Δ𝑡^{3/2𝑡𝑗+1}_𝑡𝑗𝑢_𝑡𝑡²𝑑𝑡^1/2+𝑢_𝑡^𝑗+1/2_∞⋅^𝑡𝑗+1_𝑡𝑗𝑢_𝑡²𝑑𝑡^1/2. (3.46) Then 𝐴𝑢^𝑗+1/2−𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁,𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2≤3𝐴𝑢4𝜀^𝑗+1/2−𝑃_𝑁𝐴𝑢^𝑗+1/22+𝜀3𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/22+3𝐴𝑢4𝜀^𝑗+1/2−𝐴𝑢^𝑗+1,𝑢^𝑗2+𝜀3𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/22+3𝐴𝑢4𝜀^𝑗+1,𝑢^𝑗−𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁²+𝜀3𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/22≤𝜀𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/22+3𝐴𝑢4𝜀^𝑗+1/2−𝑃_𝑁𝐴𝑢^𝑗+1/22+3𝐴𝑢4𝜀^𝑗+1/2−𝐴𝑢^𝑗+1,𝑢^𝑗2+3𝐴𝑢4𝜀^𝑗+1,𝑢^𝑗−𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁². (3.47) Direct computation yields 𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/22≤𝑀²₀+𝑀²₁𝐷²𝑒^𝑗+1/22+𝐶𝑒^𝑗+1/22,𝐴𝑢^𝑗+1/2−𝑃_𝑁𝐴𝑢^𝑗+1/22≤𝐶𝑁⁻⁸,𝐴𝑢^𝑗+1/2−𝐴𝑢^𝑗+1,𝑢^𝑗2≤𝐶𝜌^𝑗₁²+𝜌^𝑗₂²+𝜌^𝑗₃²,𝐴𝑢^𝑗+1,𝑢^𝑗−𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁²≤𝐺^𝑗₁²_∞+𝐺^𝑗₂²_∞𝑒^𝑗+12+𝑒^𝑗2+𝐶𝑁⁻⁸, (3.48) where 𝐺^𝑗₁=𝛾₂8𝑢^𝑗+1+𝑈_𝑁^𝑗+12+𝑈_𝑁^𝑗+1+𝑈^𝑗_𝑁²+𝑢^𝑗+1+𝑈^𝑗_𝑁²+𝛾₁3𝑈_𝑁^𝑗+1+𝑢^𝑗+1+𝑈^𝑗_𝑁−12,𝐺^𝑗₂=𝛾₂8𝑢^𝑗+𝑈^𝑗_𝑁²+𝑢^𝑗+𝑢^𝑗+12+𝑢^𝑗+1+𝑈^𝑗_𝑁²+𝛾₁3𝑢^𝑗+1+𝑢^𝑗+𝑈^𝑗_𝑁−12. (3.49) Applying Lemma 2.3 and Theorem 3.7, we obtain that ‖𝐺₁‖_∞≤𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂) and ‖𝐺₂‖_∞≤𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂). Taking 𝜀=𝑚₀/4(𝑀²₀+𝑀²₁) in (3.47),we have 𝐴𝑢^𝑗+1/2−𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁,𝐷𝑚^𝑗+1/2𝐷𝑒^𝑗+1/2≤𝑚₀4𝐷²𝑒^𝑗+1/22+𝐶₂𝑒^𝑗+12+𝑒^𝑗2+𝑁⁻⁸+Δ𝑡^3𝑡𝑗+1_𝑡𝑗𝑢_𝑡𝑡²𝑑𝑡+𝑢_𝑡^𝑗+1/22_∞^𝑡𝑗+1_𝑡𝑗𝑢_𝑡²𝑑𝑡, (3.50) where 𝐶₂=𝐶₂(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂)>0 is a constant.

Lemma 3.6. Assume that 𝑢 is the solution of (1.2)–(1.4). Then there exists a positive constant 𝐶=𝐶(𝑢₀,𝑚,𝛾₁,𝛾₂)>0 such that 𝐷²𝑢^𝑗+1/2−𝐴𝑢^𝑗+1/2,𝐷𝑚𝑈_𝑁^𝑗+1/2−𝑚𝑢^𝑗+1/2𝐷𝑒^𝑗+1/2≤𝑚₀4𝐷²𝑒^𝑗+1/22+𝐶𝑒^𝑗+1/22+𝑁⁻⁸+ℎ^3𝑡𝑗+1_𝑡𝑗𝑢_𝑡𝑡𝑑𝑡. (3.51)

Proof. By (2.9), we have 𝐷³𝑢^𝑗+1/2−𝐷𝐴𝑢^𝑗+1/2_∞≤𝐷³𝑢^𝑗+1/2_∞+𝐴^′𝑢^𝑗+1/2𝐷𝑢^𝑗+1/2_∞≤𝐷³𝑢^𝑗+1/2_∞+𝐴^′𝑢^𝑗+1/2_∞𝐷𝑢^𝑗+1/2_∞≤𝐶. (3.52) In the other hand, 𝑚𝑈_𝑁^𝑗+1/2−𝑚𝑢^𝑗+1/22≤𝑀²₁𝑈_𝑁^𝑗+1/2−𝑢^𝑗+1/22≤𝐶𝑈_𝑁^𝑗+1/2−𝑢^𝑗+𝑢^𝑗+12²+𝑢^𝑗+𝑢^𝑗+12−𝑢^𝑗+1/22≤𝐶𝑒^𝑗+1/22+𝜂^𝑗+1/22+ℎ³96^𝑡𝑗+1_𝑡𝑗𝑢_𝑡𝑡𝑑𝑡. (3.53) By Young inequality, we obtain 𝐷²𝑢^𝑗+1/2−𝐴𝑢^𝑗+1/2,𝐷𝑚𝑈_𝑁^𝑗+1/2−𝑚𝑢^𝑗+1/2𝐷𝑒^𝑗+1/2=−𝐷³𝑢^𝑗+1/2−𝐷𝐴𝑢^𝑗+1/2,𝑚𝑈_𝑁^𝑗+1/2−𝑚𝑢^𝑗+1/2𝐷𝑒^𝑗+1/2≤𝐷³𝑢^𝑗+1/2−𝐷𝐴𝑢^𝑗+1/2_∞𝑚𝑈_𝑁^𝑗+1/2−𝑚𝑢^𝑗+1/2𝐷𝑒^𝑗+1/2≤𝜀𝐷𝑒^𝑗+1/22+𝐶_𝜀𝑚𝑈_𝑁^𝑗+1/2−𝑚𝑢^𝑗+1/22≤𝜀𝐷²𝑒^𝑗+1/2𝑒^𝑗+1/2+𝐶_𝜀𝑚𝑈_𝑁^𝑗+1/2−𝑚𝑢^𝑗+1/22≤𝜀𝐷²𝑒^𝑗+1/22+𝐶_𝜀𝑒^𝑗+1/22+𝑁⁻⁸+ℎ^3𝑡𝑗+1_𝑡𝑗𝑢_𝑡𝑡𝑑𝑡. (3.54) Choosing 𝜀=𝑚₀/4 in the previous inequality leads to (3.51).

Finally, we obtain the main theorem of this paper.

Theorem 3.7. Assume that 𝑢(𝑥,𝑡) is the solution of (1.2)–(1.4) and satisfies that 𝐷⁴𝑢∈𝐿^∞𝑄_𝑇,𝑢_𝑡∈𝐿^∞𝑄_𝑇,𝐷²𝑢_𝑡𝑡∈𝐿²𝑄_𝑇,𝑢_𝑡𝑡𝑡∈𝐿²𝑄_𝑇. (3.55)𝑈^𝑗_𝑁∈𝑆_𝑁(𝑗=1,2,…,Λ) is the solution of the full-discretization (3.1) and (3.2). If ℎ is sufficiently small, there exists a positive constant 𝐶 such that, for any 𝑗=0,1,2,…,Λ, 𝑒^𝑗+1=𝑃_𝑁𝑢𝑡_𝑗+1−𝑈_𝑁^𝑗+1≤𝐶𝑁⁻²+𝑒⁰+𝐵ℎ², (3.56) where 𝐵=∫^𝑡𝑗+1₀(‖𝐷²𝑢_𝑡𝑡‖²+‖𝑢_𝑡𝑡‖²+‖𝑢_𝑡𝑡𝑡‖²+max_{0≤𝑙≤𝑗}{‖𝑢_𝑡^𝑙+1/2‖²_∞}⋅‖𝑢_𝑡‖²)𝑑𝑡.

Proof. By (3.27), (3.36), (3.44), and (3.51), we have 𝑒^𝑗+12≤𝑒^𝑗2+ℎ⁴320^𝑡𝑗+1_𝑡𝑗𝑢_𝑡𝑡𝑡²𝑑𝑡+ℎ𝑒^𝑗+1/22+2ℎ−𝐷²𝑢^𝑗+1/2−𝐷²𝑈_𝑁^𝑗+1+𝐷²𝑈^𝑗_𝑁2,𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2+𝐴𝑢^𝑗+1/2−𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁,𝐷𝑚𝑈_𝑁^𝑗+1/2𝐷𝑒^𝑗+1/2+𝐷²𝑢^𝑗+1/2−𝐴𝑢^𝑗+1/2,𝐷𝑚𝑈_𝑁^𝑗+1/2−𝑚𝑢^𝑗+1/2𝐷𝑒^𝑗+1/2≤𝑒^𝑗2+ℎ𝐶₁𝑁⁻⁴+𝑒^𝑗+12+𝑒^𝑗2+𝐶₂ℎ^4𝑡j+1_𝑡𝑗𝐷²𝑢_𝑡𝑡²+𝑢_𝑡𝑡²+𝑢_𝑡𝑡𝑡²+𝑢_𝑡^𝑗+1/22_∞𝑢_𝑡²𝑑𝑡, (3.57) where 𝐶₁=𝐶₁(𝑚,𝑢₀,𝑇,𝛾₁,𝛾₂)>0 and 𝐶₂=𝐶₂(𝑚,𝑢₀,𝑇,𝛾₁,𝛾₂)>0 are constants. For sufficiently small ℎ such that 𝐶₁ℎ≤1/2, denoting 𝐶=2(𝐶₁+𝐶₂), we obtain 𝑒^𝑗+12≤1+𝐶ℎ𝑒^𝑗2+𝐶ℎ𝑁⁻⁴+ℎ⁴𝐵^𝑗, (3.58) where 𝐵^𝑗=^𝑡𝑗+1_𝑡𝑗𝐷²𝑢_𝑡𝑡²+𝑢_𝑡𝑡²+𝑢_𝑡𝑡𝑡²+𝑢_𝑡^𝑗+1/22_∞𝑢_𝑡²𝑑𝑡. (3.59) By the Gronwall inequality of the discrete form, we obtain 𝑒^𝑗+12≤1+𝐶ℎ^𝑗+1𝑒⁰²+𝐶^𝑗_𝑙=01+𝐶ℎ^𝑙ℎ𝑡𝑁⁻⁴+ℎ⁴𝐵^𝑙≤exp𝐶𝑗+1ℎ𝑒⁰²+𝐶^𝑗_𝑙=0ℎ𝑁⁻⁴+ℎ⁴𝐵^𝑙≤exp𝐶𝑗+1ℎ𝑒⁰²+𝐶𝑗ℎ𝑁⁻⁴+ℎ^4𝑗_𝑙=0𝐵^𝑙. (3.60) Direct computation gives ^𝑗_𝑙=0𝐵^𝑙≤^𝑡𝑗+1₀𝐷²𝑢_𝑡𝑡²+𝑢_𝑡𝑡²+𝑢_𝑡𝑡𝑡²+max_{0≤𝑙≤𝑗}𝑢_𝑡^𝑙+1/22_∞⋅𝑢_𝑡²𝑑𝑡. (3.61) Then we complete the conclusion (3.56).

Furthermore, we get the following theorem.

Theorem 3.8. Assume 𝑢 is the solution of (1.2)–(1.4) and satisfies that 𝐷⁴𝑢∈𝐿^∞𝑄_𝑇,𝑢_𝑡∈𝐿^∞𝑄_𝑇,𝐷²𝑢_𝑡𝑡∈𝐿²𝑄_𝑇,𝑢_𝑡𝑡𝑡∈𝐿²𝑄_𝑇. (3.62)𝑈^𝑗_𝑁∈𝑆_𝑁(𝑗=1,2,…,Λ) is the solution of the full-discretization scheme (3.1)-(3.2), and 𝑈⁰ satisfies ‖𝑒⁰‖=‖𝑃_𝑁𝑢⁰−𝑈⁰‖≤𝐶𝑁⁻⁴. If ℎ is sufficiently small, there exists a constant 𝐶=𝐶(𝑢₀,𝑚,𝑇,𝛾₁,𝛾₂)>0 such that 𝑢𝑥,𝑡_𝑗−𝑈^𝑗_𝑁≤𝐶𝑁⁻²+ℎ²,𝑗=1,2,…,Λ. (3.63)

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) where 𝑚₀>0 is a constant. The full-discretization spectral method of (2.2) and (2.3) reads: find 𝑈^𝑗_𝑁=∑^𝑁_𝑙=0𝛼_𝑗𝑙cos𝑙𝜋𝑥(𝑗=1,2,…,Λ) such that 𝑈_𝑁^𝑗+1−𝑈^𝑗_𝑁ℎ+𝐷,𝑣²𝑈_𝑁^𝑗+1+𝐷²𝑈^𝑗_𝑁2−𝑃_𝑁𝐴𝑈_𝑁^𝑗+1,𝑈^𝑗_𝑁,𝐷𝑚𝐷²𝑈_𝑁^𝑗+1+𝐷²𝑈^𝑗_𝑁2𝑈𝐷𝑣=0,_0𝑁,𝑣−𝑢₀,𝑣=0. (4.2) In our computations we fix 𝑁=32 and choose five different time-step sizes ℎ_𝑘(𝑘=1,2,…,5). Let Λ_𝑘 be the integer with ℎ_𝑘Λ_𝑘=𝑇. Since we have no exact solution of (2.2) and (2.3), we take 𝑁=32 and ℎ₀=0.15625×10⁻⁴ to compute an approximating solution 𝑈^Λ0_𝑁 with ℎ₀Λ₀=𝑇 and regard this as an exact solution. we also choose five different time-step sizes ℎ_𝑘(𝑘=1,2,…,5) with ℎ_𝑘Λ_𝑘=𝑇 to obtain five approximating solutions 𝑈^Λ𝑘_𝑁(𝑘=1,2,…,5) and compute the error estimation. Define an error function: err𝑇,ℎ_𝑘=¹₀𝑈^Λ𝑘_𝑁−𝑈^Λ0_𝑁²𝑑𝑥^1/2. (4.3) This function characterizes the estimations with respect to time-step size.

4.1. Example 1

Take 𝑚₀=1 and 𝑇=0.1. We also take two different initial functions 𝑢¹₀(𝑥)=𝑥⁵(1−𝑥)⁵ and 𝑢²₀=𝑥⁵(1−𝑥)⁵𝑒^𝑥 to carry out numerical computations. Figure 1 shows the development of the solutions for time 𝑡 from 𝑡=0 to 𝑡=0.1 with fixed step-size ℎ₀=0.15625×10⁻⁴.

(a)

(b)

Figure 1 

The development of the solution of the full-discrete scheme when initial value is 𝑢0(𝑥)=𝑥5(1−𝑥)5 (a) and 𝑢0(𝑥)=𝑥5(1−𝑥)5𝑒𝑥 (b).

We also choose five different time-step sizes ℎ_𝑘 to carry out numerical computations and apply the error function in (4.3) to illustrate the estimation and convergence order in time variable 𝑡, see Table 1.

Table 1 

The error and the convergence order.

4.2. Example 2

Take 𝑚₀=0.05 and 𝑇=0.1. We also take two different initial functions 𝑢¹₀(𝑥)=𝑥⁵(1−𝑥)⁵ and 𝑢²₀=𝑥⁵(1−𝑥)⁵𝑒^𝑥 to carry out numerical computations. Figure 2 shows the development of the solutions for time 𝑡 from 𝑡=0 to 𝑡=0.1 with fixed step-sizes ℎ₀=0.15625×10⁻⁴.

(a)

(b)

Figure 2 

The development of the solution of the full-discrete scheme with initial value 𝑢0(𝑥)=𝑥5(1−𝑥)5 (a) and 𝑢0(𝑥)=𝑥5(1−𝑥)5𝑒𝑥 (b).

We also choose five different time-step sizes ℎ_𝑘 to carry out numerical computations and apply the error function in (4.3) to illustrate the estimation and convergence order in time variable 𝑡, see Table 2.

Table 2 

The error and the convergence order.

4.3. Example 3

Take 𝑚₀=0.005 and 𝑇=0.1. We also take two different initial functions 𝑢¹₀(𝑥)=𝑥⁵(1−𝑥)⁵ and 𝑢²₀=𝑥⁵(1−𝑥)⁵𝑒^𝑥 to carry out numerical computations. Figure 3 shows the development of the solutions for time 𝑡 from 𝑡=0 to 𝑡=0.1 with fixed step-sizes ℎ₀=0.15625×10⁻⁴.

(a)

(b)

Figure 3 

The development of the solution of the full-discrete scheme with initial value 𝑢0(𝑥)=𝑥5(1−𝑥)5 (a) and 𝑢0(𝑥)=𝑥5(1−𝑥)5𝑒𝑥 (b).

We also choose five different time-step sizes ℎ_𝑘 to carry out numerical computations and apply the error function in (4.3) to illustrate the estimation and convergence order in time variable 𝑡, see Table 3.

Table 3 

The error and the convergence order.

Acknowledgment

This work is supported by NSFC no. 11071102.