A weak Galerkin finite element method for nonlinear conservation laws

1.   Introduction

The nonlinear hyperbolic equation of conservation laws is considered: seeking an unknown function $u$ satisfying

with $I = (0, 1)$ and a periodic boundary condition. For simplicity, the nonlinear flux function $f(u)$ is assumed to be smooth enough.

The Runge-Kutta discontinuous Galerkin (RKDG) method has been developed for solving time-dependent nonlinear conservation laws [1,2]. The RKDG method, by name, uses DG method for spacial discretization and explicit high order Runge-Kutta method for time discretizations. The stability and error analysis of the RKDG method has been studied in [22,23] for the second and the third order explicit total variation diminishing Runge-Kutta method. A discontinuous Galerkin method with Lagrange multiplier (DGLM) has been developed in [6,7] for nonlinear conservation laws with backward Euler method for time discretization. Lagrange multipliers are introduced on each element so that they are the only globally coupled variables in the resulting system. The final global system of the DGLM has fewer numbers of coupled unknowns than the usual DG methods.

The most finite element methods for conservation law employ purely upwind or general monotone fluxes. In [11], discontinuous Galerkin methods using more general upwind-biased numerical fluxes have been investigated for time-dependent linear conservation laws. Optimal order of convergence rate has been obtained. As pointed out in [11], purely upwind fluxes may be difficult to construct for complex systems.

Weak Galerkin methods refer to general finite element techniques for partial differential equations and were first introduced in [18,19] for second order elliptic equations. Weak Galerkin methods make use of discontinuous piecewise polynomials on finite element partitions with arbitrary shape of polygons and polyhedrons. The weak Galerkin methods have been applied to solve various PDEs such as second order elliptic equations, biharmonic equations, Stokes equations, parabolic equations, second order hyperbolic equations, Maxwell's equations and singularly perturbed convection-diffusion-reaction problems [8,9,10,12,13,14,15,16,17,19,20]. A least-squares based weak Galerkin method is presented for stationary linear hyperbolic equations [21].

The objective of this work is to develop a weak Galerkin finite element method for the time-dependent nonlinear conservation laws (1.1)-(1.2), with the explicit first order Euler method and the third order explicit TVD Runge-Kutta method for time discretization. Similar to [11], this new WG formulation provides a class of finite element methods featuring two built-in parameters $\lambda_1$ and $\lambda_2$. By tuning these parameters, different schemes can be obtained for solving the problem (1.1)-(1.2) including purely upwinding scheme. However, unlike the method in [11], our new WG method can be used for the time-dependent nonlinear conservation laws. The stability is derived for the semi-discretized WG method. For the forward Euler time discrete WG method, the $L^2$ error estimate of $O(h^{k+\frac12}+\tau)$ is derived in general and convergence rate of $O(h^{k+1}+\tau)$ is obtained for some special combination of the parameters. The temporal-spatial CFL condition $\tau< Ch^2$ is necessary in the error analysis for the first order forward Euler method. If the third order explicit TVDRK time discrete scheme [3] is used, the $L^2$ error estimate of $O(h^{k+\frac12}+\tau^3)$ is proved under the CFL condition $\tau < Ch$.

The rest of the paper is organized as follows. In Section 2, a WG finite element is proposed for spatial discretization. The stability is derived for semi-discretized WG method. The forward Euler time discretized WG scheme and its error analysis are presented in Section 3. Error analysis of the WG method with third order explicit TVDRK time discretization can be found in Section 4. Numerical experiments are presented in Section 5 to support the theoretical results. We end the paper with a conclusion.

The usual notation of norms in Sobolev spaces will be used. For any integer $s\geq 0$, let $H^s(D)$ represent the well-known Sobolev space equipped with the norm $\|\cdot\|_{s,D}$, which consists of functions with (distributional) derivatives of order no more than $s$ in $L^2(\Omega)$. Next, denote by $(\cdot,\cdot)_D$ the scalar inner product on $L^2(D)$ and $\|\cdot\|_D$ denotes the associated $L^2$ norm. Furthermore, let $\|\cdot\|_{\infty, D}$ be the norm on $L^{\infty}(D)$. If $D = I$, we omit this subscript.

