Eﬃcient time discretization for local discontinuous Galerkin methods

In this paper, we explore a few eﬃcient time discretization techniques for the local discontinuous Galerkin (LDG) methods to solve partial diﬀerential equations (PDEs) with higher order spatial derivatives. The main diﬃculty is the stiﬀness of the LDG spatial discretization operator, which would require a unreasonably small time step for an explicit local time stepping method. We focus our discussion on the semi-implicit spectral deferred correction (SDC) method, and study its stability and accuracy when coupled with the LDG spatial discretization. We also discuss two other time discretization techniques, namely the additive Runge-Kutta (ARK) method and the exponential time diﬀerencing (ETD) method, coupled with the LDG spatial discretization. A comparison is made among these three time discretization techniques, to conclude that all three methods are eﬃcient when coupled with the LDG spatial discretization for solving PDEs containing higher order spatial derivatives.


Introduction
In this paper we are concerned with efficient time discretization of the local discontinuous Galerkin (LDG) method for solving time dependent partial differential equations containing higher order spatial derivatives.
The discontinuous Galerkin (DG) method is a class of finite element methods, using discontinuous, piecewise polynomials as the solution and the test space. It was first designed as a method for solving hyperbolic conservation laws containing only first order spatial derivatives, e.g. Reed and Hill [23] for solving linear equations, and Cockburn et al. [8,7,6,9] for solving nonlinear equations.
The LDG method is an extension of the DG method aimed at solving partial differential equations (PDEs) containing higher than first order spatial derivatives. The first LDG method was constructed by Cockburn and Shu in [10] for solving nonlinear convection diffusion equations containing second order spatial derivatives. Their work was motivated by the successful numerical experiments of Bassi and Rebay [4] for the compressible Navier-Stokes equations. Yan and Shu developed an LDG method for a general KdV type equation (containing third order spatial derivatives) in [29], and they generalized the LDG method to PDEs with fourth and fifth order spatial derivatives in [30]. Levy, Shu and Yan [21] developed LDG methods for nonlinear dispersive equations that have compactly supported traveling wave solutions, the so-called "compactons". More recently, Xu and Shu [25,26,27,28] further developed the LDG method to solve many nonlinear wave equations with higher order derivatives, including the general KdV-Burgers type equations, the general fifth-order KdV type equations, the fully nonlinear K(n, n, n) equations, the generalized nonlinear Schrödinger equations, the coupled nonlinear Schrödinger equations, the Kuramoto-Sivashinsky equations, the Ito-type coupled KdV equations, the Kadomtsev-Petviashvili equation, and the Zakharov-Kuznetsov equation. A common feature of these LDG methods is that stability can be proved for quite general nonlinear cases. DG and LDG methods also have several attractive properties, such as their flexibility for arbitrary h and p adaptivity and their excellent parallel efficiency.
The application of the DG or LDG method discretizes the spatial variables and generates a large coupled system of ordinary differential equations (ODEs). One would then need to use a suitable ODE solver to discretize the time variable. For hyperbolic conservation laws and convection dominated convection diffusion equations, explicit and nonlinearly stable Runge-Kutta time discretizations [24,16] are suitable choices. The resulting fully discretized scheme, termed Runge-Kutta DG (RKDG) method [8], are stable, efficient and accurate for solving such convection or convection dominated problems. However, for PDEs containing higher order spatial derivatives, especially when the coefficients in front of these higher order spatial derivative terms are not small, such explicit and local time discretization will suffer from extremely small time step restriction for stability, due to the stiffness of the spatial LDG operator which approximates both the lower and higher order spatial derivatives. Often, such a small time step is not needed for the purpose of accuracy and is purely an artifact of the explicit time discretization technique. It would therefore be desirable to use either nonlocal or implicit time discretization techniques to alleviate this problem. Implicit time discretization is already used in [30] for some of the numerical experiments. In [27,28], the exponential time differencing (ETD) method [13,18], which is an explicit but non-local time discretization method, is used for some of the test cases. A suitable time discretization technique for the LDG methods solving PDEs with higher order spatial derivatives should have good stiffness stability and should be easy to implement. Such methods are usually implicit [17].
The ODE system obtained from the LDG spatial discretization for a PDE containing higher order spatial derivatives often includes both stiff and non-stiff terms due to a disparity in time scales. An efficient time discretization technique is a semi-implicit method treating the non-stiff terms (often low-order but nonlinear) explicitly and the stiff terms (often higher order but linear) implicitly. This is much easier to implement and much less in computational cost than the fully implicit time discretization methods. If only the convection terms are nonlinear, the semi-implicit methods are as effective as the fully implicit methods in allowing large time steps. The main time discretization method explored in this paper is the spectral deferred correction (SDC) method constructed by Dutt, Greengard and Rokhlin [15]. An advantage of this method is that it is a one step method (namely, to march to time level n+1 one would only need to store the value of the solution at time level n) and can be constructed easily and systematically for any order of accuracy. This is in contrast to Runge-Kutta methods which are more difficult to construct for higher order of accuracy. The SDC method relies on the Picard integral equation and is driven iteratively by the base Euler forward (for non-stiff problems) or Euler backward (for stiff problems) method. It has been shown in [15,22] that the implicit SDC method driven by the Euler backward method is A-stable or A(α)-stable, with α close to 90 • . In [15] an L-stable implicit SDC method is also constructed, which is a combination of two different SDC methods and hence doubles the cost of the usual SDC method. We present in this paper a different way to construct an L-stable implicit SDC method with a smaller computational cost than that in [15].
Two other techniques for time discretization of the LDG methods are also explored in this paper. One is the additive Runge-Kutta method given in [19], which has nice stiffness stability. The other is the exponential time differencing (ETD) method given in [13,18], which is an explicit but nonlocal time discretization. We perform numerical experiments to compare these three different techniques of time discretization, and conclude that they are all competitive when solving ODEs from the LDG discretization of PDEs containing higher order spatial derivatives.
We remark that, for stiff problems, implicit Runge-Kutta methods possess excellent stability properties. But they are very expensive when high order accuracy is required except for the diagonally implicit Runge-Kutta methods [1,2,20,17]. The additive Runge-Kutta method [19] that we explore in this paper belongs to this class. Implicit multi-step algorithms can be easily constructed for very high order of accuracy, but they have poor stability properties.
The organization of the paper is as follows. In Section 2, we give a description of the spectral deferred correction (SDC) method. In Section 3, we give a brief introduction to the additive Runge-Kutta (ARK) and exponential time differencing (ETD) methods. Section 4 is devoted to a concise description of the local discontinuous Galerkin (LDG) spatial discretization via an example. Numerical examples are presented in Section 5, testing the performance of the three different time discretization techniques with the LDG spatial discretization for a series of PDEs with higher order spatial derivatives including the Korteweg-de Vries (KdV) equation, the Kuramoto-Sivashinsky equation and the Kawahara equation. Finally, we give concluding remarks in Section 6.
2 The spectral deferred correction method Dutt, Greengard and Rokhlin presented a new variation of the classical method of deferred correction, called the spectral deferred correction (SDC) method, in [15]. It is based on low order time integration methods, which are corrected iteratively, with the order of accuracy increased by one for each additional iteration. Attractively, we can split stiff and non-stiff terms as needed and treat them differently (implicitly for the stiff terms and explicitly for the non-stiff terms). Minion presented the semi-implicit SDC (SISDC) method in [22]. In the following we will introduce these methods. Consider the ODE system where u 0 , u(t) ∈ C n and F : R×C n −→ C n . F ∈ C 1 (R×C n ), which is sufficient to guarantee local existence and uniqueness of the solution to (2.1). Suppose now the time interval [0, T ] is divided into N intervals by the partition 0 = t 0 < t 1 < · · · < t n < · · · < t N = T . Let ∆t n ≡ t n+1 − t n and u n denotes the numerical approximation of u(t n ), with u 0 = u(0).
Rewrite (2.1) into an integral form in the subinterval [t n , t n+1 ]: Then divide the time interval [t n , t n+1 ] into P subintervals by choosing the points t n,m for m = 0, 1, · · · , P such that t n = t n,0 < t n,1 < · · · < t n,m < · · · < t n,P = t n+1 . Let ∆t n,m = t n,m+1 −t n,m and u k n,m denotes the k th order approximation to u(t n,m ). To avoid the instability of approximation at equispaced nodes for high order accuracy, the points {t n,m } P m=0 are chosen to be the Gauss-Lobatto nodes on [t n , t n+1 ]. We can also use the Gauss-Legendre, Gauss-Radau or Chebyshev nodes. Starting from u n , we give the algorithm to calculate u n+1 in the following.
Compute the initial approximation u 1 n,0 = u n . For non-stiff/stiff problems, use the forward/backward Euler method to compute a first order accurate approximate solution u 1 at the nodes {t n,m } P m=1 . For m = 0, · · · , P − 1 1. For non-stiff problems, 2. For stiff problems, Compute successive corrections For k = 1, · · · , K u k+1 n,0 = u n . For m = 0, · · · , P − 1 is the integral of the P -th degree interpolating polynomial on the P +1 points (t n,m , F (t n,m , u k n,m )) P m=0 over the subinterval [t n,m , t n,m+1 ], which is the numerical quadrature approximation of F (τ, u(τ ))dτ. (2.7) Finally we have u n+1 = u K+1 n,P .

