A Godunov-Mixed Finite Element Method on Changing Meshes for the Nonlinear Sobolev Equations

Abstract

A Godunov-mixed finite element method on changing meshes is presented to simulate the nonlinear Sobolev equations. The convection term of the nonlinear Sobolev equations is approximated by a Godunov-type procedure and the diffusion term by an expanded mixed finite element method. The method can simultaneously approximate the scalar unknown and the vector flux effectively, reducing the continuity of the finite element space. Almost optimal error estimates in 𝐿²-norm under very general changes in the mesh can be obtained. Finally, a numerical experiment is given to illustrate the efficiency of the method.

1. Introduction

We consider the following nonlinear Sobolev equations: 𝑢_𝑡−∇⋅𝑎(𝑥,𝑡,𝑢)∇𝑢_𝑡()()()(],𝑢()[],𝑢()+𝑏𝑥,𝑡,𝑢∇𝑢+𝐜𝑥,𝑡,𝑢⋅∇𝑢=𝑓𝑥,𝑡,𝑢,𝑥∈Ω,𝑡∈0,𝑇𝑥,𝑡=0,𝑥∈𝜕Ω,𝑡∈0,𝑇𝑥,0=𝑢₀(𝑥),𝑥∈Ω, (1.1) where Ω is a bounded subset of 𝑅^𝑛 (𝑛≤3) with smooth boundary 𝜕Ω, 𝑢₀(𝑥) and 𝑓(𝑥,𝑡,𝑢) are known functions, and the coefficients 𝑎(𝑥,𝑡,𝑢),𝑏(𝑥,𝑡,𝑢),𝐜(𝑥,𝑡,𝑢)=(𝑐₁(𝑥,𝑡,𝑢),𝑐₂(𝑥,𝑡,𝑢))^𝚃 satisfy the following condition: 0<𝑎₀≤𝑎(𝑥,𝑡,𝑢)≤𝑎₁,𝜕𝑎(𝑥,𝑡,𝑢)𝜕𝑡≤𝐾₁,0<𝑏₀≤𝑏(𝑥,𝑡,𝑢)≤𝑏₁,|𝐜(𝑥,𝑡,𝑢)|=𝑐²₁(𝑥,𝑡,𝑢)+𝑐²₂(𝑥,𝑡,𝑢)≤𝑐₀,𝜕𝐜(𝑥,𝑡,𝑢)𝜕𝑢≤𝐾₂,𝑎()()()(),()𝑥,𝑡,𝑢,𝑏𝑥,𝑡,𝑢,𝐜𝑥,𝑡,𝑢,𝑓𝑥,𝑡,𝑢𝜕𝐜𝑥,𝑡,𝑢𝜕𝑢areLipschitzcontinuousabout𝑢, (1.2) where 𝑎₀, 𝑎₁, 𝑏₀, 𝑏₁, 𝑐₀, 𝐾₁, and 𝐾₂ are positive constants. We assume that 𝑢(𝑥,𝑡) satisfy the smooth condition in the following analysis.

For time-changing localized phenomena, such as sharp fronts and layers, the finite element method on changing meshes [1–3] is advantageous over fixed finite element method. The reason is that the former treats the problem with the finite element method on space domain by using different meshes and different basic functions at different time levels so that it has the capability of self-adaptive local grid modification (refinement or unrefinement) to efficiently capture propagating fronts or moving layers. The work [4] had combined this method with mixed finite element method to study parabolic problems. In [5], an upwind-mixed method on changing meshes was considered for two-phase miscible flow in porous media.

Sobolev equations have important applications in many mathematical and physical problems, such as the percolation theory when the fluid flows through the cracks [6], the transfer problem of the moisture in the soil [7], and the heat conduction problem in different materials [8]. So there exists great and actual significance to research Sobolev equations. Many works had researched on numerical treatments for Sobolev equations. More attentions were paid for treating a damping term ∇⋅[𝑎∇𝑢_𝑡], which is a distinct character of Sobolev equations different from parabolic equation. For example, time stepping Galerkin method was presented for nonlinear Sobolev equations in [9, 10]. In [11, 12], nonlinear Sobolev equations with convection term were researched by using finite difference streamline-diffusion method and discontinuous Galerkin method, respectively. Two new least-squares mixed finite element procedures were formulated for solving convection-dominated Sobolev equations in [13].

Methods which combined Godunov-type schemes for advection with mixed finite elements for diffusion were introduced in [14] and had been applied to flow problems in reservoir engineering, contaminant transport, and computational fluid dynamics. Applications of these types of methods to single and two-phase flow in oil reservoirs were discussed in [15, 16]; application to the Navier-Stokes equations was given in [17]. These methods had proven useful for advective flow problems because they combined element-by-element conservation of mass with numerical stability and minimal numerical diffusion. Dawson [18] researched advective flow problems in one space dimension by high-order Godunov-mixed method. In 1993, he expanded this method to multidimensions [19] and presented three variations. In these methods, advection was approximated by a Godunov-type procedure, and diffusion was approximated by a low-order-mixed finite element method.

The object of this paper is to present a Godunov-mixed finite element method on changing meshes for the nonlinear Sobolev equations. The convection term 𝐜⋅∇𝑢 of the nonlinear Sobolev equations is approximated by a Godunov-type procedure and the diffusion term by an expanded mixed finite element method. This method can simultaneously approximate the scalar unknown and the vector flux effectively, reducing the continuity of the finite element space. We describe this method in Section 2. In Section 3, we introduce three projections and a lemma. We derive almost optimal error estimates in 𝐿²-norm under very general changes in the mesh in Section 4. In Section 5, we present results of numerical experiment, which confirm our theoretical results.

Throughout the analysis, the symbol 𝐾 will denote a generic constant, which is independent of mesh parameters Δ𝑡 and ℎ and not necessarily the same at different occurrences.

2. The Godunov-Mixed Method on Changing Meshes

At first we give some notation and basic assumptions. The usual Sobolev spaces and norms are adopted on Ω. The inner product on 𝐿²(Ω) is denoted by (𝑓,𝑔)=∫_Ω𝑓𝑔𝑑𝑥. Define the following spaces and norms: 𝐻^𝑚(div,Ω)=𝐟=𝑓_𝑥,𝑓_𝑦;𝑓_𝑥,𝑓_𝑦,∇⋅𝐟∈𝐻^𝑚(Ω‖𝐟‖),𝑚≥0,_{2𝐻𝑚(div)}=𝑓_𝑥²_𝑚+𝑓_𝑦²_𝑚+‖∇⋅𝐟‖²_𝑚,𝐻(div,Ω)=𝐻⁰()𝐿div,Ω,𝑊=²(Ω),{()}.𝜑≡constantonΩ𝑉=𝐯∈𝐻div;Ω∣𝐯⋅𝐧=0on𝜕Ω,𝐧istheunitoutwardnormto𝜕Ω (2.1)

Let Δ𝑡^𝑛>0  (𝑛=1,2,…,𝑁^∗) denote different time steps, 𝑡^𝑛=∑^𝑛_𝑘=1Δ𝑡^𝑘, 𝑇=∑^𝑁∗_𝑛=1Δ𝑡^𝑛, Δ𝑡=max_𝑛Δ𝑡^𝑛. Assume that the time steps Δ𝑡^𝑛 do not change too rapidly; that is, we assume there exist positive constants 𝑡_∗ and 𝑡^∗ which are independent of 𝑛 and Δ𝑡 such that 𝑡_∗≤Δ𝑡^𝑛Δ𝑡^𝑛−1≤𝑡^∗. (2.2) For a given function 𝑔(𝑥,𝑡), let 𝑔^𝑛=𝑔(𝑥,𝑡^𝑛).