2.   Semi-discrete weak Galerkin scheme and its stability

In this section, we introduce a weak Galerkin finite element method for solving the model problem (1.1)-(1.2).

2.1.   Semi-discrete WG scheme

Let ${{\mathcal T}}_h = \cup_{i = 1}^N I_i$ with $I_i = [x_{i-\frac12},x_{i+\frac12}]$ for $1\le i\le N$ where

Define

The weak Galerkin methods introduce a new way to define a function $v$, called weak function, which allows $v$ taking different forms in the interior and on the boundary of the element:

where ${I_i}^0$ is the interior of ${I_i}$, where ${I_i}$ is an element in ${{\mathcal T}}_h$. Since a weak function $v$ is formed by two parts $v_0$ and $v_b$, we write $v$ as $v = \{v_0,v_b\}$ in short without confusion.

Denote by $P_{k}({I_i})$ the set of polynomials on ${I_i}$ with degree no more than $k$. Let $V_h$ be a weak Galerkin finite element space consisting of weak function $v = \{v_0,v_b\}$ defined as follows for $k\ge 1$,

Denote by $v(x_{i+\frac12}^-)$ and $v(x_{i+\frac12}^+)$ the values of $v$ at $x_{i+\frac12}$ from the left element $I_i$ and the right element $I_{i+1}$, respectively. Further, the jump of $v$ at $x_{i+\frac12}$ is denoted as ${{[\![}} v{{]\!]}}_{i+\frac12} = v(x_{i+\frac12}^+)-v(x_{i+\frac12}^-)$.

For $v, w\in V_h$, we introduce some notations,

For any $v = \{v_0,v_b\}\in V_h$, a weak derivative $D_w f(v)\in P_{k}(I_i)$ for $i = 1,\cdots,N$ satisfies

where $n = -1$ at $x_{i-\frac12}$ and $n = 1$ at $x_{i+\frac12}$. Here $Dw$ means the first order derivative of $w$, i.e. $Dw = w'$.

Then we introduce a stabilizer in $V_h$ as follows:

where $\partial_+I_i = x_{i+\frac12}$, $\partial_-I_i = x_{i-\frac12}$, and $\lambda_1$ and $\lambda_2$ are two parameters. Here $v_0(x_{i+\frac12})$ and $w_0(x_{i+\frac12})$ are the left limit of $v_0(x)$ and $w_0(x)$ at $x_{i+\frac12}$ respectively, and $v_0(x_{i-\frac12})$ and $w_0(x_{i-\frac12})$ are the right limit of $v_0(x)$ and $w_0(x)$ at $x_{i-\frac12}$ respectively.

The following is the semi-discretized weak Galerkin method.

Algorithm 1 (SD-WG method). A numerical approximation for (1.1)-(1.2) can be obtained by seeking $u_h(t) = \{u_0(t),u_b(t)\}\in V_h$ satisfying $u_h(0) = Q_0\phi$ and the following equation,

where $Q_0$ is the $L^2$ projection onto $P_k(I_i)$ on each element $I_i$.

Let $I_i$ and $I_{i+1}$ be the two intervals sharing $x_{i+\frac12}$. Define the average ${{\{\{}} v {{\}\}}}$ on $x_{i+\frac12}$ by

and

Please note that the definition of average ${{\{\{}} \cdot{{\}\}}}$ above is different from the standard definition of average which is when $\lambda_1 = \lambda_2$.

Testing (2.4) by $v = \{v_0, v_b\}$ such that $v_0 = 0$ and $v_b = 1$ at $x_{i+\frac12}$ and $v_b = 0$ otherwise, we can easily obtain that at $x_{i+\frac12}$

Remark 1 (Relation to the upwinding-type DG method). If $f'(u)>0$ and taking $\lambda_2 = 0$, then $u_b = u_0^{-}$ and $s_h(u_h, v_h) = 0$. Thus, the WG scheme (2.4) reduces to