Lemma 2.1 (Local Truncation Error). The Local Truncation Error (LTE) obtained with the explicit or the implicit SDC
Proof: First, we consider the LTE of the explicit SDC method. Assume u k n,0 = u(t n ), k = 1, · · · , K + 1. Taking the difference of and (2.5), we have By induction on both m and k, we assume, when k ≤ P , Then where we have used the fact that I m+1 m (F (t, u)) is the integral of the P -th degree interpolating polynomial on (t n,m , F (u(t n,m )) P m=0 over the subinterval [t n,m , t n,m+1 ], hence it is accurate to the order O(h P +2 ). Thus we have The proof is the same for the implicit SDC K P (θ) method.
Remark 2.1. In the proof of Lemma 2.1, we find the terms θ∆t m (F (t n,m , u k+1 m )−F (t n,m , u k m )) in (2.5) and θ∆t m (F (t n,m+1 , u k+1 m+1 ) − F (t n,m+1 , u k m+1 )) in (2.6) do not affect the accuracy. They are added purely for the purpose of enhancing stability. When θ = 1, the methods are the original ones presented in [15,22].
, through the following Gauss-Lobatto quadrature where we have used the fact that the Gauss-Lobatto quadrature over the whole interval [t n , t n+1 ] is accurate of order O(h 2P +1 ). Therefore, the SDC P P (θ) scheme coupled with (2.17) has the LTE one order higher O(h P +2 ). From Remark 2.1, the usage of (2.17) is equivalent to adding one additional correction step k = P + 1. We denote this version of the scheme by SDC P +1 P (θ), which has the LTE O(h P +2 ) when P ≥ 2. This is the version of the scheme actually used in the numerical experiments in Section 5.
From Lemma 2.1, Remark 2.2 and the Lax stability for one-step time integrated methods, we have Theorem 2.2. The explicit or implicit SDC K P (θ) methods are min(K + 1, P + 2) order accurate methods for P ≥ 2. When P = 1, the SDC K P (θ) methods are min(K + 1, P + 1) order accurate.
When the right hand side F (t, u) of the ODE (2.1) can be written as the sum of a non-stiff term F N (t, u) and a stiff term F S (t, u), we have The construction of the semi-implicit SDC (SISDC) method is simply where θ 1 and θ 2 are taken to be the same number θ in our numerical experiments in Section 5. Generally, the numerical methods for stiff problems are analyzed by applying it to the equation where λ ∈ C and ∆t n = 1. The amplification factor, Am(λ), is defined by the formula where u 1 is the numerical solution of u(1). If |Am(λ)| ≤ 1, then the numerical method is said to be stable for the given value of λ. If the method is stable for all λ in the left-half plane We have plotted the stability regions of the implicit schemes SDC K P (θ) for several choices of the parameters P and K (SDC 1 , SDC 4 4 (θ) and SDC 5 4 (θ), with different values of θ (θ = 0.8, 0.9 and 1.0). These schemes are A-stable, but not L-stable.
As in [15], we have the following theorem.
Theorem 2.3. For any pair of natural numbers P , K, the amplification factor Am(λ) associated with the implicit scheme SDC K P (θ) is a rational function of λ. Moreover, there exists a real number µ(P, K, θ) such that lim λ→−∞ Am(λ) = µ(P, K, θ). (

2.24)
This theorem provides a mechanism for combining two different schemes with different P , K and θ to obtain L-stable schemes [15].
When the signs of µ(P 1 , K 1 , θ 1 ) and µ(P 2 , K 2 , θ 2 ) are opposite, Remark 2.3. In the implicit scheme SDC K P (θ), if we change θ to θ 1 only for one sub-interval in one of the correction steps, e.g. in the last step k = K, m = P − 1 u K+1 n,P = u K+1 n,P −1 + θ 1 ∆t n,P −1 (F (t n,P , u K+1 n,P ) − F (t n,P , u K n,P )) + I P P −1 (F (t, u K )) (2.26) then the corresponding µ(P, K, θ, θ 1 ) of this modified scheme is different from µ(P, K, θ). If we combine the original SDC K P (θ) scheme with this modified scheme as in Corollary 2.4, we obtain a new scheme, denoted by SDC K P (θ, θ 1 ), which is L-stable but involves much less computational cost than (2.25). In fact, (2.25) involves twice the cost of the original scheme while SDC K P (θ, θ 1 ) involves only a factor 1 + 1 P (K+1) of the original cost.
In Section 5, we give a comparison among the L-stable SDC 3,3 2,3 (1, 1) and SDC 3 2 (1, 0.8) schemes, and the non-L-stable SDC 3 2 (1) scheme, for different time steps in Example 5.2, to demonstrate the advantage of L-stable schemes for larger time steps in reducing the errors.
3 The additive Runge-Kutta and exponential time differencing methods In this section, we introduce the additive Runge-Kutta (ARK) and exponential time differencing (ETD) methods.
Following the work [3,19], the ARK methods are used to solve equation (2.19). They are given in the following form ij F S (t n + c j ∆t, u (j) ), i = 1, · · · , s (3.1) where u (0) = u n and u (i) approximates u(t n + c i ∆t). The non-stiff and stiff terms are integrated by their own (s+1)-stage Runge-Kutta methods respectively. The Butcher coefficients a and c j are constrained by order of accuracy and stability considerations. In [19], the implicit-explicit additive Runge-Kutta (ARK) methods from third-to fifth-order are presented in which the stiff terms are integrated by an L-stable, stiffly-accurate, singly diagonally implicit Runge-Kutta method while the non-stiff terms are integrated with a traditional explicit Runge-Kutta method. For a detailed description of the methods as well as their implementation and applications, we refer the readers to [19].
The ETD method was constructed to solve the equation of the following form where L is a linear term and F is nonlinear. Many interesting equations are of this form, where typically L represents the stiff part of the equation. When discretizing in space we obtain a system of ODEs of the form The exponential time differencing (ETD) methods can be described in the context of solving (3.4). Integrating the equation over a single time step from t = t n to t n+1 = t n + ∆t, we get u(t n+1 ) = e L∆t u(t n ) + e L∆t ∆t 0 e −Lτ F (t n + τ, u(t n + τ ))dτ (3.5) The equation (3.5) is exact. We denote the numerical approximation to u(t n ) by u n and write F (t n , u n ) as F n . The simplest approximation to the integral (3.5) is to take F as a constant, F = F n + O(∆t), on the interval t n ≤ t ≤ t n+1 . Then we obtain the first order scheme ETD1, given by If instead of assuming that F is constant over the interval t n ≤ t ≤ t n+1 , we use the linear approximation that Then we obtain the second order scheme ETD2, given by u n+1 = e L∆t u n + 1 ∆t L −2 ((I + ∆tL)e L∆t − I − 2∆tL) F n + 1 ∆t L −2 (−e L∆t + I + ∆tL) F n−1 .

(3.8)
A second order ETD method of Runge-Kutta type, analogous to the ETD2 method, is as follows. First, the ETD1 (3.6) is taken to give a n = e L∆t u n + L −1 (e L∆t − I) F n . (3.9) Then the approximation is applied on the interval t n ≤ t ≤ t n+1 , and is substituted into (3.5) to yield the scheme ETD2RK given by u n+1 = a n + 1 ∆t L −2 (e L∆t − I − ∆tL)(F (t n + ∆t, a n ) − F n ). (3.11) Cox and Matthews derive a set of ETD methods bases on Runge-Kutta time-stepping, which they call ETDRK schemes in [13], among general formulas for ETD-schemes to arbitrary order. In a recent paper Kassam and Trefethen [18] compare various fourth order methods for solving equations of the form (3.4), and conclude that the best by a clear margin was a modification to the ETDRK schemes. In essence, for nonlinear time dependent equations, the ETD schemes provide a systematic coupling of the explicit treatment of nonlinearities and the implicit and possibly exact integration of the stiff linear parts of the equations, while achieving high accuracy and maintaining good stability.
To overcome the vulnerability of the error cancellations in the high order ETD and ETDRK schemes, and to generalize the ETD schemes to non-diagonal problems, in [18], modified ETD schemes are proposed by using the complex contour integrals on a suitable contour Γ to evaluate the coefficients (e.g. f (L) = L −2 (e L∆t − I − ∆tL) in (3.11)) in the update formula for ETDRK. We choose the unit circle in the complex space as Γ and the trapezoidal rule to approximate (3.12) in Section 5.

The local discontinuous Galerkin method
In this section, we give a simple example to illustrate the essential ideas of the local discontinuous Galerkin method [10]. Consider the following convection diffusion equation from the right cell, I j+1 , and from the left cell, I j , respectively. We define the space V: V = {v : v is a polynomial of degree at most k for x ∈ I j , j = 1, . . . , N} (4.2) as the solution and test function space. We rewrite the equation (4.1) into a first order system and apply the discontinuous Galerkin method to the system (4.3): find the u, q ∈ V, such that for all test functions v, z ∈ V, we have The "hat" terms are the boundary terms that emerge from integration by parts. Among them, is a monotone flux, i.e. Lipschitz continuous in both arguments, consistent ( f (u, u) = f (u)), and non-decreasing in the first argument and non-increasing in the second argument. We could, for example, use the simple Lax-Friedrichs flux where the maximum is taken over a relevant range of u. The fluxes q, u can be chosen as These "numerical fluxes" should be designed based on different guiding principles for different PDEs to ensure stability. For example, upwinding should be used as a guideline for odd derivatives which correspond to waves, and eventual symmetric treatment, such as an alternating choice of the fluxes for a quantity and its derivative as in (4.8) or (4.9), should be used for even derivatives. For a detailed description of the LDG method as well as its implementation and applications, we refer the readers to the lecture notes [5], the review paper [11], and papers in the two special issues on discontinuous Galerkin methods [12,14]. In the next section we will consider the linear KdV, linear bi-harmonic, linear fifth order derivative, KdV, Kuramoto-Sivashinsky, and Kawahara equations. The corresponding LDG methods can be found in [30,29,25,28].

Numerical results
In this section, we perform numerical experiments of the LDG scheme coupled with different time integration methods (SDC, ARK and ETD) to linear and nonlinear equations containing second to fifth order spatial derivatives. We use quadruple precision and the direct linear solver in LAPACK to compute these problems on uniform spatial meshes.
Example 5.1. We show an accuracy test for the linear KdV equation (5.1) in 1D: with the initial condition u(x, 0) = sin(x) and periodic boundary conditions. The exact solution is given by u(x, t) = sin(x + t). The time steps are taken as ∆t = ∆x. When the piecewise P 2 elements are used in the LDG method, the implicit SDC 2 2 (1) scheme and the implicit part of the third order ARK scheme [19] are used in the time integration. When the piecewise P 3 elements are used in the LDG method, we use the implicit SDC 3 2 (1) scheme and the implicit part of the fourth order ARK scheme. The numerical errors and orders of accuracy can be found in Table 5.1. Apparently the errors are comparable between these two different time discretizations.
Example 5.2. We show an accuracy test for the linear bi-harmonic equation (5.2) in 1D: with the initial condition u(x, 0) = sin(x) and periodic boundary conditions. The exact solution is given by u(x, t) = e −t sin(x + t). The time steps are taken as ∆t = ∆x. When the piecewise P 2 elements are used in the LDG method, the implicit SDC 2 2 (1) scheme and the implicit part of the third order ARK scheme [19] are used in the time integration. When the    piecewise P 3 elements are used in the LDG method, we use the implicit SDC 3 2 (1) scheme and the implicit part of the fourth order ARK scheme. The numerical errors and orders of accuracy can be found in Table 5.2. Apparently the errors are again comparable between these two different time discretizations.
We also test the performance of the L-stable SDC 3,3 2,3 (1, 1) and SDC 3 2 (1, 0.8) schemes in comparison with the non-L-stable SDC 3 2 (1) scheme for various time step sizes, coupled with the P 3 LDG method in the uniform mesh with N = 160 cells, in Table 5.3. We observe the advantage of the L-stable scheme for larger time step sizes in providing smaller errors. with the initial condition u(x, 0) = sin(x) and periodic boundary conditions. The exact solution is given by u(x, t) = sin(x − t). The time steps are taken as ∆t = ∆x. When the piecewise P 2 elements are used in the LDG method, the implicit SDC 2 2 (1) scheme and the implicit part of the third order ARK scheme [19] are used in the time integration. When the piecewise P 3 elements are used in the LDG method, we use the implicit SDC 3 2 (1) scheme and the implicit part of the fourth order ARK scheme. The numerical errors and orders of accuracy can be found in Table 5.4. Again the errors are comparable between these two different time discretizations.
Example 5.4. We show an accuracy test for the KdV equation with the initial condition u(x, 0) = −2sech 2 (x). The exact solution is given by The computational domain is given by [−20, 22], and we use a periodic boundary condition which is justified since the exact solution is negligibly small at the boundary of our computational domain. The time steps are taken as ∆t = C ∆x max x |6u(x,0)| , where C is the CFL number of the convection term associated with the stability of DG discretization (see [11]). When the piecewise P 2 elements are used in the LDG method, the implicit SDC 2 2 (1) scheme, the third order ARK scheme [19], and the third order ETD3RK scheme [18] are used in the time integration. When the piecewise P 3 elements are used in the LDG method, we use the implicit SDC 3 2 (1) scheme, the fourth order ARK scheme, and the fourth order ETD4RK scheme. For the SDC and ARK methods, the first order derivative term is treated explicitly (F N in (2.19)), and the third order derivative term is treated implicitly (F S in (2.19)). For the ETD method, the first order derivative term is treated as F (t, u) in (3.4), and the third order derivative term is treated as Lu in (3.4). The numerical errors and orders of accuracy can be found in Table 5.5. The errors are basically comparable among these three different time discretizations, with ETD showing slightly larger errors.
with the exact solution given by where we take σ = 4, c = 6, k = 1 2 and x 0 = −10. The computational domain is given by [−30, 30], and we use a periodic boundary condition which is justified since the exact solution is negligibly small at the boundary of our computational domain. The time steps are taken as ∆t = C ∆x max x |u(x,0)| , where C is again the CFL number of the convection term associated with the stability of DG discretization. When the piecewise P 2 elements are used in the LDG method, the implicit SDC 2 2 (1) scheme, the third order ARK scheme [19], and the third order ETD3RK scheme [18] are used in the time integration. When the piecewise P 3 elements are used in the LDG method, we use the implicit SDC 3 2 (1) scheme, the fourth order ARK scheme, and the fourth order ETD4RK scheme. For the SDC and ARK methods, the first order derivative term is treated explicitly (F N in (2.19)), and the second to fourth order derivative terms are treated implicitly (F S in (2.19)). For the ETD method, the first order derivative term is treated as F (t, u) in (3.4), and the second to fourth order derivative terms are treated as Lu in (3.4). The computational domain is given by [−35, 35], and we use a periodic boundary condition which is justified since the exact solution is negligibly small at the boundary of our computational domain. The time steps are taken as ∆t = C ∆x max x |u(x,0)| , where C is the CFL number of the convection term associated with the stability of DG discretization. When the piecewise P 2 elements are used in the LDG method, the implicit SDC 2 2 (1) scheme, the third order ARK scheme [19], and the third order ETD3RK scheme [18] are used in the time integration. When the piecewise P 3 elements are used in the LDG method, we use the implicit SDC 3 2 (1) scheme, the fourth order ARK scheme, and the fourth order ETD4RK scheme. For the SDC and ARK methods, the first order derivative term is treated explicitly (F N in (2.19)), and the third and fifth order derivative terms are treated implicitly (F S in (2.19)). For the ETD method, the first order derivative term is treated as F (t, u) in (3.4), and the third and fifth order derivative terms are treated as Lu in (3.4). Table 5.7 gives the errors of the numerical solution at t = 1. Again, the errors are basically comparable among these three different time discretizations, with ETD showing slightly larger errors.
For the sake of comparison of efficiency, we document the CPU time of the fourth order SDC 3 2 (1), ARK and ETD4RK methods coupled with the P 3 LDG method for solving the equation (5.8) for the same uniform mesh with N = 160 cells, and the same suitably chosen ∆t which yields an L ∞ error at the level of 10 −6 for all three methods. The machine we use is an Intel(R) Pentium(R) 4, with CPU 3.20GHz, Linux system with cache size 2048kb, memory is 3111120k total, and the compiler is the Intel(R) Fortran Compiler. The results are collected in Table 5.8. Comparing the CPU times, we notice that the ARK method requires less time than the SDC and ETD methods. Notice that we have not included the start-up cost (the CPU cost to evaluate the contour integral in the ETD method, which takes about 5070 seconds, and the CPU cost for the initial LU decomposition in the SDC and ARK methods, which takes about 3.7 seconds). If the number of function evaluations is compared, the ETD method requires fewer function evaluations than the ARK method, which in turn requires fewer function evaluations than the SDC method. How, we do see that all three methods are efficient in solving this highly stiff LDG problem.  We have explored three different time discretization techniques for solving the stiff ODEs resulting from a local discontinuous Galerkin (LDG) spatial discretization to PDEs containing higher order spatial derivatives. These are the semi-implicit spectral deferred correction (SDC) method, the additive Runge-Kutta (ARK) method and the exponential time differencing (ETD) method. Numerical experiments are performed to verify that all three methods are efficient in discretizing the LDG schemes in time. In particular, the SDC method has the advantage of easy implementation for arbitrary order of accuracy. In future work we will explore efficient linear and nonlinear solvers for discretizing LDG schemes for multi-dimensional problems with possible nonlinear higher order spatial derivative terms.