Assume Ω=(0,1)×(0,1). At each time level 𝑡^𝑛, we construct a quasiuniform rectangular partition 𝐾^𝑛_ℎ={𝑒^𝑛_𝑖} of Ω: 𝛿^𝑛_𝑥∶0=𝑥^𝑛_1/2<𝑥^𝑛_3/2<⋯<𝑥^𝑛_𝑖+1/2<⋯<𝑥^𝑛_𝐼^𝑛_+1/2𝛿=1,_𝑛𝑦∶0=𝑦^𝑛_1/2<𝑦^𝑛_3/2<⋯<𝑦^𝑛_𝑗+1/2<⋯<𝑦^𝑛_𝐽𝑛+1/2=1. (2.3) Let ℎ^𝑛_𝑖,𝑥=𝑥^𝑛_𝑖+1/2−𝑥^𝑛_𝑖−1/2, ℎ^𝑛_𝑖,𝑦=𝑦^𝑛_𝑗+1/2−𝑦^𝑛_𝑗−1/2, ℎ^𝑛=max_𝑖,𝑗{ℎ^𝑛_𝑖,𝑥,ℎ^𝑛_𝑖,𝑦}, and ℎ=max_𝑛ℎ^𝑛. And 𝑥^𝑛_𝑖=(𝑥^𝑛_𝑖+1/2+𝑥^𝑛_𝑖−1/2)/2 is the midpoint of [𝑥^𝑛_𝑖−1/2,𝑥^𝑛_𝑖+1/2], ℎ^𝑛_{𝑖+1/2,𝑥}=𝑥^𝑛_𝑖+1−𝑥^𝑛_𝑖, Δ𝑡^𝑛=𝑂(ℎ^𝑛). Let 𝑦^𝑛_𝑗,ℎ^𝑛_{𝑗+1/2,𝑦} be defined analogously.

Let 𝑀^𝑘₋₁𝛿^𝑛_𝑥=𝑓∶𝑓∣_𝐵^𝑛_𝑖,𝑥∈𝑃^𝑘𝐵^𝑛_𝑖,𝑥,𝑖=1,…,𝐼^𝑛,𝑀^𝑘₀𝛿^𝑛_𝑥=𝑀^𝑘₋₁𝛿^𝑛_𝑥∩𝐶⁰([0,1]), (2.4) where 𝐵^𝑛_𝑖,𝑥=[𝑥^𝑛_𝑖−1/2,𝑥^𝑛_𝑖+1/2] and 𝑃^𝑘(𝐵^𝑛_𝑖,𝑥) is the set of all polynomials of degree less than or equal to 𝑘 defined on 𝐵^𝑛_𝑖,𝑥. Similar definitions are given to 𝑀⁰₋₁(𝛿^𝑛_𝑦) and 𝑀¹₀(𝛿^𝑛_𝑦).

Then the lowest order Raviart-Thomas spaces 𝑊^𝑛_ℎ and 𝑉^𝑛_ℎ are given by 𝑊^𝑛_ℎ=𝑀⁰₋₁𝛿^𝑛_𝑥⊗𝑀⁰₋₁𝛿^𝑛_𝑦,𝑉^𝑛_ℎ=𝑉^𝑛_ℎ,𝑥×𝑉^𝑛_ℎ,𝑦,𝑉^𝑛_ℎ,𝑥=𝑀¹₀𝛿^𝑛_𝑥⊗𝑀⁰₋₁𝛿^𝑛_𝑦,𝑉^𝑛_ℎ,𝑦=𝑀⁰₋₁𝛿^𝑛_𝑥⊗𝑀¹₀𝛿^𝑛_𝑦. (2.5) That is, the space 𝑊^𝑛_ℎ is the space of functions which are constant on each element 𝑒^𝑛_𝑖∈𝐾^𝑛_ℎ, and 𝑉^𝑛_ℎ is the space of vector valued functions whose components are continuous and linear on each element 𝑒^𝑛_𝑖∈𝐾^𝑛_ℎ. The degrees of freedom of a function 𝐯^𝑛∈𝑉^𝑛_ℎ correspond to the values of 𝐯^𝑛⋅𝜸 at the midpoints of the sides of 𝑒^𝑛_𝑖, where 𝜸 is the unit outward normal to 𝜕𝑒^𝑛_𝑖. It is easy to see that 𝑊^𝑛_ℎ×𝑉^𝑛_ℎ⊂𝑊×𝑉 and div𝑉^𝑛_ℎ=𝑊^𝑛_ℎ.

By introducing variables 𝐳=−∇𝑢,  𝐳=𝑏𝐳+𝑎𝐳_𝑡,  𝐠=𝐜𝑢=(𝑐₁𝑢,𝑐₂𝑢)^𝚃=(𝑔₁,𝑔₂)^𝚃, and 𝐜(𝑥,𝑡,𝑢)=(𝜕𝑐₁/𝜕𝑢,𝜕𝑐₂/𝜕𝑢)^𝚃, we modify the first equation in (1.1) as 𝑢_𝑡+∇⋅𝐳+(∇⋅𝐠+𝑢𝐜(𝑢)⋅𝐳)=𝑓(𝑢). (2.6)

Here we are using the so-called “expanded” mixed finite element method, proposed by Arbogast et al. [20], which gives a gradient approximation 𝐳 as well as an approximation to the diffusion term 𝐳.

The weak form of (2.6) is 𝑢_𝑡,𝑤+(∇⋅𝐳,𝑤)+(∇⋅𝐠,𝑤)+(𝑢𝐜(𝑢)⋅𝐳,𝑤)=𝑓(𝑢),𝑤, (2.7) so that 𝑢_𝑡,𝑤+(∇⋅𝐳,𝑤)+(∇⋅𝐠,𝑤)+(𝑢𝐜(𝑢)⋅)=𝑓(𝑢)(𝐳,𝑤,𝑤,∀𝑤∈𝑊,)=()()=𝑏(𝑢)𝐳,𝐯𝑢,∇⋅𝐯,∀𝐯∈𝑉,𝐳,𝐯𝐳+𝑎(𝑢)𝐳_𝑡,𝐯,∀𝐯∈𝑉. (2.8)

The Godunov-mixed method on changing meshes is presented as follows: at each time level 𝑛, find 𝑈^𝑛∈𝑊^𝑛_ℎ, 𝐙^𝑛∈𝑉^𝑛_ℎ such that 𝑈^𝑛−𝑅^𝑛𝑈^𝑛−1Δ𝑡^𝑛,𝑤^𝑛+∇⋅𝐙^𝑛,𝑤^𝑛+𝑈^𝑛𝐜𝑈^𝑛⋅𝑍^𝑛,𝑤^𝑛=𝑓𝑈^𝑛,𝑤^𝑛−∇⋅𝐆^𝑛−1,𝑤^𝑛,∀𝑤^𝑛∈𝑊^𝑛_ℎ,𝑅^𝑛𝑈^𝑛−1−𝑈^𝑛−1,𝑤^𝑛=0,∀𝑤^𝑛∈𝑊^𝑛_ℎ,𝐙^𝑛,𝐯^𝐧=𝑈^𝑛,∇⋅𝐯^𝐧,∀𝐯^𝐧∈𝑉^𝑛_ℎ,𝐙^𝑛,𝐯^𝐧=𝑏𝑈^𝑛𝐙^𝑛+𝑎𝑈^𝑛𝐙^𝑛−𝑅^𝑛𝐙^𝑛−1Δ𝑡^𝑛,𝐯^𝐧,∀𝐯^𝐧∈𝑉^𝑛_ℎ,𝑎𝑈^𝑛𝑅𝐙^{𝑛∼𝑛−1}−𝐙^∼𝑛−1,𝐯^𝐧=0,∀𝐯^𝐧∈𝑉^𝑛_ℎ. (2.9) When different finite element spaces are used at time levels 𝑡^𝑛 and 𝑡^𝑛−1, the second and fifth equations of (2.9) give the 𝐿²-projection {𝑅^𝑛𝑈^𝑛−1,𝑅𝐙^{𝑛∼𝑛−1}} of the previous approximate solution {𝑈^𝑛−1,𝐙^∼𝑛−1} into the current finite element space 𝑊^𝑛_ℎ×𝑉^𝑛_ℎ. This projection is used in the first and third equations of (2.9) as initial value to calculate {𝑈^𝑛,𝐙^∼𝑛}. If the finite element spaces are the same at time levels 𝑡^𝑛 and 𝑡^𝑛−1, we know that 𝑅^𝑛𝑈^𝑛−1=𝑈^𝑛−1, 𝑅𝐙^{𝑛∼𝑛−1}=𝐙^∼𝑛−1.