which is the classical upwinding type discontinuous Galerkin method for nonlinear conservation law.

Remark 2 (Relation to the upwinding-biased DG method). Considering a special case $f(u) = u$. Taking $\lambda_1 = \lambda_2 = \lambda$ in (2.3), and denote by

then the WG scheme (2.4) reduces to

which is the upwinding-biased DG method discussed in [11].

2.2.   Stability of the semi-discrete WG scheme

In this subsection, we will study the stability of the semi-discrete WG scheme (2.4).

Lemma 2.1. Let $\alpha_2 \le f'(s)\le \alpha_1$. If $\lambda_1> \alpha_1/2$ and $\lambda_2> -\alpha_2/2$, then for $v = \{v_0,v_b\}\in V_h$, there holds

Proof. As [5], we introduce $g'(s) = f(s)$, then

Since $f(v_b)$ and $v_b$ take single value on $\partial I_i$, the periodic boundary condition implies

Using the definition of the weak derivative (2.2), the mean value theory and the periodic boundary condition that

which implies

and

where $\alpha = \max \{|\alpha_1|,|\alpha_2|\}$. The equation (2.8) gives

We complete the proof.

Define $\|v\|^2 = \int_I v^2dx$. Then we have the following stability result.

Lemma 2.2 (Stability of the SD-WG method). Let $u_h = \{u_0,u_b\}\in V_h$ be the solution of the semi-discrete WG scheme (2.4), then

Proof. Let $v = u_h$ in the semi-discrete WG scheme (2.4). From (2.6), we have

Integrating the above inequality with respect to time between $0$ and $T$ completes the proof.

3.   Forward Euler WG method and its error analysis

We use forward Euler method for time discretization to obtain a full discrete WG finite element method.

3.1.   Forward Euler time discrete WG scheme

Let $\tau$ be a time step and $t_n = n\tau$.

Algorithm 2 (FE-WG method). Find $u_h^{n+1} = \{u_0^{n+1}, u_b^{n+1}\}\in V_h$ satisfying $u_h^0 = \{Q_0 \phi,{{\{\{}} Q_0\phi{{\}\}}}\}$ and

Testing (3.1) by $v = \{v_0, v_b\}$ such that $v_0 = 0$ and $v_b = 1$ at $x_{i+\frac12}$ and $v_b = 0$ otherwise, we can easily obtain that at $x_{i+\frac12}$

We define a projection operator $Q_hu = \{Q_0u, Q_bu\}\in V_h$, where $Q_0$ is the $L^2$ projection onto $P_k(I_i)$ on each element $I_i$ and $Q_bu = {{\{\{}} Q_0u{{\}\}}}$ on $\partial I_i$.

Define

and

Lemma 3.1. The error function $e_h^n$ defined in (3.3) satisfies the following equation,

Proof. Testing (1.1) by $v_0$ of $v = \{v_0,v_b\}\in V_h$ we arrive at

The definition of $D_w$ implies,

It follows from the definition of $Q_0$ and the above equation,

Adding $\tau s_h(Q_hu(t_n),v)$ to the both sides of the above equation implies

The difference of (3.1) and the equation above yields

Using the fact $2p(p-q) = p^2+(p-q)^2-q^2$ and letting $v = e_h^{n+1}$, (3.6) becomes

We have proved the lemma.

3.2.   Error analysis of forward Euler WG method

In this subsection we carry out an a priori error estimate for the fully discrete WG scheme with forward Euler time marching for smooth solutions. We will assume the nonlinear flux function $f(u)$ is smooth enough for simplicity.

For any function $\varphi\in H^1({ I_i})$, the following trace inequality holds true,

Lemma 3.2. Let $\rho_h^n$ and $e_h^n$ be defined in (3.4) and (3.3) respectively. Then we have

Proof. The first estimate in (3.8) is a direct result of the approximation property of the $L^2$ projection $Q_0$. The second estimate in (3.8) follows from the trace inequality (3.7) and the definitions of $\rho_b^n$ and $Q_h$,

Similarly, we can prove (3.9).

Lemma 3.3. Let $\tau\le ch^2$. Then we have

Proof. It follows from (3.7), the Cauchy-Schwarz inequality and the definition of $Q_b$ that

Using the assumption $\tau\le ch^2$, we have

which proves (3.10).

It follows from the Cauchy-Schwarz inequality that

Using (3.7), Cauchy-Schwarz inequality, Lemma 3.2 and the assumption $\tau \leq ch^2$, we have

We have proved the lemma.

Lemma 3.4. For $\tau\le ch^2$ and $\lambda_1,\lambda_2 \ge \alpha$. Then we have

Proof. It follows from definition of the weak derivative $D_w$ and the Taylor theory,

Next, we will bound all the terms above. Define $\overline{f'(u)} = \frac{1}{|I_i|}\int_{I_i}f'(u)dx$. Using the definitions of $\rho_0^n$ and $Q_0$, (3.8) and the inverse inequality, we arrive

It follows from the fact ${{\langle}} f'(u(t_n))\rho_b^n, e_b^nn{{\rangle}}_{\partial{{\mathcal T}}_h} = 0$, (3.7) and the inverse inequality,

Using the inverse inequality, the trace inequality (3.7), (3.9) and (2.9), we obtain

As [22,23], we adopt the following a priori assumption to deal with the nonlinearity of the flux $f(u)$

And the justification of such assumption will be given in Remark 3. Then the above estimate with (3.2) gives

Using (3.8), (3.15) and the inverse inequality, we have

Combining all the estimates above gives

which proves the lemma.

Theorem 3.5. Let $u_h^{n+1}\in V_h$ be the WG finite element solution of the problem (1.1)-(1.2) arising from (3.1). Then there exists a constant $C$ such that

Proof. It follows from (3.5) that

Then using Lemma 3.3 and Lemma 3.4, we have

Letting $\epsilon_1$ and $\epsilon_2$ small enough gives

Summing the above equation over $n+1$, we have

The discrete Gronwall's inequality implies

which proves the theorem.

Remark 3. The assumption (3.14) is obviously satisfied for $n = 0$ since $u_0^0 = Q_0\phi$. If (3.14) holds for a certain $n$, then it follows from the conclusion of Theorem 3.5 and $\tau\leq ch^2$ that

for $k\geq 1$. Thus the given a priori (3.14) is verified.

The result in the following theorem is a special case of Theorem 3.5. We omit the proof of the theorem since it is similar to the proof of Theorem 3.5.

Theorem 3.6. Let $u_h^{n+1}\in V_h$ be the WG finite element solution of the problem (1.1)-(1.2) arising from (3.1). Assume $f'(u)>0$ (or $f'(u)<0$) and let $\lambda_2 = 0$ (or $\lambda_1 = 0$). Then there exists a constant $C$ such that

4.   TVDRK3 WG method and its error analysis

This section discusses the fully discrete WG method coupled with the explicit TVDRK3 time-marching.

First, we set the initial value $u_h^0 = \{Q_0 \phi,{{\{\{}} Q_0\phi{{\}\}}}\}$. Then for each $n\geq 0$, the approximate solution from the time $n\tau$ to the next time $(n+1)\tau$ is defined as follows:

Algorithm 3 (RK3-WG method). Find $u_h^{n,1} = \{u_0^{n,1},u_b^{n,1}\}, u_h^{n,2} = \{u_0^{n,2},u_b^{n,2}\}$ and $u_h^{n+1} = \{u_0^{n+1},u_b^{n+1}\}$ in the finite element space $V_h$ such that, for any $v\in V_h$

Similar as (3.2), we obtain the time update for $u_b$ as

Following [22], two reference functions are defined in parallel to the TVDRK3 time discretization stages for the exact solution of the conservation law (1.1). Let $u^{(0)}(x,t) = u(x,t)$ and

Lemma 4.1 ([23]). If $\|u_{tttt}\|$ is bounded uniformly for any $t\in [0,T]$, we have

where $\mathcal{E}(x,t)$ is the local truncation error in time and $\|\mathcal{E}(x,t)\| = \mathcal{O}(\tau^4)$ uniformly for any time $t\in [0,T]$.