Let 𝐆^𝑛 in (2.9) be constructed by the Godunov method. We adopt the Godunov flux [19] 𝐻_𝑤(𝑤^𝐿,𝑤^𝑅) for a given flux function 𝑤(𝑠): 𝐻_𝑤𝑤^𝐿,𝑤^𝑅=min_{𝑤𝐿≤𝑠≤𝑤𝑅}𝑤(𝑠),if𝑤^𝐿≤𝑤^𝑅,max_{𝑤𝑅≤𝑠≤𝑤𝐿}𝑤(𝑠),otherwise, (2.10) where 𝑤^𝐿,𝑤^𝑅 are the left and right states, respectively. It is easily seen that 𝐻_𝑤 is Lipschitz in its arguments if 𝑤 is Lipschitz in 𝑠 and 𝐻_𝑤 is consistent, that is, 𝐻_𝑤(𝑠,𝑠)=𝑤(𝑠).

Let 𝐆=(𝑔₁,𝑔₂), 𝑔_𝑖=𝑐_𝑖𝑢, 𝐆^𝑛−1_{𝑖+1/2,𝑗+1/2}=(𝑔^𝑛−1_{1,𝑖+1/2,𝑗},𝑔^𝑛−1_{2,𝑖,𝑗+1/2}). Now we give the calculation of 𝐆 by 𝐻_𝑤 by the following three steps.

First Step
We construct a piecewise linear function 𝑅𝑈^𝑛−1 on each element 𝑒_𝑖^𝑛−1: 𝑅𝑈^𝑛−1∣_𝑒𝑖^𝑛−1=𝑈^𝑛−1_𝑖,𝑗+𝑥−𝑥_𝑖𝛿_𝑥𝑈^𝑛−1_𝑖,𝑗+𝑦−𝑦_𝑗𝛿_𝑦𝑈^𝑛−1_𝑖,𝑗, (2.11) where 𝛿_𝑥𝑈^𝑛−1_𝑖,𝑗 and 𝛿_𝑦𝑈^𝑛−1_𝑖,𝑗 are 𝑥,𝑦 slopes. The 𝑥 slope calculation can be performed as 𝛿_𝑥𝑈^𝑛−1_𝑖,𝑗=▵^𝑥₊𝑈^𝑛−1_𝑖,𝑗,if▵^𝑥₊𝑈^𝑛−1_𝑖,𝑗≤▵^𝑥₋𝑈^𝑛−1_𝑖,𝑗,▵^𝑥₋𝑈^𝑛−1_𝑖,𝑗,otherwise, (2.12) where ▵^𝑥₊ and ▵^𝑥₋ are forward and backward difference operators in the 𝑥 direction, respectively: ▵^𝑥₊𝑈^𝑛−1_𝑖,𝑗=𝑈^𝑛−1_𝑖+1,𝑗−𝑈^𝑛−1_𝑖,𝑗ℎ^𝑛−1_{𝑖+1/2,𝑥},▵^𝑥₋𝑈^𝑛−1_𝑖,𝑗=𝑈^𝑛−1_𝑖,𝑗−𝑈^𝑛−1_{𝑖−1,𝑗}ℎ^𝑛−1_{𝑖−1/2,𝑥}, (2.13) and 𝛿_𝑥 satisfy ℎ^𝑛−1_𝑖,𝑥2𝛿_𝑥𝑈^𝑛−1_𝑖,𝑗≤𝐶_{𝑗−𝐿≤𝑗≤𝑗+𝐿𝑖−𝑅≤𝑖≤𝑖+𝑅}𝑈^𝑛−1_𝑖,𝑗,𝛿_𝑥𝑢^𝑛−1_𝑖,𝑗=𝑢_𝑥𝑥_𝑖,𝑦_𝑗,𝑡^𝑛−1+𝑂(ℎ), (2.14) where 𝐿 and 𝑅 are positive integers independent of ℎ. The 𝑦 slope is defined analogously.

Second Step
By the Taylor expansion, for 𝑢 with smooth second derivatives, we have 𝑢𝑥_𝑖+1/2,𝑦_𝑗,𝑡^𝑛−1/2ℎ=𝑢+^{𝑛𝑖,𝑥}2𝑢_𝑥+Δ𝑡^𝑛2𝑢_𝑡+𝑂ℎ²+Δ𝑡². (2.15) By (2.6), we see that 𝑢_𝑡=𝑓−𝑢𝐜(𝑢)⋅𝑧−𝑔_1𝑢𝑢_𝑥−𝑔_2𝑦𝑢−∇⋅𝑧,_{𝑛−1/2𝑖+1/2,𝑦}ℎ=𝑢+^{𝑛𝑖,𝑥}2−Δ𝑡^𝑛2𝑔_1𝑢𝑢_𝑥−Δ𝑡^𝑛2𝑔_2𝑦+𝜎, (2.16)(2.17) where 𝜎=Δ𝑡^𝑛2𝑓−𝑢𝐜(𝑢)⋅𝑧−∇⋅𝑧+𝑂ℎ²+Δ𝑡²=𝑂ℎ²+Δ𝑡. (2.18)
Based on (2.17) and a similar expansion about (𝑥_𝑖+1,𝑦_𝑗,𝑡^𝑛−1), we define 𝑈^𝐿,𝑈^𝑅: 𝑈^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}≜𝑈^𝑛−1_𝑖,𝑗+ℎ^𝑛−1_𝑖,𝑥2−Δ𝑡^𝑛−12𝑔_1𝑢𝑈^𝑛−1_𝑖,𝑗𝛿_𝑥𝑈^𝑛−1_𝑖,𝑗−Δ𝑡^𝑛−12ℎ^𝑛−1_𝑗,𝑦Γ^𝑛−1_{𝑖,𝑗+1/2}−Γ^𝑛−1_{𝑖,𝑗−1/2},𝑈^{𝑅,𝑛−1}_{𝑖+1/2,𝑗}≜𝑈^𝑛−1_𝑖+1,𝑗+ℎ^𝑛−1_𝑖,𝑥2−Δ𝑡^𝑛−12𝑔_1𝑢𝑈^𝑛−1_𝑖+1,𝑗𝛿_𝑥𝑈^𝑛−1_𝑖+1,𝑗−Δ𝑡^𝑛−12ℎ^𝑛−1_𝑗,𝑦Γ^𝑛−1_{𝑖+1,𝑗+1/2}−Γ^𝑛−1_{𝑖+1,𝑗−1/2}, (2.19) where Γ^𝑛−1_{𝑖,𝑗+1/2}=𝐻_𝑔𝑈^{𝐿,𝑛−1}_{𝑖,𝑗+1/2},𝑈^{𝑅,𝑛−1}_{𝑖,𝑗+1/2},𝑈^{𝐿,𝑛−1}_{𝑖,𝑗+1/2}=𝑈^𝑛−1_𝑖,𝑗+ℎ^𝑛−1_𝑗,𝑦2𝛿_𝑦𝑈^𝑛−1_𝑖,𝑗,𝑈^{𝑅,𝑛−1}_{𝑖,𝑗+1/2}=𝑈^𝑛−1_𝑖,𝑗+1+ℎ^𝑛−1_𝑗+1,𝑦2𝛿_𝑦𝑈^𝑛−1_𝑖,𝑗+1. (2.20)