Denote $u^{n,i} = u^{(i)}(x,t^n)$ for any time level $n$ and $i = 0,1,2$. The error at each stage is denote by

for any $n$ and inner state $i = 0,1,2$, where $u_h^{n,0} = u_h^n$.

Lemma 4.2 (Error equation). For any $v\in V_h$, the error functions $e_h^{n,1}, e_h^{n,2}$ and $e_h^{n+1}$ satisfy the following equations

where

Proof. Let $t = t^n$ and test (4.5) by $v_0$ of $v = \{v_0, v_b\}$ in $V_h$, we arrive at

The definition of $D_w$ implies

It follows from the definition of $Q_0$ and the above equation that

The difference of (4.1) and the equation above yields (4.8).

The error equations (4.9) and (4.10) can be obtained similarly. The details are therefore omitted.

4.1.   Some basic estimates

For any $u_h = \{u_0,u_b\}$ and $v_h = \{v_0, v_b\}$ in $V_h$, it follows from the definition of weak derivative $D_w$ that

where

Since $\mathcal{R}_{b}(u,u_0)$ and $f'(u)(u-u_b)$ take single value on $\partial I_i$, the periodic boundary condition implies

Therefore, the operator $\mathcal{H}(u,u_h;v_h)$ defined in (4.11) can be rewritten as

for $u_h, v_h\in V_h$, where

Lemma 4.3. For any continuous and differentiable function $\mathcal{Z}$, there exists a positive constant $C$, independent of $n,h,\tau$, and $u_h$, such that

where $\varepsilon$ is any positive constant, and $|\mathcal{Z}|{{[\![}} v{{]\!]}}^2 = \sum_{j = 1}^N|\mathcal{Z}_{j+1/2}|{{[\![}} v{{]\!]}}_{j+1/2}^2$.

Proof. For the proof of (4.13), it follows from Cauchy-Schwartz inequality and the inverse inequality that

A simple manipulation indicates that

which implies the second conclusion (4.14).

Let $\overline{\mathcal{Z}}_j = \frac{1}{|I_j|}\int_{I_j}\mathcal{Z}\mathrm{d}x$. It follows from Cauchy-Schwartz inequality and Young's inequality that

The proof is completed.

Lemma 4.4. Let $u\in H^{k+1}(I)$ be the exact solution of (1.1). For any $u_h = \{u_0,{{\{\{}} u_0{{\}\}}}\}, v_h = \{v_0,{{\{\{}} v_0{{\}\}}}\}\in V_h$, there holds

Proof. By Taylor expansion up to the second order derivative term, we obtain

where $f''(\xi_1)$ and $f''(\xi_2)$ are the second order derivative of $f$ in the two expansion, which are both bounded.

Thus, by the triangle inequality and the inverse inequality, we have

which together with Young's inequality completes the proof.

Lemma 4.5. Let $u\in H^{k+1}(I)$ be the exact solution of (1.1). For any $u_h = \{u_0,{{\{\{}} u_0{{\}\}}}\}, v_h = \{v_0,{{\{\{}} v_0{{\}\}}}\}\in V_h$, there { hold}

where $\varepsilon$ is any positive constant, ${{[\![}} e_{0}{{]\!]}}^2: = \sum_{j = 1}^N{{[\![}} e_{0}{{]\!]}}^2_{ j+\frac{1}{2}}$, and $\vartheta = \frac{\lambda_1\lambda_2}{\lambda_1+\lambda_2}$.

Proof. Since $u$ is continuous function, there holds ${{[\![}} u{{]\!]}} = 0$. It follows from (2.3) and (4.4) that

Hence, by the triangle inequality and the inverse inequality, we have

which completed the proof of (4.16).

From (4.18) and Young's inequality, we have

The proof is completed.

Lemma 4.6. If the time step satisfies $\tau = \mathcal{O}(h)$, then we have

where $C$ is a positive constant independent of $n,h,\tau$ and $u_h$.