Third Step
With the above definitions, (𝑔^𝑛−1_{1,𝑖+1/2,𝑗},𝑔^𝑛−1_{2,𝑖,𝑗+1/2}) is calculated as follows.
If 𝑐₁(𝑈^𝑛−1_{𝑖+1/2,𝑗})>0, define 𝑔^𝑛−1_{1,𝑖+1/2,𝑗}≜𝐻_𝑔1𝑈^{𝐿,𝑛−1}_{𝑖+1/2,𝑗},𝑈^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}≜𝑐₁𝑈^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}𝑈^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}. (2.21)
If 𝑐₁(𝑈^𝑛−1_{𝑖+1/2,𝑗})≤0, define 𝑔^𝑛−1_{1,𝑖+1/2,𝑗}≜𝐻_𝑔1𝑈^{𝑅,𝑛−1}_{𝑖+1/2,𝑗},𝑈^{𝑅,𝑛−1}_{𝑖+1/2,𝑗}≜𝑐₁𝑈^{𝑅,𝑛−1}_{𝑖+1/2,𝑗}𝑈^{𝑅,𝑛−1}_{𝑖+1/2,𝑗}. (2.22) For 𝑔^𝑛−1_{2,𝑖,𝑗+1/2}, the definition is similar.
The equations (2.19)–(2.22) hold for elements at least one element away from the boundary. At the left and right boundaries, we can set 𝑈^{𝐿,𝑛−1}_1/2,𝑗=𝑈_𝐼^{𝑅,𝑛−1}_{𝑛−1+1/2,𝑗}=0,𝑗=1,2,…,𝐽^𝑛−1,𝑛=1,2,…,𝑁, (2.23) and in the slope calculation procedure, we define 𝛿_𝑥𝑈^𝑛−1_1,𝑗 and 𝛿_𝑥𝑈_𝐼^𝑛−1_{𝑛−1,𝑗} using (2.12) with ▵^𝑥₋𝑈^𝑛−1_1,𝑗=2𝑈^𝑛−1_1,𝑗ℎ^𝑛−1_1,𝑥,▵^𝑥₊𝑈_𝐼^𝑛−1_{𝑛−1,𝑗}=−2𝑈_𝐼^𝑛−1_{𝑛−1,𝑗}ℎ_𝐼^𝑛−1_{𝑛−1,𝑥}. (2.24) On the bottom and top boundaries, we set Γ^𝑛−1_𝑖,1/2=Γ^𝑛−1_{𝑖,𝐽𝑛−1+1/2}=0,𝑖=1,2,…,𝐼^𝑛−1, (2.25) in (2.19).

3. Projections and Lemma

We introduce three projections and a lemma to obtain error estimates. Let {𝑢^𝑛,𝐳^∼𝑛}∶𝐽→𝑊^𝑛_ℎ×𝑉^𝑛_ℎ be the projection of {𝑢,𝐳} in mixed finite element space such that 𝐳^∼𝑛,𝐯^𝑛−𝑢^𝑛,∇⋅𝐯^𝑛=0,∀𝐯^𝑛∈𝑉^𝑛_ℎ, (3.1) and define 𝜋𝐳^𝑛 as the 𝜋-projection of 𝐳^𝑛 such that ∇⋅𝐳^𝑛−𝜋𝐳^𝑛,𝑤^𝑛=0,∀𝑤^𝑛∈𝑊^𝑛_ℎ. (3.2)

These projections satisfy ([21, 22]): 𝑢^𝑛−𝑢^𝑛+𝑧^𝑛−𝐳^∼𝑛_𝐻(div)≤𝐾ℎ^𝑛,𝜕𝑢𝜕𝑡^𝑛−𝑢^𝑛+𝜕𝜕𝑡𝐳^𝑛−𝐳^∼𝑛_𝐻(div)≤𝐾ℎ^𝑛, (3.3) and it is easy to see that ‖𝑢^𝑛‖^∞_𝐿≤𝐾.

Let 𝐠^∼∈𝑉^𝑛_ℎ satisfies that the value of 𝐠^∼⋅𝜸 at the midpoint of the boundaries is equal to the average value of 𝐠⋅𝜸 on boundaries, that is, ∇⋅𝐠^∼𝑛−𝐠^𝑛,𝑤^𝑛=0,∀𝑤^𝑛∈𝑊^𝑛_ℎ. (3.4) Let 𝐠^∼=(𝑔₁,𝑔₂). From (3.4), we can define 𝑔₁𝑥_𝑖+1/2,𝑦_𝑗,𝑡^𝑛≜1ℎ^𝑛_𝑗,𝑦^𝑦𝑗+1/2_{𝑦𝑗−1/2}𝑐₁𝑥_𝑖+1/2,𝑦_𝑗,𝑡^𝑛,𝑢𝑢𝑥_𝑖+1/2,𝑦,𝑡^𝑛≜1𝑑𝑦ℎ^𝑛_𝑗,𝑦^𝑦𝑗+1/2_{𝑦𝑗−1/2}𝑔₁𝑢𝑥_𝑖+1/2,𝑦,𝑡^𝑛𝑑𝑦. (3.5) Then, 𝑔₂(𝑥_𝑖,𝑦_𝑗+1/2,𝑡^𝑛) can be defined similarly.

Lemma 3.1. Let 𝜉=𝑈−𝑢 and assume that 𝑢 is sufficiently smooth. Then 𝐆^𝑛−1⋅,𝑈^𝑛−1−𝐠^∼𝑛−1≤𝐾𝜉^𝑛−1+ℎ+Δ𝑡. (3.6)

Proof. Firstly, construct 𝑢^{𝐿,𝑛−1}_{𝑖+1/2,𝑗},𝑢^{𝑅,𝑛−1}_{𝑖+1/2,𝑗} from 𝑢 given in (3.1) and (3.2) by 𝑢^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}=𝑢^𝑛−1_𝑖,𝑗+ℎ^𝑛−1_𝑖,𝑥2−Δ𝑡^𝑛2𝑔_1𝑢𝑢^𝑛−1_𝑖,𝑗𝛿_𝑥𝑢^𝑛−1_𝑖,𝑗−Δ𝑡^𝑛2ℎ^𝑛−1_𝑗,𝑦Γ^𝑛−1_{𝑖,𝑗+1/2}(𝑢)−Γ^𝑛−1_{𝑖,𝑗−1/2}(𝑢), (3.7) at interior edges, where 𝛿_𝑥𝑢^𝑛−1_𝑖,𝑗=▵^𝑥₊𝑢^𝑛−1_𝑖,𝑗,if𝛿_𝑥𝑈^𝑛−1_𝑖,𝑗=▵^𝑥₊𝑈^𝑛−1_𝑖,𝑗,▵^𝑥₋𝑢^𝑛−1_𝑖,𝑗Γ,otherwise,_{𝑛−1𝑖,𝑗+1/2}(𝑢)=𝐻_𝑔𝑢^{𝐿,𝑛−1}_{𝑖,𝑗+1/2},𝑢^{𝑅,𝑛−1}_{𝑖,𝑗+1/2},𝑢^{𝐿,𝑛−1}_{𝑖,𝑗+1/2}=𝑢^𝑛−1_𝑖,𝑗+ℎ^𝑛−1_𝑗,𝑦2𝛿_𝑦𝑢^𝑛−1_𝑖,𝑗,𝑢^{𝑅,𝑛−1}_{𝑖,𝑗+1/2}=𝑢^𝑛−1_𝑖,𝑗+1+ℎ^𝑛−1_𝑗+1,𝑦2𝛿_𝑦𝑢^𝑛−1_𝑖,𝑗+1. (3.8) And define 𝑢^{𝑅,𝑛−1}_{𝑖+1/2,𝑗} similarly.
Then, assuming 𝑔₁ is twice differentiable and Lipschitz continuous, and using the Lipschitz continuity of the Godunov flux and (2.14), it can be shown that 𝑈^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}−𝑢^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}≤𝐾^𝑖+2_{𝑙=𝑖−1}^𝑗+2_{𝑚=𝑗−1}𝑈^𝑛−1_𝑙,𝑚−𝑢^𝑛−1_𝑙,𝑚. (3.9) A similar bound holds for |𝑈^{𝑅,𝑛−1}_{𝑖+1/2,𝑗}−𝑢^{𝑅,𝑛−1}_{𝑖+1/2,𝑗}|.
At the boundaries, we follow (2.23) and (2.25) and define 𝑢^{𝐿,𝑛−1}_1/2,𝑗=𝑢_𝐼^{𝑅,𝑛−1}_{𝑛−1+1/2,𝑗}=0,𝑗=1,2,…,𝐽^𝑛−1,Γ^𝑛−1_𝑖,1/2(𝑢)=Γ^𝑛−1_{𝑖,𝐽𝑛−1+1/2}(𝑢)=0,𝑖=1,2,…,𝐼^𝑛−1. (3.10) Define 𝑢^{𝐿,𝑛−1}_{𝑖+1/2,𝑗},𝑢^{𝑅,𝑛−1}_{𝑖+1/2,𝑗} analogous to 𝑢^{𝐿,𝑛−1}_{𝑖+1/2,𝑗},𝑢^{𝑅,𝑛−1}_{𝑖+1/2,𝑗}.
If 𝑐₁(𝑈^𝑛−1_{𝑖+1/2,𝑗})>0, we consider 𝑔^𝑛−1_{1,𝑖+1/2,𝑗}−𝑔₁𝑥_𝑖+1/2,𝑦_𝑗,𝑡^𝑛−1=𝐻_𝑔1𝑈^{𝐿,𝑛−1}_{𝑖+1/2,𝑗},𝑈^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}−1ℎ^𝑛−1_𝑗,𝑦^{𝑦𝑗+1/2𝑦}_𝑗−1/2𝑔₁𝑢𝑥_𝑖+1/2,𝑦,𝑡^𝑛−1=𝐻𝑑𝑦_𝑔1𝑈^{𝐿,𝑛−1}_{𝑖+1/2,𝑗},𝑈^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}−𝐻_𝑔1𝑢^{𝐿,𝑛−1}_{𝑖+1/2,𝑗},𝑢^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}+𝐻_𝑔1𝑢^{𝐿,𝑛−1}_{𝑖+1/2,𝑗},𝑢^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}−𝐻_𝑔1𝑢^{𝐿,𝑛−1}_{𝑖+1/2,𝑗},𝑢^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}+𝐻_𝑔1𝑢^{𝐿,𝑛−1}_{𝑖+1/2,𝑗},𝑢^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}−𝑔₁𝑢^𝑛−1/2_{𝑖+1/2,𝑗}+𝑔₁𝑢^𝑛−1/2_{𝑖+1/2,𝑗}−1ℎ^𝑛−1_𝑗,𝑦^{𝑦𝑗+1/2𝑦}_𝑗−1/2𝑔₁𝑢𝑥_𝑖+1/2,𝑦,𝑡^𝑛−1𝑑𝑦=𝑅₁+𝑅₂+𝑅₃+𝑅₄. (3.11)
By the Lipschitz continuity of 𝐻_𝑔1 and (3.9), we have 𝑅₁≤𝐾𝑈^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}−𝑢^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}+𝑈^{𝑅,𝑛−1}_{𝑖+1/2,𝑗}−𝑢^{𝑅,𝑛−1}_{𝑖+1/2,𝑗}≤𝐾^𝑖+2_{𝑙=𝑖−1}^𝑗+2_{𝑚=𝑗−1}𝑈^𝑛−1_𝑙,𝑚−𝑢^𝑛−1_𝑙,𝑚. (3.12) The similar arguments are applied to 𝑅₂ to get 𝑅₂≤𝐾^𝑖+2_{𝑙=𝑖−1}^𝑗+2_{𝑚=𝑗−1}𝑢^𝑛−1_𝑙,𝑚−𝑢^𝑛−1_{𝑖+1/2,𝑗}. (3.13) By the consistency of 𝐻_𝑔1, we see 𝑅₃≤𝐾𝑢^{𝐿,𝑛−1}_{𝑖+1/2,𝑗}−𝑢^𝑛−1/2_{𝑖+1/2,𝑗}=𝑂Δ𝑡+ℎ². (3.14) For 𝑅₄, it is easy to have 𝑅₄≤𝐾Δ𝑡+ℎ². (3.15)
Using (3.3), (3.12)–(3.15), and equivalence of norms, we have 𝑔₁^𝑛−1−𝑔₁^𝑛−1≤𝐾_𝑗𝑖𝑔^𝑛−1_{1,𝑖+1/2,𝑗}−𝑔^𝑛−1_{1,𝑖+1/2,𝑗}²ℎ^𝑛−1_{𝑖+1/2,𝑥}ℎ^𝑛−1_𝑗,𝑦≤𝐾𝜉^𝑛−1+Δ𝑡+ℎ. (3.16)
The similar approach can be used when 𝑐₁(𝑈^𝑛−1_{𝑖+1/2,𝑗})≤0. Defining 𝑔₂ and 𝑔₂ similarly and following analogous arguments yields 𝑔₂^𝑛−1−𝑔₂^𝑛−1≤𝐶𝜉^𝑛−1+Δ𝑡+ℎ. (3.17) Thus, Lemma 3.1 holds.

4. Error Estimates

At time level 𝑡^𝑛, for all 𝑤^𝑛∈𝑊^𝑛_ℎ, 𝐯^𝐧∈𝑉^𝑛_ℎ, the exact solutions satisfy 𝑢^𝑛−𝑢^𝑛−1Δ𝑡^𝑛,𝑤^𝑛+∇⋅𝐳^𝑛,𝑤^𝑛+𝑢^𝑛𝐜𝑢^𝑛⋅𝐳^𝑛,𝑤^𝑛=𝑓𝑢^𝑛,𝑤^𝑛−∇⋅𝐠^𝑛,𝑤^𝑛−𝜌^𝑛,𝑤^𝑛,𝐳,𝐯^𝐧=𝑢^𝑛,∇⋅𝐯^𝐧,𝐳^𝑛,𝐯^𝐧=𝑏𝑢^𝑛𝐳^𝑛+𝑎𝑢^𝑛𝐳^𝑛−𝐳^𝑛−1Δ𝑡^𝑛,𝐯^𝐧+𝐫^𝑛,𝐯^𝐧, (4.1) where 𝜌^𝑛=𝑢^𝑛_𝑡−(𝑢^𝑛−𝑢^𝑛−1)/Δ𝑡^𝑛, 𝐫^𝑛=𝑎(𝑢^𝑛)(𝐳^𝑛_𝑡−(𝐳^𝑛−𝐳^𝑛−1)/Δ𝑡^𝑛).

Let 𝜉_𝑢=𝑈−𝑢,𝛽_𝑢=𝑢−𝑢,𝝃_𝐳=𝐙−𝐳^∼,𝜷_𝐳=𝐳−𝐳^∼,𝝃_𝐳=𝐙−𝜋𝐳,𝜷_𝐳=𝐳−𝜋𝐳. (4.2)