Proof. Taking $v = e_h^{n,1}$ in the error equation (4.8), we get

Firstly, we consider the first term of $\mathcal{H}(u^n, u_h^n;e_h^{n,1})$. It follows from the definition of the weak derivative $D_w$ that

From Cauchy-Schwarz inequality, the inverse inequality and the approximation property of $Q_0$, we have

where the assumption $f'(u)$ is bounded has been used in the second inequality above.

Now we turn to the term $T_2$. Since periodic boundary condition is considered, $T_2$ can be rewritten as

Then, by the triangle inequality and the inverse inequality, we obtain

Collecting (4.22), (4.23) and (4.24), we can conclude

which together with the inequality (4.16) of Lemma 4.5 and (4.21) yields

Cancelling $\|e_0^{n,1}\|$ on both sides of the above equation, and noting that $\tau = \mathcal{O}(h)$, we have

which implies (4.19).

In a similar way, we can obtain the conclusion (4.20). The proof is completed.

By taking the test function $v = e_h^n,4e_h^{n,1}$, and $6e_h^{n,2}$ in the error equations (4.8), (4.9) and (4.10), respectively, we obtain the energy equation for $e_h^n$ in the form

where

4.2.   The estimate for $\Xi_1$

It follows from (4.14) and (4.15) of Lemma 4.3 that

Denote by $\mathbb{B} = \max_{s\in\mathbb{R}} |f'(s)|$. In the above inequality, taking $\varepsilon$ small enough such that $\varepsilon \mathbb{B}\leq 0.1\vartheta$, then we have

From Lemma 4.4, we have

Using (4.17) of Lemma 4.5 with $\varepsilon = 0.1$, we have

Denote by

where $C$ is a general positive constant independent of $h,\tau,u$ and $u_h$. Collecting (4.26), (4.27) and (4.28), we obtain

Similarly, we have

for $i = 1,2$.

Since $\mathcal{E} = \mathcal{O}(\tau^4)$, by Cauchy-Schwartz inequality, we have

Collecting the three above estimates, we obtain the estimate of $\Xi_1$ as follows

4.3.   The estimate for $\Xi_2$

Following [23], we introduce the following notations

Obviously,

From the above equalities and the error equations (4.8)-(4.10), we easily obtain the following results.

Lemma 4.7. For the fully discrete WG method (4.1)-(4.3) with the explicit TVDRK3 time marching, we have the following equations

for any $v_h = \{v_0, v_b\}\in V_h$.

Note that $e_0^{n+1}-e_0^n = E_0^n+F_0^n+G_0^n$, there holds

Let $v_h = F_h^n$ in (4.31) and $v_h = E_h^n$ in (4.32), and by (4.12), we have

where

Lemma 4.8 (Estimate of $\mathcal{M}_1$). If the time step satisfies $\tau = \mathcal{O}(h)$, then we have

Proof. Since $u^{n,i}-u_h^{n,i} = \rho_h^{n,i}+e_h^{n,i} (i = 0,1,2)$, then $\mathcal{M}_1$ can be rewritten as

where

for $w = \rho_h$ or $e_h$.

Denote by

and collecting the terms with the same speed $f'(u^n)$, the operator $\mathcal{L}(w)$ can be expressed as

Now we estimate $\mathcal{L}_i(e_h), i = 1,\cdots,5$ firstly. Due to periodic boundary condition, and using integration by parts, we have

Consequently, by Cauchy-Schwartz inequality, we can conclude that

Since $|z^{n,1}| = \mathcal{O}(\tau) = \mathcal{O}(h)$, it follows from the equation (4.13) of Lemma 4.3 that

In a similar way, we obtain

Collecting the estimate of $\mathcal{L}_i(e_h), i = 1,\cdots,5$ and using Cauchy-Schwarz inequality yields

Next, we turn to the term $\mathcal{L}_i(\rho_h), i = 1,\cdots,5$. Note that $f'(u)$ is bounded, it follows from the Cauchy-Schwartz inequality and inverse equality that

Since

it follows from the approximation property of $L^2$ projector $Q_0$ that

which together with (4.37) yields

Similarly we can estimate the second term as $|\mathcal{L}_2(\rho_h)|\leq C(\|E_0^n\|^2+h^{2k}\tau^2)$.

Using $|z^{n,i}| = \mathcal{O}(\tau) = \mathcal{O}(h)$ and similar as the proof of $\mathcal{L}_1(\rho_h)$, it is easy to obtain

Since $\tau = \mathcal{O}(h)$, finally we have

By the triangle inequality, we have

which together with (4.36) and (4.38) completes the proof.

Lemma 4.9 (Estimate of $\mathcal{M}_2$). There exists a positive constant $C$, independent of $n,h,\tau$, and $u_h$, such that

Proof. It follows from Lemma 4.4 that

for $w_h\in V_h$ and $i = 0,1,2$.

Let $w_h = E_h^n$ and $F_h^n$ in the above inequality, respectively. By the triangle inequality and the definition of $\mathcal{M}_2$, we have

which completes the proof.

Lemma 4.10 (Estimate of $\mathcal{M}_3$). There exists a positive constant $C$, independent of $n,h,\tau$, and $u_h$, such that

Proof. Using (4.17) of Lemma 4.5 with $\varepsilon = 0.1$, we have

Combining the above two estimates completes the proof.

Since $\mathcal{E} = \mathcal{O}(\tau^4)$, by Cauchy-Schwartz inequality, we have

which together with (4.34) and Lemma 4.8, 4.9 and 4.10 yields

Now we consider the upper bound of $J_3$ and $J_4$. Since $J_3 = 3(G_0^n, F_0^n)$ and $J_4 = 3(G_0^n,G_0^n)$, we estimate $3(G_0^n, v_0)$ in general. From (4.12) and (4.32), we have

where

Denote by

and collecting the terms with the same speed $f'(u^n)$, we rewrite the operator $\mathcal{I}_1(v_0)$ as

It follows from (4.13) of Lemma 4.3, Cauchy-Schwartz inequality and $|z^{n,i}| = \mathcal{O}(\tau) = \mathcal{O}(h)$ that

Note that $f'(u)$ is bounded, it follows from the Cauchy-Schwartz inequality and inverse equality that

Since

and

by the approximation property of $L^2$ projector $Q_0$, we have

Plugging (4.45) into (4.44), and using $2ab\leq a^2+b^2$, we get

Similar to the proof of (4.44), and note that $|z^{n,i}| = \mathcal{O}(\tau)$, we have

Combining (4.41), (4.42), (4.43), (4.46) and (4.47), and assume $\tau = \mathcal{O}(h)$, we have

Thanks to Lemma 4.4, we have

It follows from (4.18) and Young's inequality that

which together with (4.45) and inverse inequality yields

Using the fact $\mathcal{E} = \mathcal{O}(\tau^4)$ and Cauchy-Schwartz inequality, we get

which combining with (4.48), (4.49), (4.50) and (4.40) yields

Taking $v_0 = F_0^n$ in (4.51) and using the assumption $\tau \leq \gamma h$, we have

Taking $v_0 = G_0^n$ in (4.51), we have

then

Therefore,

which combining with (4.52) yields

In the above inequality, taking $\tau$ small enough such that $C\gamma = 1/4$, then we get

which together with (4.39) and (4.33) yields

4.4.   Main results

In this subsection, we will give the $L^2$ norm error estimate between the exact solution $u(t^n)$ and the WG solution $u_0^n$.

Theorem 4.11. Let $u_h$ be the numerical solution of the fully discrete WG scheme (4.1)-(4.4) with the explicit TVDRK3 time marching. Let $u$ be the exact solution of problem (1.1)-(1.2), where $f(u)$ is smooth enough. If $u$ and its spatial derivatives up to the second order are all continuous in $I = (0,1)$, and $\|u\|_{k+1}$, $\|u_t\|_{k+1}$, $\|u_{tt}\|_{k+1}$ and $\|u_{ttt}\|_{k+1}$ are bounded uniformly for any time $t\in (0,T]$, then the following error estimate holds