Using projections (3.1) and (3.2), we subtract (2.9) from (4.1) to get 𝜉^𝑛_𝑢−𝜉_𝑢^𝑛−1Δ𝑡^𝑛,𝑤^𝑛+∇⋅𝝃^𝑛_𝐳,𝑤^𝑛+𝑈^𝑛𝐜𝑈^𝑛⋅𝐙^𝑛−𝑢^𝑛𝐜𝑢^𝑛⋅𝐳^𝑛,𝑤^𝑛=𝑓𝑈^𝑛−𝑓𝑢^𝑛,𝑤^𝑛+𝜌^𝑛,𝑤^𝑛+∇⋅𝐠^∼𝑛−𝐆^𝑛−1,𝑤^𝑛+𝛽^𝑛_𝑢−𝛽_𝑢^𝑛−1Δ𝑡^𝑛,𝑤^𝑛,𝝃^𝑛_𝐳,𝐯^𝐧=𝜉^𝑛_𝑢,∇⋅𝐯^𝐧,𝝃^𝑛_𝐳,𝐯^𝐧=𝑏𝑈^𝑛𝐙^𝐧−𝑏𝑢^𝑛𝐳^𝐧,𝐯^𝐧−𝐫^𝑛,𝐯^𝐧−𝜷^𝑛_𝐳,𝐯^𝐧+𝑎𝑈^𝑛𝐙^𝐧−𝐙^𝐧−𝟏Δ𝑡^𝑛−𝑎𝑢^𝑛𝐳^𝐧−𝐳^𝐧−𝟏Δ𝑡^𝑛,𝐯^𝐧, (4.3) where the last term of the first equation in (4.3) is related to the changing meshes. If the meshes do not change, this term will be zero.

Taking 𝑤=𝜉^𝑛_𝑢, 𝐯=𝝃^𝑛_𝐳 in the second equation and 𝐯=𝝃^𝑛_𝐳 in the third, then adding them together, we have 𝜉^𝑛_𝑢−𝜉_𝑢^𝑛−1Δ𝑡^𝑛,𝜉^𝑛_𝑢+𝑏𝑈^𝑛𝐙^𝑛−𝑏𝑢^𝑛𝐳^𝑛,𝝃^𝑛_𝐳+𝑎𝑈^𝑛𝐙^𝑛−𝐙^𝑛−1Δ𝑡^𝑛,𝝃^𝑛_𝐳−𝑎𝑢^𝑛𝐳^𝑛−𝐳^𝑛−1Δ𝑡^𝑛,𝝃^𝑛_𝐳=𝑓𝑈^𝑛−𝑓𝑢^𝑛,𝜉^𝑛_𝑢+𝜌^𝑛,𝜉^𝑛_𝑢+∇⋅𝐠^∼𝑛−𝐆^𝑛−1,𝜉^𝑛_𝑢−𝑈^𝑛𝐜𝑈^𝑛⋅𝐙^𝑛−𝑢^𝑛𝐜𝑢^𝑛⋅𝐳^𝑛,𝜉^𝑛_𝑢+𝐫^𝑛,𝝃^𝑛_𝐳−𝜷^𝑛_𝐳,𝝃^𝑛_𝐳+𝛽^𝑛_𝑢−𝛽_𝑢^𝑛−1Δ𝑡^𝑛,𝜉^𝑛_𝑢. (4.4) For (4.4), we see

(I) 𝑏𝑈^𝑛𝐙^𝑛−𝑏𝑢^𝑛𝐳^𝑛,𝝃^𝑛_𝐳=𝑏𝑈^𝑛𝝃^𝑛_𝐳,𝝃^𝑛_𝐳−𝑏𝑈^𝑛𝜷^𝑛_𝐳,𝝃^𝑛_𝐳+𝑏𝑈^𝑛−𝑏𝑢^𝑛𝐳^𝑛,𝝃^𝑛_𝐳, (4.5)

(II) 𝑎𝑈^𝑛𝐙^𝑛−𝐙^𝑛−1Δ𝑡^𝑛,𝝃^𝑛_𝐳−𝑎𝑢^𝑛𝐳^𝑛−𝐳^𝑛−1Δ𝑡^𝑛,𝝃^𝑛_𝐳=𝑎𝑈^𝑛𝝃^𝑛_𝐳−𝝃_𝐳^𝑛−1Δ𝑡^𝑛,𝝃^𝑛_𝐳−𝑎𝑈^𝑛𝜷^𝑛_𝐳−𝜷_𝐳^𝑛−1Δ𝑡^𝑛,𝝃^𝑛_𝐳+𝑎𝑈^𝑛−𝑎𝑢^𝑛𝐳^𝑛−𝐳^𝑛−1Δ𝑡^𝑛,𝝃^𝑛_𝐳, (4.6)

(III) 𝑈^𝑛𝐜𝑈^𝑛⋅𝐙^𝑛−𝑢^𝑛𝐜𝑢^𝑛⋅𝐳^𝑛,𝜉^𝑛_𝑢=𝜉^𝑛_𝑢−𝛽^𝑛_𝑢𝐜𝑈^𝑛⋅𝐙^𝑛,𝜉^𝑛_𝑢+𝑢^𝑛𝐜𝑈^𝑛⋅𝝃^𝑛_𝐳−𝜷^𝑛_𝐳,𝜉^𝑛_𝑢+𝑢^𝑛𝐜𝑈^𝑛−𝐜𝑢^𝑛⋅𝐳^𝑛,𝜉^𝑛_𝑢. (4.7) The last second term in (4.6) is related to the changing meshes. If the meshes do not change, this term will be zero.

Substitute (4.5)–(4.7) into (4.4) to get 𝜉^𝑛_𝑢−𝜉_𝑢^𝑛−1Δ𝑡^𝑛,𝜉^𝑛_𝑢+𝑏𝑈^𝑛𝝃^𝑛_𝐳,𝝃^𝑛_𝐳+𝑎𝑈^𝑛𝝃^𝑛_𝐳−𝝃_𝐳^𝑛−1Δ𝑡^𝑛,𝝃^𝑛_𝐳=𝑓𝑈^𝑛−𝑓𝑢^𝑛,𝜉^𝑛_𝑢+𝜌^𝑛,𝜉^𝑛_𝑢+∇⋅𝐠^∼𝑛−𝐆^𝑛−1,𝜉^𝑛_𝑢+𝐫^𝑛,𝝃^𝑛_𝐳−𝜷^𝑛_𝐳,𝝃^𝑛_𝐳−𝜉^𝑛_𝑢−𝛽^𝑛_𝑢𝐜𝑈^𝑛⋅𝐙^𝑛,𝜉^𝑛_𝑢−𝑢^𝑛𝐜𝑈^𝑛⋅𝝃^𝑛_𝐳−𝜷^𝑛_𝐳,𝜉^𝑛_𝑢−𝑢^𝑛𝐜𝑈^𝑛−𝐜𝑢^𝑛⋅𝐳^𝑛,𝜉^𝑛_𝑢+𝛽^𝑛_𝑢−𝛽_𝑢^𝑛−1Δ𝑡^𝑛,𝜉^𝑛_𝑢+𝑏𝑈^𝑛𝜷^𝑛_𝐳,𝝃^𝑛_𝐳+𝑏𝑢^𝑛−𝑏𝑈^𝑛𝐳^𝑛,𝝃^𝑛_𝐳+𝑎𝑢^𝑛−𝑎𝑈^𝑛𝐳^𝑛−𝐳^𝑛−1Δ𝑡^𝑛,𝝃^𝑛_𝐳+𝑎𝑈^𝑛𝜷^𝑛_𝑧−𝜷_𝑧^𝑛−1Δ𝑡^𝑛,𝝃^𝑛_𝐳=𝑇₁+𝑇₂+⋯+𝑇₁₃. (4.8) Furthermore, for the first and third terms on the left-hand side of (4.8), we have 𝜉^𝑛_𝑢−𝜉_𝑢^𝑛−1Δ𝑡^𝑛,𝜉^𝑛_𝑢=12Δ𝑡^𝑛𝜉^𝑛_𝑢,𝜉^𝑛_𝑢−𝜉_𝑢^𝑛−1,𝜉_𝑢^𝑛−1+12Δ𝑡^𝑛𝜉^𝑛_𝑢−𝜉_𝑢^𝑛−12,𝑎𝑈^𝑛𝝃^𝑛_𝐳−𝝃_𝐳^𝑛−1Δ𝑡^𝑛,𝝃^𝑛_𝐳≥12Δ𝑡^𝑛𝑎𝑈^𝑛𝝃^𝑛_𝐳,𝝃^𝑛_𝐳−𝑎𝑈^𝑛−1𝝃_𝐳^𝑛−1,𝝃_𝐳^𝑛−1−𝑎₂2𝝃_𝐳^𝑛−12+𝑎₀2Δ𝑡^𝑛𝝃^𝑛_𝐳−𝝃_𝐳^𝑛−12. (4.9) By (4.9), we modify (4.8) as 12Δ𝑡^𝑛𝜉^𝑛_𝑢,𝜉^𝑛_𝑢−𝜉_𝑢^𝑛−1,𝜉_𝑢^𝑛−1+12Δ𝑡^𝑛𝜉^𝑛_𝑢−𝜉_𝑢^𝑛−12+𝑏𝑈^𝑛𝝃^𝑛_𝐳,𝝃^𝑛_𝐳+12Δ𝑡^𝑛𝑎𝑈^𝑛𝝃^𝑛_𝐳,𝝃^𝑛_𝐳−𝑎𝑈^𝑛−1𝝃_𝐳^𝑛−1,𝝃_𝐳^𝑛−1+𝑎₀2Δ𝑡^𝑛𝝃^𝑛_𝐳−𝝃_𝐳^𝑛−12≤𝑎₂2𝝃_𝐳^𝑛−12+𝑇₁+𝑇₂+⋯+𝑇₁₃. (4.10)