under the condition $\tau\leq \gamma h$ with a fixed constant $\gamma >0$. Here $C$ is a positive constant independent of $h,\tau$ and $u_h$.

Proof. Denote by

then from (4.25), (4.29) and (4.53), we obtain

As [22,23], we adopt the following a priori assumption for $m\, (m\tau<T)$ to deal with the nonlinearity of the flux $f(u)$

The justification of such assumption will be given later.

It follows from the assumption (4.56) and the approximation property of $L^2$ operator $Q_0$

implies that

Therefore, we have

which together with (4.55), (4.19) and (4.20) yields

where $C$ is a positive constant independent of $m$, $n$, $h$ and $\tau$. Then an application of the discrete Gronwall inequality yields

which together with

yields that

Now we turn to verify the reasonableness of the a priori assumption (4.56). Since $e_0^0 = 0$, by (4.19) and (4.20), we easily obtain

for $i = 0,1,2$ and $k\geq 1$. Supposing (4.56) holds for $m$, we can show that this assumption is also true for $m+1$. Inequality (4.57) as well as (4.19) and (4.20) imply that

for $k\geq 1$ and $i = 1,2$. Thus the assumption (4.56) is reasonable, and hence the above error estimate holds for any $m\,(m\tau<T)$. The proof is completed.

5.   Numerical examples

In this section, we present the numerical examples to verify our theoretical findings. In our numerical experiments, we shall use piecewise uniform meshes which are constructed by equally dividing spatial domain into $N$ subintervals. The main purpose of this paper is to investigate how to use WG method to discretize spacial variables in conservation law. In order to reduce the time errors, we use the third order explicit TVDRK time discrete scheme in time, and take time step $\tau = CFL\cdot h$, where $CFL$ is a constant independent of the mesh size $h$. WG scheme with $P_k (k = 1,2,3)$ element is used for spatial discretization.

Example 1 (A smooth solution of a linear equation). We solve the following linear problem (1.1) with periodic boundary condition,

The exact solution to this problem is

Let time step size $\tau = 0.05h$ in this example. WG scheme with $P_k (k = 1,2,3)$ element is used for spatial discretization. The $L^2$ errors and the order of convergence, at $T = 2\pi$, are reported in Table 1. It is observed that the order of convergence of the $L^2$ error achieves $(k+1)$-th order of accuracy.

Example 2 (A blow-up solution of a nonlinear equation). We solve the following nonlinear problem (1.1) with periodic boundary condition,

with initial value function

The exact solution $u(x,t)$ is obtained by solving characteristic relation $u = \phi(x-u\cdot t)$, which is smooth up to $t = 1/\pi$. Please find more details in [4].

Let time step size $\tau = 0.1h$ in this example. In Table 2 we present the results at $T = 0.2$ when the solution is still smooth. Observe that there is at least $(k+1/2)$-th order of convergence rate in $L^2$ norm. $u_b$ of WG solution $u_h = \{u_0,u_b\}$ at $T = 1/\pi$ is plotted in Figure 1, where $P_2$ element is used.

Example 3 (A discontinuous solution of a linear equation). We solve the following linear equation (1.1) with periodic boundary condition,

with initial value function

The exact solution is $u(x,t) = \phi(x-t)$.

In this problem the exact solution is non-smooth. It aims to show our WG method is also stable and efficient to the problems with non-smooth solutions. Let time step size $\tau = 0.1h$ in this example. The errors and the order of convergence are reported in Table 3. Observe that the method is stable. The WG solution using $P_1$ element with $(\lambda_1,\lambda_2) = (1.2,0.1)$ and the DG solution with $P_1$ element are plotted in Figure 2 (a) and (b) respectively.

6.   Conclusion

A weak Galerkin finite element method is proposed for nonlinear conservation laws. This method provides a class of finite element schemes by tuning the built-in parameters including purely upwind scheme. This makes the WG method highly flexible. Compared with the DG method [11], the present method can be applied for solving nonlinear conservation laws. Error estimations are given. The convenience of the proposed method is validated by the numerical examples.