By the Lipschitz continuity of 𝑓(𝑢), 𝑎(𝑢), 𝑏(𝑢), 𝜕𝐜/𝜕𝑢, and (3.3), the following terms on the right-hand side of (4.10) can be obtained: 𝑇₁≤𝐾𝑈^𝑛−𝑢^𝑛2+𝐾𝜉^𝑛_𝑢²≤𝐾𝜉^𝑛_𝑢²+𝐾ℎ²,𝑇₂≤𝐾Δ𝑡^𝑛𝜕²𝑢𝜕𝑡²²_{𝐿2(𝑡𝑛−1,𝑡𝑛;𝐻1)}+𝐾𝜉^𝑛_𝑢²,𝑇₄≤𝐾Δ𝑡^𝑛𝜕²𝑢𝜕𝑡²²_{𝐿2(𝑡𝑛−1,𝑡𝑛;𝐻1)}+𝐾𝝃^𝑛_𝐳²,𝑇₅≤𝐾𝝃^𝑛_𝐳²+𝐾ℎ²,𝑇₆≤𝐾𝜉^𝑛_𝑢²+𝐾ℎ²,𝑇₇≤𝐾𝜉^𝑛_𝑢²+𝐾𝝃^𝑛_𝐳²+𝐾ℎ²,𝑇₈≤𝐾𝑈^𝑛−𝑢^𝑛2+𝐾𝜉^𝑛_𝑢²≤𝐾𝜉^𝑛_𝑢²+𝐾ℎ²,𝑇₁₀≤𝐾𝑏₁ℎ²+𝐾𝝃^𝑛_𝐳²,𝑇₁₁≤𝐾𝝃^𝑛_𝐳²+𝐾𝜉^𝑛_𝑢²+𝐾ℎ²,𝑇₁₂≤𝐾𝝃^𝑛_𝐳²+𝐾𝑈^𝑛−𝑢^𝑛2+𝐾𝐳^𝑛−𝐳^𝑛−1Δ𝑡^𝑛2≤𝐾𝝃^𝑛_𝐳²+𝐾𝜉^𝑛_𝑢²+𝐾Δ𝑡^𝑛−1𝜕𝑢𝜕𝑡²_𝐿²_(𝑡^𝑛−1_,𝑡^𝑛_;𝐻¹₎+𝐾ℎ². (4.11)

We turn to consider 𝑇₃. Letting 𝐯=𝐠^∼𝑛−𝐆^𝑛−1 in the second equation of (4.10), then we have 𝝃^𝑛_𝐳,𝐠^∼𝑛−𝐆^𝑛−1=𝜉^𝑛_𝑢,∇⋅𝐠^∼𝑛−𝐆^𝑛−1. (4.12) By Lemma 3.1, we have 𝑇₃≤12𝑏𝑈^𝑛𝝃^𝑛_𝐳,𝝃^𝑛_𝐳+𝐾𝑏𝐠₀^−1∼𝑛−𝐆^𝑛−12≤12𝑏𝑈^𝑛𝝃^𝑛_𝐳,𝝃^𝑛_𝐳+𝐾𝑏𝐠₀^−1∼𝑛−𝐠^∼𝑛−12+𝐠^∼𝑛−1−𝐆^𝑛−12≤12𝑏𝑈^𝑛𝝃^𝑛_𝐳,𝝃^𝑛_𝐳+𝐾𝑏₀⁻¹𝜉_𝑢^𝑛−12+ℎ²+Δ𝑡². (4.13)

Substituting (4.11) and (4.13) into (4.10), we have 12Δ𝑡^𝑛𝜉^𝑛_𝑢,𝜉^𝑛_𝑢−𝜉_𝑢^𝑛−1,𝜉_𝑢^𝑛−1+12Δ𝑡^𝑛𝜉^𝑛_𝑢−𝜉_𝑢^𝑛−12+12𝑏𝑈^𝑛𝝃^𝑛_𝐳,𝝃^𝑛_𝐳+𝑎₀2Δ𝑡^𝑛𝝃^𝑛_𝐳−𝝃_𝐳^𝑛−12+12Δ𝑡^𝑛𝑎^𝑛𝝃^𝑛_𝐳,𝝃^𝑛_𝐳−𝑎^𝑛−1𝝃_𝐳^𝑛−1,𝝃_𝐳^𝑛−1≤𝐾𝝃^𝑛_𝐳²+𝝃_𝐳^𝑛−12+𝐾𝜉^𝑛_𝑢²+𝜉_𝑢^𝑛−12+Δ𝑡²+ℎ²+𝐾Δ𝑡^𝑛𝜕𝑢𝜕𝑡²_𝐿²_(𝑡^𝑛−1_,𝑡^𝑛_;𝐻¹₎+𝜕²𝑢𝜕𝑡²²_{𝐿2(𝑡𝑛−1,𝑡𝑛;𝐻1)}+𝑇₉+𝑇₁₃. (4.14)

Let 𝑁 be the time step at which ‖𝜉^𝑛_𝑢‖ is maximum, that is, 𝜉^𝑁_𝑢²=max_{1≤𝑛≤𝑁∗}𝜉^𝑛_𝑢². (4.15) Multiplying (4.14) by 2Δ𝑡^𝑛 and summing on 𝑛 from 1 to 𝑁, we obtain 𝜉^𝑁_𝑢²+^𝑁_𝑛=1𝜉^𝑛_𝑢−𝜉_𝑢^𝑛−12+^𝑁_𝑛=1𝑏𝑈^𝑛𝝃^𝑛_𝐳,𝝃^𝑛_𝐳Δ𝑡^𝑛+𝑎₀^𝑁_𝑛=1𝝃^𝑛_𝐳−𝝃_𝐳^𝑛−12+𝑎₀𝝃^𝑁_𝐳²≤𝐾Δ𝑡²+ℎ²+𝐾^𝑁_𝑛=1𝜉^𝑛_𝑢²Δ𝑡^𝑛+𝜉⁰_𝑢²+𝑎₁𝝃⁰_𝐳²+𝐾^𝑁_𝑛=1𝝃^𝑛_𝐳²Δ𝑡^𝑛+^𝑁_𝑛=1𝛽^𝑛_𝑢−𝛽_𝑢^𝑛−1Δ𝑡^𝑛,𝜉^𝑛_𝑢Δ𝑡^𝑛+^𝑁_𝑛=1𝑎𝑈^𝑛𝜷^𝑛_𝐳−𝜷_𝐳^𝑛−1Δ𝑡^𝑛,𝝃^𝑛_𝐳Δ𝑡^𝑛. (4.16)

Assuming that the mesh is modified at most 𝑀 times, and 𝑀≤𝑀^∗, where 𝑀^∗ is independent of ℎ and Δ𝑡, we get [3]: ^𝑁_𝑛=1𝛽^𝑛_𝑢−𝛽_𝑢^𝑛−1Δ𝑡^𝑛,𝜉^𝑛_𝑢Δ𝑡^𝑛=^𝑁_𝑛=1𝛽^𝑛_𝑢,𝜉^𝑛_𝑢−𝛽_𝑢^𝑛−1,𝜉^𝑛_𝑢≤𝐾𝑀^∗ℎ²+14𝜉^𝑁_𝑢²,_𝑁𝑛=1𝑎𝑈^𝑛𝜷^𝑛_𝐳−𝜷_𝐳^𝑛−1Δ𝑡^𝑛,𝝃^𝑛_𝐳Δ𝑡^𝑛≤𝐾ℎ²+𝜀^𝑁_𝑛=1𝝃^𝑛_𝐳²Δ𝑡^𝑛, (4.17) where 𝜀 is a sufficiently small positive constant.

Substituting (4.17) into (4.16), we find 𝜉^𝑁_𝑢²+^𝑁_𝑛=1𝜉^𝑛_𝑢−𝜉_𝑢^𝑛−12+^𝑁_𝑛=1𝑏𝑈^𝑛𝝃^𝑛_𝐳,𝝃^𝑛_𝐳Δ𝑡^𝑛+𝑎₀^𝑁_𝑛=1𝝃^𝑛_𝐳−𝝃_𝐳^𝑛−12+𝑎₀𝝃^𝑁_𝐳²≤𝐾Δ𝑡²+ℎ²+𝑀^∗ℎ²+𝜉⁰_𝑢²+𝑎₁𝝃⁰_𝐳²+𝐾^𝑁_𝑛=1𝜉^𝑛_𝑢²Δ𝑡^𝑛+𝐾^𝑁_𝑛=1𝝃^𝑛_𝐳²Δ𝑡^𝑛+𝜀^𝑁_𝑛=1𝝃^𝑛_𝐳²Δ𝑡^𝑛+14𝜉^𝑁_𝑢². (4.18) Using the discrete Gronwall's lemma, we have max_𝑛𝜉^𝑛_𝑢²≤𝐾ℎ²+Δ𝑡²+𝑀^∗ℎ². (4.19) Similarly as (4.14), we can derive max_𝑛𝝃^𝑛_𝐳²≤𝐾ℎ²+Δ𝑡²+𝑀^∗ℎ². (4.20)

By the triangle inequality and (3.3), we can obtain the following result.

Theorem 4.1. Assuming that the coefficients satisfy condition (1.2), the mesh is modified at most 𝑀 times, 𝑀≤𝑀^∗, and 𝑢∈𝐿²(𝐻¹), 𝑢_𝑡∈𝐿²(𝐻¹), 𝑢_𝑡𝑡∈𝐿²(𝐻¹), then we have max_𝑛𝑢^𝑛−𝑈^𝑛≤𝐾(ℎ+Δ𝑡)+𝐾𝑀^∗ℎ,max_𝑛𝐳^𝑛−𝐙^𝑛≤𝐾(ℎ+Δ𝑡)+𝐾𝑀^∗ℎ, (4.21) where 𝐾 and 𝑀^∗ are positive constants independent of ℎ and Δ𝑡.

5. Numerical Example

In this section, we give numerical experiment of the following nonlinear equations to illustrate the efficiency of the Godunov-mixed finite element method on changing meshes: 𝑢_𝑡=𝑎(𝑥,𝑡)𝑢_𝑥𝑡+𝑏(𝑥,𝑡,𝑢)𝑢_𝑥𝑥−𝑐(𝑥,𝑡)𝑢_𝑥(),()∈()×(],𝑢()()[],𝑢()[],+𝑓𝑥,𝑡,𝑢𝑥,𝑡0,10,10,𝑡=𝑢1,𝑡=0,𝑡∈0,1⋅,0=sin𝜋𝑥,𝑥∈0,1 (5.1) where (𝑥,𝑡)∈[0,1]×[0,1], 𝑎(𝑥,𝑡)=𝑡²(𝑥+0.25)+0.25, 𝑏(𝑥,𝑡,𝑢)=0.005𝑢𝑡²(𝑥+0.05), 𝑐(𝑥,𝑡)=0.05𝑡(𝑥+1)+1, and 𝑓(𝑥,𝑡,𝑢) is chosen properly so that the exact solution is 𝑢=𝑒^−𝑡sin𝜋𝑥.

In the numerical example, ℎ and Δ𝑡 change as the following four cases.

Case 1. 𝑡∈[0,0.4], set ℎ=Δ𝑡=0.1; 𝑡∈(0.4,0.6], set ℎ=Δ𝑡=0.05; 𝑡∈(0.6,1.0], set ℎ=Δ𝑡=0.1 and calculate 40 steps at every time interval.

Case 2. 𝑡∈[0,0.4], set ℎ=Δ𝑡=0.05; 𝑡∈(0.4,0.6], set ℎ=Δ𝑡=0.025; 𝑡∈(0.6,1.0], set ℎ=Δ𝑡=0.05 and calculate 80 steps at every time interval.

Case 3. 𝑡∈[0,0.4], set ℎ=Δ𝑡=0.025; 𝑡∈(0.4,0.6], set ℎ=Δ𝑡=0.0125; 𝑡∈(0.6,1.0], set ℎ=Δ𝑡=0.025 and calculate 160 steps at every time interval.

Case 4. 𝑡∈[0,0.4], set ℎ=Δ𝑡=0.0125; 𝑡∈(0.4,0.6], set ℎ=Δ𝑡=0.00625; 𝑡∈(0.6,1.0], set ℎ=Δ𝑡=0.0125 and calculate 320 steps at every time interval.
The numerical solutions 𝑈^𝑛, 𝐙^𝑛 are calculated and the 𝐿² error estimates of 𝑈^𝑛−𝑢^𝑛,𝐙^𝑛−𝐳^𝑛 are obtained; see Tables 1 and 2.

Table 1 

𝐿²-norm estimates of 𝑈^𝑛−𝑢^𝑛.

Table 2 

𝐿²-norm estimates of 𝐙^𝑛−𝐳^𝑛.

In Table 3, convergence rates of 𝑈^𝑛−𝑢^𝑛,𝐙^𝑛−𝐳^𝑛 are given. The figures of convergence rate are shown by Figure 1. The convergence rate of 𝑈^𝑛−𝑢^𝑛 in Figure 1(a) is one order, and in Figure 1(b) the convergence rate of 𝐙^𝑛−𝐳^𝑛 is a little smaller than one order at the beginning. But the convergence rate 𝐙^𝑛−𝐳^𝑛 will be close to one order when ℎ minishes, which is consistent with the analysis in this paper.

Table 3 

Convergence rate of 𝐿²-norm.

(a)  Convergence rates of 𝑢

(b) Convergence rates of 𝐳∼

Figure 1 

Convergence rate figures.

For Case 4, we compare the exact solution {𝑢,𝐳^∼} with the approximate solution {𝑈,𝐙^∼} at time 𝑡=0.25, 0.5, and 0.75, respectively; see Figure 2. From Figures 2(a) and 2(b), we can see that the approximate solutions are very close to the exact solutions.

(a) Comparing 𝑢 with 𝑈

(b) Comparing 𝐳∼ with 𝐙∼

Figure 2 

Compare figures for Case 4 at time 𝑡=0.25, 0.5, and 0.75.

Acknowledgments

The author would like to thank the referees for their constructive comments leading to an improved presentation of this paper. This work was supported by the NSF of China (no. 11271231), the NSF of Shandong Province (no. ZR2010AQ019), the Scientific Research Award Fund for Excellent Middle-Aged and Young Scientists of Shandong Province (no. BS2009NJ003), and Science and Technology Development Planning Project of Shandong Province (no. 2011GGH20118), China.