Reducing polynomial degree by one for inner-stage operators affects neither stability nor accuracy of the Runge--Kutta discontinuous Galerkin method
Zheng Sun























./figures/ 































definition
exampleExample[section]


 
 

tabularscale=1 
	tabular





*equationsection


Reducing polynomial degree for the RKDG method]Reducing polynomial degree by one for inner-stage operators affects neither stability nor accuracy of the Runge–Kutta discontinuous Galerkin method

Z. Sun]Zheng Sun
Department of Mathematics,
	The University of Alabama,
	Box 870350,
	Tuscaloosa, AL 35487, USA
zsun30@ua.edu

































	The Runge–Kutta (RK) discontinuous Galerkin (DG) method is a mainstream numerical algorithm  for solving hyperbolic equations. In this paper, we use the linear advection equation in one and two dimensions as a model problem to prove the following results: For an arbitrarily high-order RKDG scheme in Butcher form, as long as we use the P^k approximation in the final stage, even if we drop the kth-order polynomial modes and use the P^k-1 approximation for the DG operators at all inner RK stages, the resulting numerical method still maintains the same type of stability and convergence rate as those of the original RKDG method. Numerical examples are provided to validate the analysis. The numerical method analyzed in this paper is a special case of the Class A RKDG method with stage-dependent polynomial spaces proposed in <cit.>. Our analysis provides theoretical justifications for employing cost-effective and low-order spatial discretization at specific RK stages for developing more efficient DG schemes without affecting stability and accuracy of the original method.

[2020]65M12, 65M15


[
    [
    
=====





§ INTRODUCTION


In this paper, we perform stability analysis and optimal error estimates of a class of Runge–Kutta (RK) discontinuous Galerkin (DG) schemes for solving a linear hyperbolic equation. In these schemes, the DG operators at all inner RK stages adopt polynomial approximations of one degree lower than those of the standard RKDG schemes. The numerical method analyzed in this paper belongs to Class A of the RKDG method with stage-dependent polynomial spaces (sd-RKDG method), which we proposed in <cit.>. Hence, we will refer to it as the sdA-RKDG method. Its general formulation is given in (<ref>) and (<ref>).



 §.§ RKDG method and its variant.
 The DG method is a class of finite element methods that utilize discontinuous piecewise polynomial spaces. It was initially introduced by Reed and Hill in the 1970s for solving steady-state transport equations <cit.>. Later, in a series of works by Cockburn et al. <cit.>, the DG method was integrated with the Runge–Kutta (RK) time discretization method for solving time-dependent hyperbolic conservation laws. The resulting fully discrete method is referred to as the RKDG method. The RKDG method features many favorable properties, such as good stability, high-order accuracy, flexibility for handling complex geometries, and capability for sharp capturing of shocks. It has become one of the mainstream numerical methods for computational fluid dynamics <cit.>. 

The construction of the RKDG scheme follows the classical method-of-lines framework. For example, consider the initial value problem of the linear advection equation

    u_t + β·∇ u = 0,    u = u(x,t),   x∈Ω⊂ℝ^d,    d = 1,2,    0<t<T,

where β∈ℝ^d is a constant or a constant vector. To avoid unnecessary technicalities, we will consider a periodic rectangular domain with Cartesian meshes. If one wishes to construct a second-order RKDG scheme, the first step is to apply the DG spatial discretization to the linear advection equation (<ref>), which yields a semi-discrete scheme in the form of an ordinary differential equation (ODE) system

    / t u_h = L_h u_h,

where u_h is the semi-discrete solution, and L_h is the DG operator. To obtain the fully discrete scheme, a second-order RK method is used for time discretization. For example, using the explicit midpoint rule with a time step size of τ, we have

	
    u_h^(2) =    u_h^n +τ/2 L_h u_h^n,
    
    		u_h^n+1 =    u_h^n + τ L_h u_h^(2).


This formulation results in a stable second-order RKDG scheme.

However, in our recent work <cit.>, we observe that even if one uses the first-order P^0 approximation in the first stage (<ref>), the resulting numerical scheme is still stable and admits optimal accuracy. Let _h be the P^0 approximation of L_h. In other words, it is the P^1-DG operator without the evaluation of linear modes. See (<ref>) for the rigorous setup. Then the scheme

	
    u_h^(2) =    u_h^n +τ/2_h u_h^n,
    
    		u_h^n+1 =    u_h^n + τ L_h u_h^(2)


is still stable under the usual Courant–Friedrichs–Lewy (CFL) condition and maintains second-order accuracy. (<ref>) is an example of the sdA-RKDG scheme (<ref>).

In the implementation of the RKDG scheme (<ref>), one needs to evaluate the second-order P^1-DG operator L_h twice, while in the implementation of the sdA-RKDG scheme (<ref>), one only needs to evaluate a first-order P^0-DG operator once and a second-order P^1-DG operator once. The ratio of the total numbers of polynomial coefficients one needs to evaluate in (<ref>) and (<ref>) is d+2:2d+2, which is 75% in one dimension (1D, also for one-dimensional) and about 66.7% in two dimensions (2D, also for two-dimensional). As one can see, in (<ref>), by breaking the method-of-lines structure and allowing DG operators with reduced polynomial order, we can obtain a more efficient numerical scheme without sacrificing stability and optimal order of accuracy.



 §.§ An ODE argument for optimal accuracy.
 The stability and accuracy of (<ref>) were analyzed through the von Neumann approach in <cit.>. For an intuitive explanation of the optimal order of accuracy, although not rigorous, the following hand-waving argument may offer some insights[This argument, applied within the context of fully discrete truncation error analysis, helps explain the order of accuracy for finite difference schemes with perturbed stages. However, it may have issues when applied to DG schemes, as discussed in Section <ref>.]. Indeed, one can consider _h as a first-order perturbation of the second-order operator L_h, namely, _h = L_h + 𝒪(h), with h being the spatial mesh size. Substituting this expression into (<ref>), one obtains

    u_h^n+1 = (I+τ L_h + τ^2/2L_h^2)u_h^n + 𝒪(τ^2h).

Note that without the term 𝒪(τ^2h), (<ref>) gives back the scheme (<ref>). The perturbation in the first stage operator, after being multiplied by τ twice, contributes only to a third-order local truncation error 𝒪(τ^2h) under the usual time step constraint τ≤ Ch, which does not affect the global-in-time second-order convergence.

Indeed, the above analysis of the low-order (or low-precision) perturbation of RK stages has been systematically studied in the context of ODEs by Grant in <cit.>. See also <cit.>. As a special case of Grant's work, we have the following observation: For an RK method in Butcher form of rth order, if one perturbs all inner-stage operators by 𝒪(ε) and leaves the final stage untouched, the resulting perturbed RK method admits the local truncation error 𝒪(τ^r+1)+𝒪(ετ^2). Thanks to the extra τ and the fact τ∼ h, one can take the perturbation as large as ε∼ h^k at inner stages and still achieve optimal accuracy in the global error 𝒪(τ^r + h^k+1). In other words, for high-order RK schemes in Butcher form, reducing the accuracy of the inner-stage operators by one may not affect the overall convergence rate. In <cit.>, we used this idea to develop an RKDG method with compact stencils. In <cit.>, we further explored this idea and developed the sd-RKDG method, and Grant's theory may relate to the performance of the Class A schemes of the sd-RKDG method.



 §.§ Limitations of the existing analysis.
 Although the von Neumann analysis in <cit.> and the local truncation error analysis in <cit.> shed light on the stability and accuracy of the sdA-RKDG schemes, these analyses have essential limitations and should not be considered as a complete theoretical justification for the sdA-RKDG method. The von Neumann analysis only applies to linear problems with constant coefficients on periodic domains. Moreover, this approach requires to analyze matrix eigenvalues, which are difficult to obtain analytically for high-order schemes or in multiple dimensions. The local truncation error analysis of the perturbed RK method relies on the assumption that the ODE operator has bounded derivatives. However, in the context of the fully discrete schemes for partial differential equations (PDEs), the operators L_h and _h scale inversely proportional to the mesh size h, which violates the essential assumptions in numerical ODE analysis as the limit h→ 0 is approached. Indeed, in the design of sd-RKDG schemes <cit.>, we have encountered cases where the local truncation error analysis mismatches the actual convergence rate for the sd-RKDG schemes. It may indicate that the ODE approach in <cit.> may not fully explain the numerical performance of the PDE solvers. In fact, for the error estimates in this paper, a very different argument will be used for the proof.



 §.§ Literature review.
 The classical approach for analyzing finite element schemes is through an energy-type argument, in which one exploits the weak formulation of the numerical scheme and derives norm estimates by taking appropriate test functions <cit.>. In terms of the DG method for hyperbolic problems, although there are many earlier works on the analysis of the DG schemes for the steady-state transport equations <cit.> or semi-discrete DG schemes for the time-dependent hyperbolic conservation laws <cit.>, a systematic stability and error analysis was not readily available until recently. In <cit.>, Zhang and Shu analyzed the stability and error estimates of the second- and third-order RKDG schemes. These works also relate to <cit.> on the study of the stabilized finite element methods by Burman et al. and also to <cit.> on the stability analysis of the third-order RK method for generic linear semi-negative problems by Tadmor. The analysis of RKDG schemes of higher order was hindered by the stability analysis of the fourth-order RK scheme for some time, which was then resolved in <cit.>. With the experience on the fourth-order RK method <cit.>, the stability analysis for both the RK method for generic linear semi-negative problems and for the RKDG method was extended to very high order <cit.>. The stability theory also facilitated the error estimates of the RKDG schemes. See <cit.>.



 §.§ Challenges, main ideas, and main results.
 In this paper, we perform stability analysis and optimal error estimates of the sdA-RKDG method for the linear advection equation in both 1D and 2D. The analysis is conducted for schemes of arbitrarily high order.  

The major challenge in the analysis is that the sdA-RKDG method goes beyond the classical method-of-lines framework and involves two different spatial operators within each time step. While in most previous works, the fully discrete analysis typically only involves the same spatial operator across all stages. From this perspective, the challenge in the analysis of the sdA-RKDG method shares similarities with the analysis of the Lax–Wendroff DG method <cit.> and the RKDG method with stage-dependent numerical fluxes <cit.>. However, they are also very different in the sense that <cit.> mainly involve different flux terms, while the sdA-RKDG method involves operators associated with different polynomial spaces. This difference presents significant challenges: For the sdA-RKDG method, even the solution u_h^n itself may not be in the test function space that defines the low-order operator _h in  (<ref>).

The main idea to overcome these difficulties is the following. We will work on the strong form, instead of the weak bilinear form, of the sdA-RKDG scheme. We note that _h = (I-)L_h, where (I-) is the L^2 projection to the P^k-1 space (and  is the projection to the orthogonal complement of the P^k-1 space in the P^k space). Therefore, the _h operator can be formally embedded into the usual P^k-DG space, which avoids the complication of working with different spaces. For the stability analysis, we reformulate the sdA-RKDG scheme as the standard RKDG scheme plus some perturbation terms consisting of compositions of L_h and . The main results for the stability analysis are summarized below.

	
  * In each time step, the L^2 energy of the sdA-RKDG solution can be bounded by that of the standard RKDG solution along with some vanishing jump terms  (Lemma <ref>). Through this energy estimate, one can deduce that the sdA-RKDG method retains the same type of stability as the standard RKDG method (Theorem <ref>).

For the error estimates, we derive the error equation in a compact form. We design a class of special projection operators so that certain remainder terms in the error equation can be appropriately controlled or eliminated. These projection operators differ from those found in classical finite element literature, as their construction also incorporates the time step size τ. The main results for the error analysis are summarized below.


	
  * We construct a class of projection operators that are well-defined and admit the optimal approximation estimate under the usual CFL condition (Lemmas <ref> and <ref>). With these new projection operators, we prove the optimal error estimates of the sdA-RKDG schemes (Theorems <ref> and <ref>).


We want to remark that the error estimates in this paper are written in a different language from that in previous works <cit.>, where the analysis is conducted using the multi-stage weak form of the RKDG schemes. To leverage the stability estimates in Section <ref>, the error estimates in this paper are performed under the compact strong form of the sdA-RKDG method with all stages combined. The language in this paper is more similar to that in the classical work by Brenner et al.  <cit.>. However, despite the differences in notations at first glance, the two approaches are indeed closely interrelated. See Appendix <ref>. 



 §.§ Significance and novelty.
 In this work, we utilize the RKDG method for the linear advection equation as an illustrative example. Through systematic and rigorous analysis, we demonstrate the feasibility of employing cost-effective and coarse spatial discretization at specific RK stages for  developing a numerical scheme that is highly efficient while maintaining the stability and accuracy of the original method. Our findings provide valuable insights, laying the groundwork for future research of fully discrete RK solvers for time-dependent PDEs involving low-order perturbed spatial operators. The ideas and techniques presented in this paper may extend beyond DG spatial discretization, encompassing other methods such as the continuous finite element method and the spectral method.

Below, we summarize our major theoretical novelties and contributions.

	
  * This paper exemplifies a novel framework for conducting stability and error analysis for numerical PDE schemes beyond the method of lines.
	
  * In the stability analysis, we establish energy estimates bounding the L^2 norm of the high-order DG modes by jump terms. See Lemmas <ref>, <ref>, and <ref>. These estimates indicate that filtering out high-order DG modes may have effects comparable to imposing jump penalties. The modal filter and jump penalty are two very commonly used approaches for stabilizing DG schemes. Our estimates reveal the inherent connections between the two stabilization approaches and demonstrate the stability effects in the fully discrete setting.
	
  * In the error estimates, by adopting the compact strong form in the analysis, we gain further insights into the error accumulation through the RK iteration.
	
		
  ∘ Firstly, through the error equation (<ref>) and its norm estimates (<ref>), one can directly see how the norm increment or decrement of the RK scheme can affect the growth of the projected error ξ_h^n, which bridges the error estimates of the RKDG scheme with recent works on the norm estimates of the RK method for homogeneous linear semi-negative problems <cit.>.
		
  ∘ Secondly, by writing the error equation in the compact strong form, one can conveniently identify the remainder terms after a complete time step. It facilitates the design of novel projection operators (and typically, global projection operators) with information both in space and time for optimal error estimates.
	





 §.§ Organization of the paper.
 As we are performing a systematic and unified analysis for sdA-RKDG schemes of arbitrarily high order, the paper unavoidably becomes lengthy and technical. Below, we explain the organization and provide tips on reading the paper.

	
  * Since there are no significant differences in the 1D and 2D analyses, we use unified notations for both cases and conduct the analysis simultaneously. However, to avoid being distracted by the complications in the 2D notations, it would be helpful for readers to concentrate on the 1D case when verifying propositions and lemmas.
	
  * While the scope of the paper focuses on generic high-order sdA-RKDG schemes, to enhance readability, we begin each section on stability and error analysis with the second-order case (<ref>). The lemmas introduced in the analysis of the second-order scheme will serve as the base case for extension to higher-order schemes. To avoid being confused by the general high-order cases, we suggest the readers to skip Sections <ref>, <ref>, <ref>, and <ref> in the first-round reading.
	
  * The key to the stability analysis is to acquire estimates in Lemmas <ref> and <ref>, and the key to the optimal error estimates is to construct projection operators in Lemmas <ref> and <ref>. The proofs of these lemmas are technical. To avoid disrupting the flow of the paper, we postpone these proofs to Section <ref> near the end of the paper. 


The rest of the paper is organized as follows. In Section <ref>, we introduce the notations and preliminaries. Sections <ref> and <ref> provide stability analysis and optimal error estimates, respectively. In each of these sections, we begin with the second-order case before presenting results for higher-order cases. Section <ref> presents numerical examples to verify our theoretical analysis, while Section <ref> contains proofs of several technical lemmas. Finally, conclusions are given in Section <ref>.






§ NOTATIONS AND PRELIMINARIES




 §.§ DG spatial discretization






  §.§.§ 1D case
 

 

Consider the 1D linear advection equation

    u_t + β u_x = 0,    u = u(x,t),    x∈Ω⊂ℝ,    0<t<T.

We assume the periodic boundary condition, and β > 0 is a positive constant. 

Let Ω_h = {I_i}_i = 1^N, with I_i = (x_i-1/2, x_i+1/2) be a quasi-uniform mesh partition of Ω. We denote by h_i = x_i+1/2 - x_i-1/2 and h = max_i = 1,⋯, N h_i. The finite element space of the DG method is chosen as 

    V_h^ℓ = {v∈ L^2(Ω): v|_I_i∈P^ℓ(I_i)  ∀ i = 1,⋯,N},

where P^ℓ(I_i) is the linear space for polynomials on I_i of degree less than or equal to ℓ. 

Let wv = ∑_i=1^Nwv = ∑_i=1^N ∫_I_i w v  x be the L^2 inner product and v = √(vv) be the L^2 norm on Ω_h. We denote by L_h the DG approximation of the operator L = -β∂_x with the upwind flux, where L_h: V_h^k →V_h^k is defined such that

    L_h w_hv_h = βw_h(v_h)_x -β∑_i=1^N ((w_h^- v_h^-)_i+ - (w_h^- v_h^+)_i-)     ∀ v_h ∈V_h^k.

Here (v_h^±)_i+1/2 = lim_ε→ 0^± v_h(x_i+1/2+ε) is the limit of v_h from the left or the right of x_i+1/2. 


Moreover, we introduce the following notations associated with jumps and traces.

    [v_h] = v_h^+- v_h^-,   w_hv_h = β∑_i=1^N ([w_h][v_h])_i+1/2,   v_h= v_hv_h^1/2
 
and 

    v_h = (β∑_i=1^N ((v_h^-)_i+1/2^2+(v_h^+)_i+1/2^2))^1/2.





  §.§.§ 2D case


 

Consider the 2D linear advection equation

    u_t + β^x u_x + β^y u_y = 0,    u = u(x,y,t),    (x,y) ∈Ω⊂ℝ^2,    0<t<T.

To avoid unnecessary technicalities, we consider a rectangular domain with periodic boundary conditions.  β^x>0 and β^y>0 are positive constants. 

Let Ω_h = {K_ij} be a partition of the domain, where K_ij = (x_i-1/2,x_i+1/2)× (y_j-1/2,y_j+1/2), i = 1,⋯, N_x and j = 1,⋯, N_y. We denote by h_x^i = x_i+1/2 - x_i-1/2 and h_y^j = y_j+1/2 - y_j-1/2.  Then we set

    V_h^ℓ = {v∈ L^2(Ω): v|_K_ij∈P^ℓ(K_ij)  ∀ i,j}

where P^ℓ(K_ij) is the linear space for polynomials on K_ij of degree less than or equal to ℓ. 

Let wv = ∑_i=1^N_x∑_j=1^N_ywv_K_ij = ∑_i=1^N_x∑_j=1^N_y∫_K_ij w v  x  y be the L^2 inner product and v = √(vv) be the L^2 norm on Ω_h. We denote by L_h the DG approximation of the operator L = -β^x ∂_x-β^y ∂_y, where L_h: V_h^k → V_h^k is defined such that for all v_h ∈ V_h^k, 

    L_h w_hv_h  = w_hβ^x(v_h)_x+β^y(v_h)_y
       -∑_i=1^N_x∑_j=1^N_y∫_x_i-^x_i+β^y (w_h(x,y_j+^-) v_h(x,y_j+^-)-w_h(x,y_j-^-) v_h(x,y_j-^+)) x
       -∑_i=1^N_x∑_j=1^N_y∫_y_j-^y_j+β^x (w_h(x_i+^-,y) v_h(x_i+^-,y)-w_h(x_i-^-,y) v_h(x_i-^+,y)) y.

Here the superscript - (+) indicates the limit from the left or below (right or above). 

Moreover, we introduce the following notations associated with jumps and traces.

    [v_h]_i,j+ =    v_h(x,y_j+^+) -v_h(x,y_j+^-), 
    
    	[v_h]_i+,j =    v_h(x_i+^+,y) -v_h(x_i+^-,y),


    w_hv_h
    	=   ∑_i=1^N_x∑_j=1^N_y(∫_x_i-^x_i+β^y([w_h][v_h])_i,j+ x + ∫_y_j-^y_j+β^x 
    	([w_h][v_h])_i+,j y),
    v_h =   v_hv_h^1/2,

	
    v_h
    =   (∑_i=1^N_x∑_j=1^N_yβ^y∫_x_i-^x_i+(v_h(x,y_j+^-))^2 + (v_h(x,y_j+^+))^2 x . 
       .+ ∑_i=1^N_x∑_j=1^N_yβ^x∫_y_j-^y_j+(v_h(x_i+^-,y))^2 + (v_h(x_i+^+,y))^2
     y)^.


 



  §.§.§ Notations for both 1D and 2D cases.


 

Throughout the paper, we use C, possibly with subscripts, for a generic constant independent of the time step size τ and the spatial mesh size h. 

Consider a mesh partition of the domain Ω_h = {K}. The DG space is set as

    V_h^ℓ = {v∈ L^2(Ω): v|_K ∈ P^ℓ(K)  ∀ K∈Ω_h}.

The P^k-DG operator L_h: V_h^k → V_h^k is defined in (<ref>) and (<ref>). It can be shown that L_h is a semi-negative operator, as explained in Proposition <ref>. For more details, see <cit.>.

[Negative semi-definiteness of the DG operator] In both 1D and 2D, 

    w_hv_h = -(L_h + L_h^⊤)w_hv_h  ∀ w_h, v_h ∈ V_h^k.

and v_h^2 = -(L_h+L_h^⊤)v_hv_h≥ 0 for all v_h ∈ V_h^k. 


The semi-discrete P^k-DG scheme can be written as a semi-negative ODE system

    / t u_h = L_h u_h.



In the sdA-RKDG scheme, we also need the P^k-1-DG operator, which is defined as _h: V_h^k→ V_h^k-1 such that (please pay attention to the test function space)

    _h w_hv_h = L_h w_hv_h  ∀ v_h ∈V_h^k-1.

Let V_h^⊥ be the orthogonal complement of V_h^k-1 in V_h^k. We introduce the operators

    Π: L^2(Ω) → V_h^k,   : L^2(Ω) → V_h^k-1, : L^2(Ω) → V_h^⊥
 for the L^2 projection to V_h^k, V_h^k-1, and V_h^⊥, respectively. It can be seen that the P^k-1-DG operator _h:V_h^k→ V_h^k-1 satisfies 

    _h =  L_h = (I-)L_h

and can be formally regarded as an operator on V_h^k, namely, _h : V_h^k→ V_h^k. 

With the inverse estimate, it can be seen that

    v_h≤ Cv_h≤ Ch^-v_h  ∀ v∈ V_h^k.

Then using the definition of L_h, we have L_h v_h≤ Ch^-1v_h for all v_h ∈ V_h^k, or equivalently, L_h≤ Ch^-1. Let us denote by λ =  τ/h. Hence we have τ_h = τ L_h≤τL_h≤ Cλ. Without loss of generality, we assume λ≤ 1 throughout the paper.  





 §.§ RKDG and sdA-RKDG schemes


 

The standard RKDG scheme in Butcher form is written as


    u_h^(i) =    u_h^n + τ∑_j = 1^i-1 a_ijL_h u_h^(j),     i = 1,2,⋯, s,
    
    u_h^n+1 =    u_h^n + τ∑_i = 1^s b_iL_h u_h^(i).


We will only consider explicit RK schemes in this paper and assume a_ij = 0 if j≥ i. In this case, for an rth-order RK method, the RK operator for the linear autonomous problem (<ref>) should match its stability polynomial, and can hence be written in a compact form as

    u_h^n+1 = R_h u_h^n :=∑_i = 0^s α_i (τ L_h)^i u_h^n = (I + τ L_h ∑_i = 1^s α_i (τ L_h)^i-1) u_h^n,

where α_i = (i!)^-1 for  i = 0,1,⋯, r.

In the sdA-RKDG scheme, we replace the P^k-DG operator at inner stages by the P^k-1-DG operator _h. Then the scheme can be written as


    u_h^(i) =    u_h^n + τ∑_j = 1^i-1 a_ij_h u_h^(j),     i = 1,2,⋯, s,
    
    u_h^n+1 =    u_h^n + τ∑_i = 1^s b_iL_h u_h^(i).


It yields the compact form

    u_h^n+1 = _h u_h^n:=(I + τ L_h ∑_i = 1^s α_i (τ_h)^i-1) u_h^n,

where α_i = (i!)^-1 for  i = 0,1,⋯, r. In other words, we only use the P^k-DG operator L_h on the very left, and we use the P^k-1-DG operator _h at all other places. 


When there is no confusion, we will refer to (<ref>) as an RKrDGk scheme, and (<ref>) as an sdA-RKrDGk scheme. 

[Efficiency]
	Note that for linear problems, both (<ref>) and (<ref>) can be efficiently evaluated using Horner's method as:
	
    u_h^n+1 = (I + τ L_h(α_1 I + Z_h( α_2 I + Z_h(α_3 I +⋯ + Z_h (α_s-1I + α_s Z_h)))))u_h^n,

	where Z_h = τ L_h for the RKDG scheme and Z_h = τ_h for the sdA-RKDG scheme. Both methods can be evaluated with s matrix multiplications. Typically, for d≥ 2 and upwind DG schemes, (<ref>) requires even fewer floating-point multiplications than directly multiplying R_h or _h assembled beforehand.
	
	Since the ratio of the total degrees of freedom for P^k-1 and P^k polynomials is k:d+k, the ratio of computational operations involved in _h and L_h is approximately k:d+k. Hence, in each time step, the computational operations involved in the sdA-RKDG and the sd-RKDG methods are approximately (s-1)k + d + k : s(d+k) = sk + d : s(k+d). In the extreme case of solving high-dimensional problems with d≫ k, this ratio is approximately 1:s.
	
	In Section <ref>, we compare the efficiency of the sdA-RKDG and the RKDG schemes in 2D by analyzing the number of nonzero entries in the matrices involved in each time iteration. We present results in Table <ref>. The computational cost of the sdA-RKDG scheme ranges from approximately 70.6% to 86.1% of that of the RKDG method for 2≤ s = r = k+1≤ 9.


[RK method in Butcher form]
The analysis in this paper only considers the RK method in Butcher form and may not be directly applicable to the RK method in Shu–Osher form (also known as the strong-stability-preserving form <cit.>), unless a transformation is performed to rewrite a Shu–Osher form RK method back into the Butcher form. This is because our analysis relies on the following crucial fact: The operator on the very left of (<ref>) is the full DG operator L_h, and does not involve the reduced operator _h. Altering the inner stages of an RKDG scheme in Shu–Osher form could result in having _h present in the final stage, thereby disrupting the structure we have discussed. Such a method could yield a Class B sd-RKDG method, as described in <cit.>, which maintains optimal convergence for problems without sonic points but could suffer from accuracy degeneracy when sonic points occur.

We note that the RK method in Butcher form is also widely used for numerical solutions for hyperbolic conservation laws. See, for example, <cit.>.



§ STABILITY ANALYSIS




 §.§ Stability analysis: the second-order case




  §.§.§ Main result


 


In this section, we prove the stability of the sdA-RK2DG1 scheme (<ref>).
[Stability of sdA-RK2DG1]
The sdA-RK2DG1 scheme (<ref>) is monotonically stable, namely,  u_h^n+1≤u_h^n for λ≪ 1. 





  §.§.§ Preliminaries: stability of the RK2DG1 scheme


 

Before proceeding to prove Theorem <ref>, we first recall the stability analysis of the standard RK2DG1 scheme (<ref>), which can be written in the compact form as

    u_h^n+1 = (I + τ L_h + τ^2/2L_h^2)u_h^n:= R_h,2 u_h^n.
 

The energy estimate of the RK2 scheme can be found in many previous works. See, for example, <cit.>.
[Energy identity of RK2]
	
    R_h,2 w_h^2  = w_h^2 + τ^4/4L_h^2 w_h^2 - τw_h^2 - τ^2w_hL_h w_h - τ^3/2L_h w_h^2.


In general, the right hand side of (<ref>) may not be bounded by w_h^2 appropriately. Indeed, it is known that the RK2DGk method is unstable under the usual CFL condition if k≥ 2 <cit.>. But when coupled with P^1 elements, the RK2DG1 method is stable due to the following norm estimate, whose proof can be found in <cit.>. 
[Bounding high-order DG derivatives hy jumps]
	With P^1 elements, 
	
    τ^4L_h^2 w_h^2 ≤ C(τλ^3w_h^2 + τ^3λL_h w_h^2).


Simple algebra shows that τ^2w_hL_hw_h≤
3/4·τw_h^2 + 1/3·τ^3L_hw_h^2, and hence this cross term can be absorbed by the negative jump terms in (<ref>). Thus, Propositions <ref> and <ref> together imply the following stability results for the RK2DG1 scheme (<ref>). 
[Monotonicity stability of RK2DG1] Suppose λ≪ 1. Then there exists a positive constant ε>0 such that  
	
    R_h,2 w_h^2  ≤w_h^2 - ε(τw_h^2 + τ^3L_h w_h^2).






  §.§.§ Proof of Theorem <ref>


 

Now we prove Theorem <ref>. The main idea for the proof is to write the sdA-RK2DG1 scheme (<ref>) as a perturbation of the standard RK2DG1 scheme (<ref>). Recall that _h = (I-)L_h. Hence we have

    u_h^n+1 =    u_h^n + τ L_hu_h^n + τ^2/2L_h_h u_h^n = R_h,2 u_h^n - τ^2/2 L_h  L_h u_h^n.

Note that each term is in V_h^k. We can take the square on both sides of (<ref>) to get

    u_h^n+1^2 =   R_h,2 u_h^n^2 + τ^4/4L_h  L_h u_h^n^2 - τ^2L_h L_h u_h^nR_h,2u_h^n.

The estimate of R_h,2 w_h^2 has already been shown in Proposition <ref>. It remains to estimate L_h  L_h w_h^2 and terms of the form L_h L_h w_hv_h. 

To estimate such terms, we need to introduce Lemma <ref>, which basically states that the highest-order polynomial modes produced by the DG operator can be bounded by the jump term. Its proof is postponed to Section <ref>. 
[Bounding highest-order DG modes by jumps] 
	
		
    L_h w_hv_h≤ C h^-w_hv_h,
     L_h^⊤ w_hv_h≤ C  h^-w_hv_h.

	

With Lemma <ref>, we can estimate |L_h L_h w_hv_h| and L_h  L_h w_h in Lemmas <ref> and <ref>, respectively. 

	
    |L_h L_h w_hv_h| ≤ Ch^-1w_hv_h.



	Using the definition of the adjoint operator and the property of the L^2 projection, one can get
	
    L_h L_h w_hv_h =    L_h w_hL_h^⊤ v_h  = 	 L_h w_h L_h^⊤ v_h.

	From Lemma <ref>, we have L_h w_h≤ C h^-w_h and L_h^⊤ v_h≤ C h^-v_h. Therefore, applying the Cauchy–Schwarz inequality to the right side of (<ref>), we have 
	
    |L_h L_h w_hv_h|
    		≤ Ch^-w_h· Ch^-v_h≤  Ch^-1w_hv_h.



	
    L_h L_h w_h≤  Ch^-3/2w_h.



	Note that (<ref>) of Lemma <ref> implies
	L_h w_h≤ C h^-w_h. 	Recall that L_h v_h≤ Ch^-1v_h. Thus L_h L_h w_h≤ Ch^-1 L_h w_h≤ C h^-3/2w_h.


Now we continue to prove Theorem <ref>. Applying Lemma <ref>, we have

    τ^4/4L_h L_h u_h^n^2 ≤ C λ^3 τu_h^n^2.

Using Lemma <ref>, the triangle inequality, and Young's inequality ab≤ C(a^2+b^2), we have

    τ^2|L_h L_h u_h^nR_h,2u_h^n| ≤    Cλτu_h^nR_h,2u_h^n
    ≤    Cλτu_h^n( u_h^n + τL_h u_h^n + τ^2/2L_h^2 u_h^n)
    ≤    Cλτ(u_h^n^2 + τ^2 L_hu_h^n^2 + τ^4 L_h^2u_h^n^2).

Substituting (<ref>) and (<ref>) into (<ref>), it gives

    u_h^n+1^2 = _h,2 u_h^n^2≤R_h,2 u_h^n^2 +Cλ(τu_h^n^2 +τ^3 L_h u_h^n^2 + τ^5 L_h^2 u_h^n^2).

	Recall (<ref>). We have τv_h^2 ≤ Cλv_h^2. Thus, with Proposition <ref>, we get

    τ^5 L_h^2u_h^n^2 ≤ C λτ^4  L_h^2 u_h^n^2 ≤ C (λ^4 τu_h^n^2 + λ^2 τ^3 L_h u_h^n^2).

Substituting (<ref>) into (<ref>), one can get 
    u_h^n+1^2 ≤R_h,2 u_h^n^2 +Cλ(τu_h^n^2 +Cλτ^3 L_h u_h^n^2).
 
Applying the stability estimate of R_h,2u_h^n^2 in Proposition <ref>, it yields

    u_h^n+1^2 ≤u_h^n^2 +(Cλ-ε)τu_h^n^2 +(Cλ-ε)τ^3 L_h u_h^n^2.

Recall that ε is a fixed positive constant. Therefore, when λ≪ 1, we have Cλ-ε<0, which implies u_h^n+1^2 ≤u_h^n^2. This completes the proof of Theorem <ref>. 




 §.§ Stability analysis: general cases




  §.§.§ Main results


 

The main result of our energy estimate is given in Lemma <ref>, which states that the L^2 energy of the sdA-RKDG solution can be bounded by that of the standard RKDG solution along with some diminishing jump terms as λ approaches zero. The proof of Lemma <ref> is given in Section <ref>.
[Generalization of (<ref>)]
	_h w_h^2 ≤R_h w_h^2 + Cλ∑_i = 0^2s-2τ^2i+1L_h^i w_h^2. 


With Lemma <ref>, we can prove that the sdA-RKDG schemes share the same stability property as the standard RKDG schemes, as explained in Theorem <ref>. A sketch proof of the theorem is given in Section <ref>.
[Stability of sdA-RKDG schemes]
	Stability results in Theorem <ref> and Theorem <ref> remain valid after replacing R_h= ∑_i = 0^s α_i (τ L_h)^i by _h = I + τ L_h ∑_i = 1^s α_i (τ_h)^i-1.




  §.§.§ Preliminaries: stability of generic RKDG schemes


 

The RK schemes admit the following energy expansion <cit.>. 
[Energy identity of RK]
	Let R_h w_h  = ∑_i = 0^sα_i (τ L_h)^i w_h. Then
	
    R_h w_h^2 = ∑_i=0^s β_i τ^2iL_h^i w_h^2 +  ∑_i,j = 0^s-1γ_ijτ^i+j+1L_h^i w_hL_h^jw_h.

	Here 
	
    β_i = ∑_ℓ = max{0,2i-s}^min{2i,s}α_ℓα_2i-ℓ(-1)^i-ℓγ_ij = ∑_ℓ=max{0,i+j+1-s}^min{i,j}(-1)^min{i,j}+1-ℓα_ℓα_i+j+1-ℓ.



In general, the stability of the RK scheme is determined by the sign of the first nonzero β_i and the size of the largest negative-definite leading principal minor of (γ_ij). The general stability criteria with arbitrary coefficients α_i can be found in <cit.>. Below, we list the stability of r-stage r-th order RK scheme R_h = ∑_i = 0^r(i!)^-1 (τ L_h)^i. 

[Different types of stability] 
	
		
  * The method is weakly(γ) stable if 
		R_h w_h≤ (1+ Cλ^γ)w_h. This implies u_h^n≤ e^Ct^nu_h^0 if λ≤λ_0 h^1/(γ-1).  
		
  * The method is strongly(m_⋆) stable if R_h^m w_h≤w_h for m≥ m_⋆. This implies u_h^n+m≤u_h^n if λ≤λ_0 for some constant λ_0.  
		
  * The method is monotonically stable if it is strongly(1) stable. 
	

[Stability of the r-stage RKrDGk]
	Suppose L_h+L_h^⊤≤ 0 is negative semi-definite. Let R_h = ∑_i = 0^r(i!)^-1 (τ L_h)^i. For all k≥ 0, we have the following stability results.
	
		
  * If r ≡ 1 4, then the method is weakly(r+1) stable.  
		
  * If r ≡ 2 4, then the method is weakly(r+2) stable.
		
  * If r ≡ 3 4, then the method is monotonically stable.
		
  * If r ≡ 0 4, then the method is strongly(2) stable.
	

Additionally, similar to Proposition <ref> for the RK2DG1 case, by noting that L_h^i w_h can be bounded by jump terms L_h^jw_h when the polynomial degree k is small compared to i, the following improved stability results can be acquired <cit.>. 
[Stability of RKrDGk with low-order polynomials]
	Consider the RKrDGk method R_h=∑_i = 1^r(i!)^-1 (τ L_h)^i. We have the following stability results.
	
		
  * The method is monotonically stable if k≤⌊ r/2⌋ - 1. 
		
  * The method is strongly(m_⋆) stable for some m_⋆≥ 1 if k≤⌊ r/2⌋. If m_⋆ = 1, then the method is monotonically stable.
	 


[“Strong stability"]
	In <cit.>, we referred to strong(1) stability or monotonicity stability as “strong stability," while in <cit.>, the authors referred to the stability implied by strong(m_⋆) stability as the “strong stability." 	Here we avoid using the terminology “strong stability" to prevent ambiguity, as it has been referred to differently in previous works. 






  §.§.§ Proof of Lemma <ref>


 

To explain the motivation, we start by taking a glance at the third-order sdA-RKDG scheme u_h^n+1 = _h,3 u_h^n, where

    _h,3w_h = w_h + τ L_h w_h + τ^2/2L_h_h w_h + τ^3/6L_h_h^2w_h.

With _h = (I-)L_h, (<ref>) can be reformulated as 

    _h,3w_h =   (I + τ L_h + τ^2/2L_h^2 + τ^3/6L_h^3)w_h -τ^2/2 L_h  L_h w_h - τ^3/6L_h(L_h^2-_h^2)w_h
    
    		=   R_h,3 w_h -τ^2/2 L_h  L_h w_h - τ^3/6(L_h^2 L_h +L_h L_h^2 -L_h( L_h)^2)w_h.

From (<ref>), it can be seen that in general, we need to deal with operators in the form of a noncommunicative product of L_h^i and . To this end, we need to introduce some notations. 
We denote by = (i_1,i_2,⋯, i_[]) a vector with positive integer entries. Let []≥2 be the length of the vector  and || = i_1 + i_2 + ⋯ + i_[] be the ℓ^1 norm of the vector . Define 

    L_h^ = L_h^i_1 L_h^i_2⋯ L_h^i_[].


Reformulation of the sdA-RKDG scheme.  Subtracting (<ref>) from (<ref>) and replacing u_h^n by w_h, we have

    _h w_h = R_h w_h + ∑_j=1^s α_j τ^j L_h(_h^j-1-L_h^j-1) w_h.

We can expand L_h(_h^i-1 - L_h^i-1) to rewrite the sdA-RKDG scheme in the following form. 

    _h w_h = R_h w_h + ∑_: ||≤ sμ_τ^|| L_h^ w_h.

Here the real numbers {μ_} are known and depend only on {α_i}. As their specific values are not essential, we omit the characterizations of the coefficients μ_. 


Energy estimates. By squaring both sides of (<ref>) and using the elementary inequality (∑_i|a_i|)^2≤ C∑_i|a_i|^2, one can see that the sdA-RKDG scheme admits the following energy estimate. 

    _hw_h^2 =   R_h w_h^2 +  ∑_: ||≤ sμ_τ^||L_h^ w_h^2 + 2∑_: ||≤ sμ_τ^||L_h^ w_hR_h w_h
    ≤   R_h w_h^2 +  C∑_: ||≤ sτ^2||L_h^ w_h^2 + C∑_: ||≤ sτ^|||L_h^ w_hR_h w_h|.


Similar to the second-order case, we need to estimate terms of the forms L_h^ w_h^2 and L_h^ w_hv_h. To this end, we need the following Lemmas <ref> and <ref>.  The proof of Lemma <ref> is technical and is  therefore postponed to Section <ref>. 
[Generalization of Lemma <ref>]
	
    |L_h^ w_hv_h|
    			≤   C∑_j=0^i_1-1h^-||+i_[]+jL_h^i_[]-1 w_hL_h^j v_h .
	


[Generalization of Lemma <ref>]
	
    L_h^ w_h≤   C∑_j=0^i_1-1h^-||+i_[]-L_h^i_[]-1 w_h.
	


We use the inverse estimate (<ref>) and the fact L_h≤ Ch^-1 to get 

    L_h^j v_h≤ Ch^-L_h^j v_h≤ Ch^-L_h^j v_h≤ Ch^-j-v_h.
 
Hence Lemma <ref> implies

    |L_h^ w_hv_h|
    ≤   C∑_j=0^i_1-1h^-||+i_[]-L_h^i_[]-1 w_hv_h .
	
The proof of Lemma <ref> can be completed after taking v_h = L_h^ w_h in (<ref>) and dividing by L_h^ w_h on both sides. 


Now we continue to prove Lemma <ref>. To estimate the second term in (<ref>), we apply Lemma <ref> to get

    τ^2||L_h^ w_h^2
    ≤   C∑_j=0^i_1-1λ^2||-2i_[]+1τ^2i_[]-1L_h^i_[]-1 w_h^2≤ Cλτ∑_j=0^i_1-1(τ L_h)^i_[]-1 w_h^2.
	
To estimate the last term in (<ref>), we set v_h = R_h w_h in Lemma <ref>. Then

    |τ^|| L_h^ w_hR_h w_h|
    		≤   Cτ∑_j=0^i_1-1λ^||-i_[]-j(τ L_h)^i_[]-1 w_h(τ L_h)^j R_hw_h.
	
With the triangle inequality, we have (τ L_h)^j R_hw_h≤∑_ℓ=0^s|α_ℓ| (τ L_h)^ℓ+j w_h.	Therefore, 

    |τ^|| L_h^ w_hR_h w_h|
    		≤   Cτ∑_j=0^i_1-1∑_ℓ = 0^sλ^||-i_[]-j(τ L_h)^i_[]-1 w_h(τ L_h)^ℓ + j w_h.
	
Note that ||-i_[]-j ≥ i_1 - j ≥ 1. Hence λ^||-i_[]-j≤λ. Applying Young's inequality ab≤ C(a^2+b^2) gives

    |τ^|| L_h^ w_hR_h w_h| ≤ Cλτ∑_j=0^s+i_1-1(τ L_h)^j w_h^2.

Substitute (<ref>) and (<ref>) into 
(<ref>). Recalling that 1≤ i_1, i_[]≤ s-1, one can get

    u_h^n+1^2 = _h w_h^2 ≤R_h w_h^2 + Cλτ∑_j = 0^2s-2(τ L_h)^j w_h^2.
 
This completes the proof of Lemma <ref>. 

[Stability of other Class A sd-RKDG ] In addition to (<ref>), we have also considered other sd-RKDG schemes of Class A in <cit.>. In such schemes, each DG operator in (<ref>) does not have to be the P^k-1-DG operator _h and can instead be chosen as either a P^k-DG operator L_h or a P^k-1-DG operator _h. It is important to note that such schemes can also be expressed in the generic form of (<ref>), and all previous analyses can proceed accordingly. Hence, Lemma <ref> holds for all Class A sd-RKDG schemes in <cit.>. Consequently, with the analysis in Section <ref>, one can show that all Class A sd-RKDG schemes will inherit the same stability property as their parent RKDG schemes.





  §.§.§ Sketch proof of Theorem <ref>


 

	From Lemma <ref>, one can see that the stability of the RKDG and sdA-RKDG schemes is determined by the sign of the leading term of β_i and the size of the largest negative-definite leading principal minor of (γ_ij). Note that when λ is sufficiently small, the additional term Cλ∑_i=0^2s-2τ^2i+1L_h^i w_h^2 in Lemma <ref> neither affects β_i nor changes the negative definiteness of the leading principal minors in (γ_ij). Hence, all the stability results of the standard RKDG methods will be retained. 
	
	Below, as an example, we give the proof of the third results in Theorem <ref>. The proof of other cases is omitted. 
	
	Let ζ denote the index such that β_i = 0 for i < ζ and β_ζ≠ 0. According to <cit.>, for r ≡ 3 4, it follows that ζ = (r + 1)/2, β_ζ = 2(-1)^ζ+1/(r+1)! < 0, and Γ_ζ = (r_ij)_0≤ i,j≤ζ < 0 is negative definite. By following the proof outlined in <cit.>, one can derive <cit.>. After simplifying the coefficients there, it can be seen that there exists a positive constant ε>0 such that
	
    R_h w_h^2 ≤w_h^2 + (β_ζ + Cλ) τ^2ζL_h^ζ w_h^2 - ε∑_i = 0^ζ-1τ^2i+1L_h^i w_h^2.
 
	Note that for i≥ζ, we can use the inverse estimate v_h≤ Ch^-1/2v_h  to get 
	
    τ^2i+1L_h^i w_h^2 ≤ Cλτ^2iL_h^i w_h^2 ≤ Cλ^2i-2ζ+1τ^2ζL_h^ζ w_h^2 ≤ Cλτ^2ζL_h^ζ w_h^2.

	Therefore, 
	
    Cλ∑_i = 0^2s-2τ^2i+1L_h^i w_h^2 ≤ Cλ∑_i = 0^ζ-1τ^2i+1L_h^i w_h^2 + Cλτ^2ζL_h^ζ w_h^2.

	Then substituting estimates (<ref>) and (<ref>) into Lemma <ref>, we have 
	
    _h w_h^2 ≤w_h^2 + (β_ζ + Cλ) τ^2ζL_h^ζ w_h^2 - (ε-Cλ) ∑_i = 0^ζ-1τ^2i+1L_h^i w_h^2.

	Since β_ζ < 0 and ε>0, we have _h w_h^2 ≤w_h^2 when λ≪ 1, which completes the proof of monotonicity stability of _h for r ≡ 3 4. 
	


§ OPTIMAL ERROR ESTIMATES




 §.§ Preliminaries: projection operators
 



  §.§.§ 1D case
 

 

We denote by H^ℓ(Ω_h) = {v∈ L^2(Ω): v|_I_i∈ H^ℓ(I_i) ∀ i = 1,⋯,N} and v_H^ℓ = (∑_i=1^N v_H^ℓ(I_i)^2)^1/2.
The following Gauss–Radau projection <cit.> is used for the optimal error estimate of the DG method in 1D. : H^1(Ω_h) → V_h^k is an operator such that for w∈ H^1(Ω_h), the projection error η =  w - w satisfies


    ηv_h =     0   ∀ v_h ∈ P^k-1(I_i)   ∀ i = 1,⋯, N,
    η_i+^- =     0  ∀ i  = 1,⋯, N.



The Gauss–Radau projection  satisfies the following approximation estimate and superconvergence property. 
[Gauss–Radau projection]
	: H^1(Ω_h)→ V_h^k in (<ref>) is well-defined. Furthermore, when ♭≥ 1, 
	
		
  * for w∈  H^♭(Ω_h), we have the approximation estimate
		
    (I- ) w≤ C w_H^min(ℓ,♭)h^min(ℓ,♭)    ∀ 0≤ℓ≤ k+1;

		
  * for w∈ C(Ω)∩ H^1(Ω_h),
		we have the superconvergence property
		
    ( L - L_h)wv_h = 0.

	 
	





  §.§.§ 2D case
 

 

Let H^ℓ(Ω_h) = {v∈ L^2(Ω): v|_K_ij∈ H^ℓ(K_ij) ∀ i,j} and  v_H^ℓ = (∑_i,jv_H^ℓ(K_ij)^2)^1/2. In addition, we assume uniform meshes for which h_x^i ≡ h_x and h_y^j ≡ h_y are constants. The following projection <cit.> is used for the optimal error estimate of the P^k-DG method in 2D. : H^1(Ω_h) → V_h^k is an operator such that for w∈ H^1(Ω_h), the projection error η =  w - w satisfies


    η1_K_ij =  0   ∀  i,j,
       ηβ^x(v_h)_x+β^y(v_h)_y_K_ij
    	-∫_x_i-^x_i+β^y η(x,y_j+^-)( v_h(x,y_j+^-)- v_h(x,y_j-^+)) x
       -∫_y_j-^y_j+β^x η(x_i+^-,y) ( v_h(x_i+^-,y)-v_h(x_i-^+,y)) y = 0   ∀ v_h ∈ P^k(K_ij) ∀ i,j.



The special projection  satisfies the following approximation and superconvergence properties, as explained in Lemma <ref>. The case ℓ = k+2 for (<ref>) is proved in <cit.>. We have supplemented a sketched proof in Section <ref> for the case 1≤ℓ≤ k+1. 
[Liu–Shu–Zhang projection]
	Let Ω_h be a uniform partition in 2D. : H^1(Ω_h) → V_h^k in (<ref>) is well-defined. Furthermore, when ♭≥ 1, 
	
		
  * for w∈ H^♭(Ω_h), we have the approximation estimate
		
    (I- ) w≤ C w_H^min(ℓ,♭)h^min(ℓ,♭)    ∀ 0≤ℓ≤ k+1;

		
  * for w∈ C(Ω)∩ H^♭(Ω_h), we have the superconvergence property
		
    ( L - L_h)wv_h≤ C w_H^min(ℓ,♭)h^min(ℓ,♭)-1v_h    ∀ 1≤ℓ≤ k+2.

	





  §.§.§ Notations for both 1D and 2D cases


 

In order to obtain the optimal error estimates for both the 1D and 2D cases in one go, in this section, we unify the notations and write properties in Proposition <ref> and Lemma <ref> as a single proposition. Note that the regularity assumptions for the superconvergence properties in 1D and 2D are different, and we have adopted the stronger assumption below for (<ref>). 

[Special projection ]
	On both 1D quasi-uniform meshes and 2D uniform meshes, there exists a projection operator : H^1(Ω_h)→ V_h^k, such that when ♭≥ 1, 
	
		
  * for w∈ H^♭(Ω_h), we have the approximation estimate
		
    (I- ) w≤ C w_H^min(ℓ,♭)h^min(ℓ,♭)    ∀ 0≤ℓ≤ k+1;

		
  * for w∈ C(Ω)∩ H^♭(Ω_h), we have the superconvergence property
		
    ( L - L_h)wv_h≤ C_d w_H^min(ℓ,♭)h^min(ℓ,♭)-1v_h    ∀ 1≤ℓ≤ k+2.

	


[Error estimates in other settings]
The optimal error estimates in later sections can be extended to other settings, as long as one can construct a projection operator with the properties outlined in Proposition <ref>. For example, using the Generalized Gauss–Radau projection in <cit.>, one can prove optimal error estimates for sdA-RKDG schemes with upwind-biased fluxes in 1D; with the projections from <cit.>, one can prove optimal error estimates for the sdA-RKDG schemes with upwind or upwind-biased fluxes on special simplex meshes satisfying the so-called flow conditions in two and three dimensions.


In our error estimates, we need to deal with functions involving both spatial and temporal variables. We define

    L^∞(0,T;H^♭(Ω_h)) = {u(·,t):[0,T]→ H^♭(Ω_h)|sup_0≤ t≤ Tu(·,t)_H^♭ <+∞}
 
and denote by u_L^∞(H^♭) = sup_0≤ t≤ Tu(·,t)_H^♭. Here the “sup" function should be interpreted as the essential supreme. We will abuse the notation and formally identify u(x,t): Ω× [0,T]→ℝ with u(·,t): [0,T]→ H^♭(Ω_h), and use u∈ C^♭-1(Ω× [0,T])∩ L^∞(0,T;H^♭(Ω_h)) to represent u(x,t)∈ C^♭-1(Ω× [0,T]) and u(·,t)∈ L^∞(0,T;H^♭(Ω_h)).



  §.§.§ A few useful estimates


 

The following lemma, derived from Proposition <ref>, will be frequently used. Its proof is postponed to Section <ref>.

	On both 1D quasi-uniform meshes and 2D uniform meshes, we have the following results. 
	
		
  * If w ∈ H^♭(Ω_h) and ♭≥ 1, then
		
		(I-) w ≤   C w_H^min(ℓ,♭)h^min(ℓ,♭),        = 
		
		Π, ,     ∀0≤ℓ≤k+1,

		,	      ∀0≤ℓ≤k.
		 
		
		
  * If w ∈ H^♭(Ω_h) and ♭≥ 2, then
			
		(I-) L w ≤    C w_H^min(ℓ,♭)h^min(ℓ,♭)-1,        = 
		
		Π, ,      ∀1≤ℓ≤k+2,

		,      ∀1≤ℓ≤k+1.
		
		
		
  * If w ∈ C^♭-1(Ω)∩ H^♭(Ω_h) and ♭≥ 1, then
		
    ( L - L_h )L^iw≤ C_d w_H^min(ℓ,♭)h^min(ℓ,♭)-i-1    ∀1≤ℓ≤ k+2,

		where 0≤ i ≤min(ℓ,♭)-1.
		
  * If w ∈ C^♭-1(Ω)∩ H^♭(Ω_h) and ♭≥ 1, then 
		
    ( L- _h ) w≤ Cw_H^min(ℓ,♭)h^min(ℓ,♭)-1    ∀1≤ℓ≤ k+1.

	





 §.§ Optimal error estimates: the second-order case




  §.§.§ Main results


 

In this section, we prove the optimal error estimate of the sdA-RK2DG1 scheme (<ref>), or equivalently, (<ref>). 
[Optimal error estimate of sdA-RK2DG1]
	Consider the sdA-RK2DG1 scheme (<ref>) on either 1D quasi-uniform meshes or 2D uniform meshes. Suppose the exact solution u to (<ref>) is sufficiently smooth, for example, with u ∈ C^2(Ω×[0,T])∩ L^∞(0,T;H^3(Ω_h)), then the sdA-RK2DG1 scheme (<ref>) admits the optimal error estimate
	
    u^n - u_h^n≤u^0 - u_h^0+ C(t^n+1)u_L^∞(H^3)(τ^2 + h^2).

	Here t^n = nτ≤ T, u^n = u(x,t^n), and u^0 = u(x,0). 


The key to the proof of Theorem <ref> is to construct the special projection in Lemma <ref>, whose proof is postponed to Section <ref>.

[Projection for sdA-RK2DG1]
	Suppose w∈ C^♭-1(Ω)∩ H^♭(Ω_h) and ♭≥ 2. Then
	
    w = (I+τ/2_h)^-1(w+τ/2 L w)

	is well-defined when λ≪ 1. Moreover, it admits the approximation estimate
	
    (I- ) w≤ C w_H^min(ℓ,♭)h^min(ℓ,♭)    ∀1≤ℓ≤ k+1.

	
	
	
	
	




  §.§.§ Proof of Theorem <ref>


 

Reference function. Using the Taylor expansion and the equation u_t = Lu, it yields

    u^n+1 = u^n + τ u_t^n + τ^2/2u_tt^n + τρ^n = u^n + τ L u^n + τ^2/2 L^2 u^n +τρ^n,

where ρ^n is the remainder term in the Taylor expansion defined as

    ρ^n(x) = 1/2τ∫_0^τ u_ttt(x,t^n+κ)(τ-κ)^2κ.

Here one can show that u_ttt = ∂_t L^2 u = L^2 u_t = L^3 u is well-defined. The exchange of L^2 and ∂_t in the second equality can be justified by using the definition of the weak derivatives, integrating by parts, and exchanging the partial derivatives at the smooth test function.  

Applying the projection Π to both sides of (<ref>), it gives

    u^n+1 =  u^n + τ L u^n + τ^2/2 L^2 u^n +τρ^n,

which, after adding and subtracting terms, can be rewritten as 

    u^n+1 =  u^n + τ L_h  u^n + τ^2/2 L_h_h  u^n +τ(ρ^n + ν_h^n+ ζ_h^n).

Here

    ν_h^n =   τ^-1(-)(u^n+1-u^n):=(-)φ^n,
    ζ_h^n =    L(I+τ/2L )u^n - L_h(I+τ/2_h) u^n.

We will use (<ref>) as a reference function in our error estimate.  

Estimate of ρ^n. To estimate (<ref>), note that u_ttt = L^3 u. One can get

    ρ^n≤ρ^n≤ Csup_0≤ t≤ t^nu_ttt(·,t)τ^2  = Csup_0≤ t≤ t^nL^3u(·,t)τ^2≤ Cu_L^∞(H^3)τ^2.
 

Estimate of ν_h^n. Note that 

    φ^n = τ^-1(u^n+1-u^n) = τ^-1∫_0^τ u_t(x,t^n+κ)κ = τ^-1∫_0^τ L u(x,t^n+κ)κ.

One can apply Jensen's inequality to see that φ^n_H^ℓ≤sup_0≤κ≤τL u(·,t^n+κ)_H^ℓ≤u_L^∞(H^ℓ+1). Therefore, u∈ L^∞(0,T;H^3(Ω_h)) implies φ^n ∈ H^2(Ω_h). Using the triangle inequality, (<ref>), and (<ref>), with ♭ = 2 and ℓ = 2, we have

    ν_h^n≤(I-)φ^n+(I-Π)φ^n≤  Cφ^n_H^2h^2≤ Cu_L^∞(H^3)h^2.


Estimate of ζ_h^n. According to the definition of  in Lemma <ref>, we have 

    ζ_h^n =    L(I+τ/2L )u^n - L_h(I+τ/2 L)u^n
    		= ( L - L_h )(I+τ/2 L )u^n.

Apply the triangle inequality and use (<ref>) in Lemma <ref> with i = 0,1, ℓ = k+2 = 3, and ♭ = 3. Then

    ζ_h^n≤   ( L - L_h ) (u^n + τ/2L u^n)
    ≤   ( L - L_h ) u^n + τ/2( L - L_h ) L u^n
    ≤    Cu_L^∞(H^3)(h^2+τ h) ≤  C u_L^∞(H^3)h^2.




Estimate of the numerical error. Subtract (<ref>) from (<ref>), and define ξ_h^n =  u^n - u_h^n. It gives

    ξ_h^n+1 =   ξ_h^n + τ L_h ξ_h^n + τ^2/2 L_h_h ξ_h^n  + τ (ρ^n +ν_h^n+ ζ_h^n )
    
    	=   _h,2ξ_h^n + τ (ρ^n +ν_h^n+ ζ_h^n ).

Applying the triangle inequality  and the monotonicity stability _h,2ξ_h^n≤ξ_h^n in Theorem <ref>, it gives

    ξ_h^n+1≤ξ_h^n +  τ(ρ^n +ν_h^n+ ζ_h^n).

With the estimates for ρ^n in (<ref>), for ν_h^n in (<ref>), and for ζ_h^n in (<ref>), we have

    ξ_h^n+1≤ξ_h^n + Cτu_L^∞(H^3)(τ^2 + h^2).

Note that τ n = t^n. Repeatedly applying (<ref>) for n times gives 

    ξ_h^n≤ξ_h^0 + Ct^nu_L^∞(H^3)(τ^2 + h^2).

Recalling the approximation estimate (<ref>) in Lemma <ref>, we have  u^n- u^n≤ Cu^n_H^2h^2 and u^0- u^0≤ Cu^0_H^2h^2. With the triangle inequality, it gives

    u^n - u_h^n≤   u^n- u^n+ξ_h^n≤ξ_h^0+ C(t^n+1)u_L^∞(H^3)(τ^2 + h^2)
    ≤   u^0 -  u^0  + u^0 - u_h^0+ C(t^n+1)u_L^∞(H^3)(τ^2 + h^2)
    ≤   u^0 - u_h^0+ C(t^n+1)u_L^∞(H^3)(τ^2 + h^2).


[Idea behind the definition of ]
	The relationship (<ref>) and the cancellation in (<ref>) motivate the definition of  in (<ref>). In order to get rid of (I+τ_h/2) in ζ_h^n in (<ref>), we introduce the term (I+τ_h/2)^-1 in the definition of . Then in order to pull out the operator L - L_h in (<ref>), we impose the term (I+τ L/2).



[Alternative error estimates for RK2DG1]
	When we set _h = L_h, the sdA-RK2DG1 scheme (<ref>) formally retrieves the RK2DG1 scheme. It is worth noting that in this case, all the analysis in this section still holds. Consequently, it provides an alternative proof of the optimal error estimates for the standard RK2DG1 scheme. Further discussions on this matter can be found in Appendix <ref>.






 §.§ Optimal error estimates: general cases




  §.§.§ Main results


 

Theorem <ref> presents the optimal error estimates of the sdA-RKrDGk scheme for generic values of r and k. Its proof is given in Section <ref>. 

[Optimal error estimates of sdA-RKrDGk]
	Consider the sdA-RKrDGk scheme (<ref>) (or (<ref>)) on either 1D quasi-uniform meshes or 2D uniform meshes. Suppose the exact solution u to (<ref>) is sufficiently regular, for example, u∈ C^♭-1([0,T]×Ω)∩ L^∞(0,T;H^♭(Ω_h)). Let q = min(r,k+1,♭-1). Then we have
	
    u^n - u_h^n≤ C_1(t^n)u^0 - u_h^0+ C_2(t^n)u_L^∞(H^q+1)(τ^q+ h^q),

	where t^n = nτ≤ T, u^n = u(x,t^n), u^0 = u(x,0), and
	
		
  * C_1(t^n) = C_2(t^n) = e^Ct^n if _h is weakly(γ) stable and λ≤λ_0 h^1/(γ-1);
		
  * C_1(t^n) = 1 and C_2(t^n) = C(t^n+1) if _h is monotonically stable and λ≪ 1;
		
  * C_1(t^n) = C and C_2(t^n) = C(t^n+1) if _h is strongly(m_⋆) stable and λ≪ 1.
	


The key step for proving Theorem <ref> is to construct the special projection in Lemma <ref>, whose proof is postponed to Section <ref>.

[Projection for sdA-RKrDGk]
	For w ∈ C^♭-1(Ω)∩ H^♭(Ω_h) and q ≤min(r,k+1,♭-1), the projection 

    w = (∑_i = 1^s α_i (τ_h)^i-1)^-1((∑_i = 1^q (i!)^-1 (τ L)^i-1)w + τ^q∑_i = q+1^s α_i (τ_h)^i-1-q L^q w).

is well-defined. Moreover, it satisfies the approximation estimate

    (I-) w≤ Cw_H^qh^q.



[Regularity assumption required by ]
	Following Lemma <ref> for the second-order case, one may naturally consider defining 
	
    w = (∑_i = 1^s α_i (τ_h)^i-1)^-1((∑_i = 1^s α_i (τ L)^i-1)w).

	However, this definition imposes a restrictive regularity assumption on w ∈ H^s(Ω_h), where s≥ r. For example, for r-stage rth-order RK schemes that are strongly(m_⋆) stable, as we need to combine m_⋆ steps into a single step in the analysis, we could have s = m_⋆ r. However, typically, we anticipate that H^r regularity should be adequate to achieve the desired accuracy. To circumvent the need for this additional regularity requirement, we modify high-order derivatives L^i-1w with q+1≤ i≤ s as _h^i-1-q L^q w. This modification shares a similar flavor to the cutting-off technique in the error analysis of the standard RKDG method <cit.>.




  §.§.§ Proof of Theorem <ref>


 

Reference function. Since u∈ C^♭-1([0,T]×Ω) and q≤♭-1, we have ∂_t^i u^n = L^i u^n for 1≤ i≤ q. Thus, using the Taylor expansion, we have

    u^n+1 = ∑_i = 0^q(i!)^-1τ^i∂_t^i u^n + τρ^n = ∑_i = 0^q(i!)^-1τ^iL^i u^n + τρ^n.

The remainder ρ^n is defined as 

    ρ^n(x) = 1/τ q!∫_0^τ∂_t^(q+1)u(x, t^n+κ)(τ-κ)^qκ.

Recalling that q≤♭ -1, one can show that ∂_t^(q+1)u = L^q+1 u is well-defined under the given regularity assumption. 


We apply Π to both sides of (<ref>) to get

    u^n+1 =  u^n + ∑_i = 1^q(i!)^-1τ^i L^i u^n + τρ^n.

By adding and subtracting terms, it can be rewritten as 

    u^n+1 =  u^n + ∑_i = 1^s α_i τ^i L_h _h^i-1 u^n + τ(ρ^n+ ν_h^n + ζ_h^n).

where 

    ν_h^n =   τ^-1(-)(u^n+1-u^n):=(-Π)φ^n,
    ζ_h^n =   ∑_i=1^q (i!)^-1τ^i-1 L^iu^n - L_h(∑_i=1^s α_i (τ_h)^i-1) u^n.

(<ref>) will be used as the reference function for our error estimate. 

Estimates of ρ^n and ν_h^n. The estimates of ρ^n and ν_h^n are similar to those in the second-order case. Further details are omitted.

    ρ^n≤ Csup_0≤ t≤ t^n∂_t^(q+1) u(·,t)τ^q = Csup_0≤ t≤ t^nL^q+1 u(·,t)τ^q≤ Cu_L^∞(H^q+1)τ^q.


    ν_h^n≤(I-)φ^n+(I-Π)φ^n≤  Cφ^n_H^qh^q≤ Cu_L^∞(H^q+1)h^q.

 

Estimate of ζ_h^n. We then estimate the term ζ_h^n. 


	When λ≪ 1, we have 
	
    ζ_h^n≤ C u_L^∞(H^q+1)(τ^q + h^q).



For ζ_h^n, we can use the definition of  to get

    ζ_h^n 
    	=   ( L - L_h )∑_i=1^q (i!)^-1 (τ L)^i-1u^n + τ^q L_h∑_i = q+1^s α_i (τ_h)^i-1-q L^q u^n
    
    	=   ζ_h^n,1 + ζ_h^n,2.


For ζ_h^n,1, we can apply a similar estimate as that for the second-order case. One can use the triangle inequality and (<ref>) in Lemma <ref> with ℓ = q+1≤♭ to get

    ζ_h^n,1≤   ∑_i=1^q (i!)^-1τ^i-1( L - L_h )L^i-1u^n
    ≤    C_d∑_i=1^q τ^i-1u^n_H^q+1h^q+1-i≤  C_du_L^∞(H^q+1)h^q.

For ζ_h^n,2, we separate the first term in the summation to get

    ζ_h^n,2 = τ^q (α_q+1 L_h L^q u^n + (τ L_h)∑_i = q+2^s α_i (τ_h)^i-2-q_h L^q u^n).

Then we apply the triangle inequality and the estimates τ_h≤τL_h≤ Cλ to get

    ζ_h^n,2≤    C τ^q (L_h L^q u^n + _h L^q u^n)
    ≤    C τ^q ((L_h -  L) L^q u^n +  (_h- L) L^q u^n+2 L^q+1 u^n).

We can apply (<ref>) in Lemma <ref> with w = u^n and ℓ = q+1≤♭ to get 

    (L_h -  L) L^q u^n≤ C_d u^n_H^q+1≤  C_du_L^∞(H^q+1).

Apply (<ref>) in Lemma <ref> with ℓ = 1 and w = L^qu^n ∈ C^♭ -q-1(Ω)∩ H^♭-q(Ω_h)⊂ C(Ω)∩ H^1(Ω_h). It gives

    (_h- L) L^q u^n≤ CL^q u^n_H^1≤ Cu_L^∞(H^q+1).

Moreover, 

    L^q+1 u^n≤ Cu^n_H^q+1≤ Cu_L^∞(H^q+1).
 
Substituting these estimates into (<ref>) gives

    ζ_h^n,2≤ C u_L^∞(H^q+1)τ^q.

Combining (<ref>) and (<ref>), we complete the proof of (<ref>). 


Estimate of the numerical error. Let ξ_h =  u - u_h. Subtracting (<ref>) from (<ref>) gives 

    ξ_h^n+1 = _h ξ_h^n + τ (ρ^n + ν_h^n+ ζ_h^n).

Apply the triangle inequality and the estimates in (<ref>), (<ref>), and (<ref>). Then we get

    ξ_h^n+1≤   _h ξ_h^n + τ(ρ^n+ν_h^n +ζ_h^n) 
    ≤   _h ξ_h^n + Cu_L^∞(H^q+1)τ(τ^q + h^q).


Below, we discuss the cases with different types of stability. 

	
  * Suppose _h is weakly(γ) stable. Then _h ξ_h^n≤ (1+Cλ^γ)ξ_h^n, and we have

    ξ_h^n+1≤  (1+ Cλ^γ)ξ_h^n + Cu_L^∞(H^q+1)τ(τ^q + h^q).

Assuming λ≤λ_0 h^1/(γ-1), we can apply discrete Gronwall's inequality to get

    ξ_h^n≤ e^Ct^n(ξ_h^0 + u_L^∞(H^q+1)(τ^q + h^q)).




  * Suppose _h is monotonically stable. Then _h ξ_h^n≤ξ_h^n, and we have

    ξ_h^n≤ξ_h^0 + Ct^nu_L^∞(H^q+1)(τ^q + h^q).




  * Suppose _h is strongly(m_⋆) stable, with m_⋆>1. Then _h^m_⋆w_h≤w_h. We need to replace _h by _h^m_⋆, n+1 by n+m_⋆, and τ by m_⋆τ in (<ref>), to get

    ξ_h^n+m_⋆≤ξ_h^n + Cm_⋆τu_L^∞(H^q+1)(τ^q + h^q).

From which we deduce that

    ξ_h^n_⋆≤ξ_h^0 + Ct^n_⋆u_L^∞(H^q+1)(τ^q + h^q),      n_⋆ =⌊ n/m_⋆⌋ m_⋆.

We need to discuss the last n - n_⋆ steps separately. Use the fact τ_h≤τL_h≤ C λ, we know that _h≤ C. Thus, (<ref>) implies 

    ξ_h^n+1≤ Cξ_h^n +  Cu_L^∞(H^q+1)τ(τ^q + h^q).

Therefore, with (<ref>) and (<ref>), we have

    ξ_h^n≤     Cξ_h^n_⋆ +Cu_L^∞(H^q+1)(n-n_⋆)τ(τ^q + h^q)
    ≤    Cξ_h^0 + C(t^n_⋆+(n-n_⋆)τ)u_L^∞(H^q+1)(τ^q + h^q)
    ≤    Cξ_h^0 + Ct^nu_L^∞(H^q+1)(τ^q + h^q).




The rest of the proof can be completed by applying the triangle inequality to (<ref>), (<ref>), and (<ref>), respectively. The derivation is similar to the second-order case, and details are omitted. Different constants can be attributed to the different types of stability.

[Alternative error estimates for RKrDGk]
	Similar to the second-order case, in the case _h = L_h, the analysis above gives an alternative proof for the optimal error estimates of the RKrDGk method.





§ NUMERICAL VERIFICATION


In this section, we verify our stability and error analysis presented in the previous sections through numerical experiments.

For the 1D tests, we consider the linear advection equation u_t + u_x = 0 on the interval x∈ [0,1], discretized with N mesh cells (h = 1/N). For the 2D tests, we consider the equation u_t + u_x + u_y = 0 on the domain (x,y)∈ [0,1]× [0,1], discretized with N× N mesh cells (h_x = h_y = 1/N). Both cases are subject to periodic boundary conditions.

In this study, we focus solely on linear advection equations. Additional numerical tests involving problems such as the nonlinear Burgers equation can be found in <cit.>.

We employ r-stage rth-order RK schemes with r = 2,3,4,5 <cit.>. For linear autonomous problems with a fixed r, all r-stage rth-order RK schemes are equivalent to the rth-order Taylor method, given by R_h = ∑_i = 0^r (i!)^-1τ^i L_h^i. Correspondingly, _h = I + τ L_h∑_i = 1^r(i!)^-1τ^i-1_h^i-1.



 §.§ Stability tests


 

	In this section, we verify different types of stability of the RKDG and sdA-RKDG schemes. To this end, we follow <cit.> to compute 
	
    δ = max((K_h)^m^2-1,10^-16)  with  
    	K_h = 
    	
    	R_h   for RKDG schemes,
    _h   for sdA-RKDG schemes.
 
	Here the operator norm can be evaluated by computing the 2-norm of the (sdA-)RKDG matrices assembled under the orthonormal Legendre basis. With τ = CFL/(dN) for the d-dimensional test (d =1,2), we note that δ is a function of the CFL number, N, and m.
	
	Below, we present plots showing the dependence of δ on the CFL number for various N values across different (sdA-)RKrDGk schemes. The 1D results, with N = 16, 32, 64, are shown in Figure <ref>. While the 2D results, with N = 4, 8, 16, are shown in Figure <ref>.
	
	Consistent with our theoretical analysis, the stability behavior appears to be independent of dimensionality, with 1D and 2D results exhibiting similar trends. Additionally, when r and k are fixed, both RKDG and sdA-RKDG schemes demonstrate similar patterns of δ regardless of the number of mesh cells. These results verify different types of stability outlined in Theorems <ref>, and <ref> (and also in Theorem <ref>). See Table <ref> for a summary.
	
	
	
	
	
	
	


	



 §.§ Accuracy tests with smooth solutions


 

	In this test, we consider the linear advection equation with sufficiently smooth solutions. For the 1D problem, we set u(x,0) = sin(2π x), and the exact solution is u(x,t) = sin(2π(x-t)). In the 2D problem, we set u(x,y,0) = sin(2π(x+y)), and the exact solution is u(x,y,t) = sin(2π(x+y-2t)). The final time is set as t^n = 1. We solve the problem using RKrDGk and sdA-RKrDGk methods with r = k+1 = 2,3,4,5. The time step size is set as τ = 0.1/(dN) for r = 2,3,4, and τ = 0.1/(dN^6/5) for r = 5.
	
	The L^2 errors of the numerical solutions on 1D uniform meshes, 1D randomly perturbed nonuniform meshes, and 2D uniform meshes are provided in Tables <ref>, <ref>, and <ref>, respectively. As observed, in all cases, the sdA-RKDG schemes achieve the optimal convergence rate, with numerical errors comparable to the corresponding RKDG schemes.








 §.§ Accuracy tests with solutions of limited regularity


 

	When r = k+1, Theorem <ref> indicates that the optimal rth-order convergence rate for RKrDGk schemes can only be attained if ♭≥ r+1 = k+2. To verify the sharpness of this regularity number, we consider the linear advection equation with the initial condition
	
    u(x,0) = (sin(2π x))^♭ -1/3

	for the 1D problem.  Note that the solution is in C^♭-1(Ω)∩ H^♭(Ω_h) but not in C^♭(Ω)∩ H^♭+1(Ω_h) <cit.>. When we set ♭ = r = k+1, observations from Tables <ref> and <ref> show that both the RKrDGk and the sdA-RKrDGk methods exhibit suboptimal convergence for r = k+1 = 2,3,4,5.  Here we have chosen a large final time to make the order degeneracy observable for r = k+1 = 4,5. Conversely, under the same settings with ♭ = r + 1 = k+2, all schemes achieve optimal convergence rates. A similar test has been conducted for the 2D case with the initial condition u(x,y,0) = (sin(2π (x+y)))^♭ -1/3, where r = k + 1 = 2,3. See Table <ref>.











 §.§ Tests of matrix sparsity


 

The efficiency of implementing RKDG and sdA-RKDG schemes can be assessed by counting the total number of nonzero entries in the matrices in the final matrix equation. Below, we consider the 2D problem with s = r = k+1. The DG spatial operator is implemented under the orthonormal Legendre basis. We consider four implementations: Horner's implementation of sdA-RKDG schemes, Horner's implementation of RKDG schemes, implementation of sdA-RKDG schemes by assembling _h, and implementation of RKDG schemes by assembling R_h. In Table <ref>, we summarize the number of nonzero matrix entries per mesh cell in a single time step (excluding identity matrices). The results indicate that for 2D problems, sdA-RKDG schemes may offer greater efficiency compared to standard RKDG schemes. For 2≤ s = r = k+1≤ 9, the computational cost of sdA-RKDG schemes ranges from 70.6% to 86.1% compared to that of RKDG schemes.









§ PROOFS OF TECHNICAL LEMMAS





 §.§ Proof of Lemma <ref>


 

We will consider the 1D and 2D cases separately. 

1D case. Using the self-adjoint property of , the definition of L_h, and the integration by parts, one can get

    L_h w_hv_h
    		=   L_h w_h v_h
    
    		=   βw_h( v_h)_x -β∑_i=1^N ((w_h^-( v_h)^-)_i+ - (w_h^-( v_h)^+)_i-)
    
    		=    -β(w_h)_x v_h - β∑_i=1^N  ([w_h]( v_h)^+)_i-.

Note that (w_h)_x ∈ V_h^k-1 and v_h ⊥ V_h^k-1. Thus, (w_h)_x v_h = 0 and 

    L_h w_hv_h = - β∑_i=1^N  ([w_h]( v_h)^+)_i-≤ Cw_h v_h.

Here we have used the Cauchy–Schwarz inequality. With the inverse estimate (<ref>), we get 

    v_h≤ Ch^-1/2 v_h≤ Ch^-1/2v_h.
 
Substituting (<ref>) into (<ref>) completes the proof of (<ref>) in the 1D case. 

For the proof of (<ref>), note we have <cit.>

    L_h^⊤ w_hv_h = -βw_h(v_h)_x +β∑_i=1^N ((w_h^+ v_h^-)_i+ - (w_h^+ v_h^+)_i-).
 
Comparing (<ref>) to (<ref>), L_h^⊤ can be defined by replacing the numerical flux w_h^- in (<ref>) by w_h^+ and changing the sign of each term. One can take advantage of this similar structure to prove (<ref>) by following similar lines as those in the proof of (<ref>). 

    L_h^⊤ w_hv_h
    = L_h^⊤ w_h v_h
    
    =    -βw_h( v_h)_x +β∑_i=1^N ((w_h^+( v_h)^-)_i+ - (w_h^+( v_h)^+)_i-)
    
    =   β(w_h)_x v_h + β∑_i=1^N  ([w_h]( v_h)^-)_i+
    = β∑_i=1^N  ([w_h]( v_h)^-)_i+
    ≤    C w_h v_h≤ Ch^-w_h v_h≤ Ch^-w_hv_h.



2D case. Similar to the 1D case, we use the self-adjoint property of , the definition of L_h, and the integration by parts. It yields


    L_h w_hv_h
    		= L_h w_h v_h
    
    		=   w_hβ^x( v_h)_x+β^y( v_h)_y
       -∑_i=1^N_x∑_j=1^N_y∫_x_i-^x_i+β^y (w_h(x,y_j+^-) ( v_h)(x,y_j+^-)-w_h(x,y_j-^-) (  v_h)(x,y_j-^+)) x
       -∑_i=1^N_x∑_j=1^N_y∫_y_j-^y_j+β^x (w_h(x_i+^-,y) ( v_h)(x_i+^-,y)-w_h(x_i-^-,y) ( v_h)(x_i-^+,y)) y
    
    		=    - β^x(w_h)_x + β^y(w_h)_y v_h - ∑_i=1^N_x∑_j=1^N_y∫_x_i-^x_i+β^y [w_h]_i,j-(  v_h)(x,y_j-^+) x
       -∑_i=1^N_x∑_j=1^N_y∫_y_j-^y_j+β^x [w_h]_i-,j ( v_h)(x_i-^+,y) y.

Note that the first term is 0 since (β^x∂_x + β^y∂_y)w_h∈ V_h^k-1 and v_h ⊥ V_h^k-1. As we have done for the 1D case, we can then prove the lemma by applying the Cauchy–Schwarz inequality, the inverse estimate, and the fact v_h≤v_h.  The 2D case of (<ref>) can be proved similarly by rewriting L_h^⊤ in a similar form as (<ref>). Details are omitted. 




 §.§ Proof of Lemma <ref>


 

To prove Lemma <ref>, we need a few preparatory lemmas. 



[Discrete integration by parts]
	
    L_h^i w_hv_h = ∑_j=0^i-1 (-1)^j+1L_h^i-j-1 w_hL_h^j v_h +(-1)^iw_hL_h^iv_h.
	


	Recall Proposition <ref>. The lemma can be proved by inductively applying L_h w_hv_h = -w_hv_h -w_hL_h v_h. Details are omitted. 	




	
    L_h^i L_h w_h≤ C h^-i-1w_h    ∀ i≥ 0.



	With the inverse estimate v_h≤ Ch^-v_h and the property L_h≤ C h^-1, we have
	
    L_h^i L_h w_h≤ Ch^-L_h^i L_h w_h≤ Ch^-L_h^i L_h w_h≤ Ch^-i- L_hw_h.

	Note from (<ref>), we have L_hw_h≤ Ch^-1/2w_h. Applying this estimate to (<ref>) completes the proof of (<ref>).  



	
	







Now we prove Lemma <ref>.


	We define z_h such that L_h^w_h = L_h^i_1 L_h z_h. In other words, 
	
    z_h =    L_h^i_2-1 L_h^i_3⋯ L_h^i_||w_h 
    
    			=   (L_h^i_2-1 L_h)
    			(L_h^i_3-1 L_h)⋯(L_h^i_||-1-1 L_h ) L_h^i_||-1w_h.

	First, with Lemma <ref>, the triangle inequality, and the Cauchy–Schwarz inequality, we have
	
    |L_h^ w_hv_h|
    
    			 =   |L_h^i_1-1(L_h L_h z_h)v_h|
    
    			=   |∑_j=0^i_1-2 (-1)^j+1L_h^i_1-j-1 L_h z_hL_h^j v_h +(-1)^i_1-1L_h L_h z_hL_h^i_1-1v_h|
    ≤   ∑_j=0^i_1-2L_h^i_1-j-1 L_h z_hL_h^j v_h+|L_h L_h z_hL_h^i_1-1v_h|.

	Then we apply Lemmas <ref> and <ref> to get
	
    |L_h^ w_hv_h|
    			≤   C∑_j=0^i_1-2h^-(i_1-j)z_hL_h^j v_h +Ch^-1 z_hL_h^i_1-1v_h
    
    			=   C∑_j=0^i_1-1h^-(i_1-j)z_hL_h^j v_h .
	
	Next, recall the definition of z_h in (<ref>) and apply Lemma <ref> for []-2 times. We have
	
    z_h≤ C h^-i_2-i_3-⋯ - i_[]-1L_h^i_[]-1w_h≤ Ch^-||+i_1+i_[]L_h^i_[]-1w_h.

	Finally, one can complete the proof of Lemma <ref> by substituting (<ref>) into (<ref>). 




 §.§ Proof of Lemma <ref>



 

	First, the well-definedness of  was prove in <cit.>. Then, (<ref>) follows from the standard approximation theorem <cit.>. See also <cit.>. 
	Next, for (<ref>), the case ℓ = k+2 is a Cauchy–Schwarz version of <cit.>. Finally, for the case 1≤ℓ≤ k+1, one can follow a similar argument as the standard suboptimal error estimate of the DG scheme.  To be more specific, with η =  w - w, we have
	
    ( L - L_h)wv_h = Lwv_h -  L_h wv_h = ηβ^x(v_h)_x+β^y(v_h)_y
       -∑_i=1^N_x∑_j=1^N_y∫_x_i-^x_i+β^y (η(x,y_j+^-) v_h(x,y_j+^-)-η(x,y_j-^-) v_h(x,y_j-^+)) x
       -∑_i=1^N_x∑_j=1^N_y∫_y_j-^y_j+β^x (η(x_i+^-,y) v_h(x_i+^-,y)-η(x_i-^-,y) v_h(x_i-^+,y)) y
       ≤ Cηv_h_H^1 + Cηv_h≤ C(h^-1η+h^-η)v_h.

	Note with the approximation property of η and the trace inequality, it can be seen that η≤ Cw_H^min(ℓ,♭)h^min(ℓ,♭) and η≤ Cη^1/2η_H^1^1/2≤ Cw_H^min(ℓ,♭)h^min(ℓ,♭)-1/2. Substituting these estimates into (<ref>) completes the proof of (<ref>) for the case 1≤ℓ≤ k+1. 





 §.§ Proof of Lemma <ref>


 

1. (<ref>) follows from (<ref>) and the standard approximation theory <cit.>. 
	
2. To prove (<ref>) for = Π,, we consider 
	(<ref>) with w = L w̃, ℓ = ℓ̃-1, and ♭ = ♭̃ -1. Note that, in the case 
	w̃∈ H^♭̃(Ω_h), 1≤ℓ̃≤ k+2 and ♭̃≥ 2, we have w∈ H^♭(Ω_h), 0≤ℓ≤ k+1 and ♭≥ 1. Hence (<ref>) can be applied, which gives
	
    (I-) L w̃≤ C L w̃_H^min(ℓ̃-1,♭̃ -1)h^min(ℓ̃-1,♭̃-1)≤ C w̃_H^min(ℓ̃,♭̃)h^min(ℓ̃,♭̃)-1.

	We can complete the proof of (<ref>) after removing the tilde symbol. The proof for = is similar. 
	
3. Similarly, for (<ref>), we consider (<ref>) with w = L^i w̃, ℓ =  ℓ̃-i, ♭ = ♭̃ -i, and v_h = ( L - L_h )L^iw̃^n. Note that, in the case w̃∈ C^♭̃-1(Ω) ∩ H^♭̃(Ω_h), 1≤ℓ̃≤ k+2, and 0≤ i ≤min(ℓ̃,♭̃)-1, we have w∈ C(Ω)∩ H^♭(Ω_h) and 1≤ℓ≤ k+2. Thus, (<ref>) holds and 
	
    ( L - L_h )L^iw̃≤    C_d L^iw̃_H^min(ℓ̃-i,♭̃-i)h^min(ℓ̃-i,♭̃-i)-1
    ≤    C_d w̃_H^min(ℓ̃,♭̃)h^min(ℓ̃,♭̃)-i-1.

	The proof is completed after removing the tilde symbol.


4. Since _h = (I-)L_h, adding and subtracting terms gives

    ( L- _h ) w
    =   ( L - L_h ) w +  L_h  w  
    
    =   ( L - L_h ) w + ( L_h  - L) w +  L w  
    
    =    (I-)( L - L_h ) w - (I-Π) Lw +  L w.

For 1≤ℓ≤ k+1, with (<ref>) in Lemma <ref>, it can be seen that

    (I-Π) Lw ≤(I-Π)Lw ≤ Cw_H^min(ℓ,♭)h^min(ℓ,♭)-1,
        L w = (I-)Lw≤ Cw_H^min(ℓ,♭)h^min(ℓ,♭)-1.
 
With (<ref>) in Lemma <ref>, we have 

    (I-)( L - L_h ) w≤( L - L_h ) w≤ C_d w_H^min(ℓ,♭)h^min(ℓ,♭)-1.

Applying the triangle inequality to (<ref>), and then substituting in the estimates above, we obtain (<ref>). 



 §.§ Proof of Lemma <ref>


 

[Invertibility of I+τ_h/2]
	(I+τ_h/2)^-1 is well-defined when λ≪1. Moreover, (I+τ_h/2)^-1w_h≤ Cw_h.


	Consider 
	
    (I+τ/2_h)v_h = w_h.

	We apply the triangle inequality to get
	
    v_h - τ/2_hv_h≤w_h.

	Recall that τ_h≤ C λ. When λ≪ 1, we have 1-τ_h/2≥ 1-Cλ/2 > 0. Hence
	
    v_h≤(1-τ/2_h)^-1w_h≤(1-Cλ/2)^-1w_h≤ Cw_h.
 
	When w_h = 0, (<ref>) forces v_h = 0, and thus v_h = 0. Therefore I+τ_h/2 is an injection. 	Note that I+τ_h/2: V_h^k → V_h^k is an operator on a finite-dimensional space. I + τ_h/2 is an injection implies that I + τ_h/2 is invertible. By rewriting v_h = (I+τ_h/2)^-1w_h in (<ref>), we obtain the desired estimate. 



[Commutation error of _h]
	If w ∈ C^♭-1(Ω)∩ H^♭(Ω_h) and ♭≥ 2, then 
	
    ( L- _h ) w≤ Cw_H^min(ℓ,♭)h^min(ℓ,♭)-1     1≤ℓ≤ k+1.



	By triangle inequality, we have
	
    ( L- _h ) w≤( L- _h ) w+ (I - ) Lw+(I - ) Lw

	Since w∈ C^♭-1(Ω)∩ H^♭(Ω_h) and ♭≥ 2, the regularity assumptions in both (<ref>) and (<ref>) of Lemma <ref> hold. Applying these estimates to (<ref>) would complete the proof of Lemma <ref>. 


Now we prove Lemma <ref>. 

	
	Since for w ∈ C^♭-1(Ω)∩ H^♭(Ω_h) with ♭≥ 2, we have w + τ L w/2 ∈ H^1(Ω_h). Thus (w + τ L w/2) is well-defined. By Lemma <ref>, w = (I+τ_h/2)^-1 (w + τ L w/2) is also well-defined.  
	
	To prove the approximation estimate (<ref>), note that
	
    (  - ) w =   (I+τ/2_h)^-1(I+τ/2 L)w -  w
    
    			=   (I+τ/2_h)^-1(I+τ/2 L)w - (I+τ/2_h)^-1( w+τ/2_h  w) 
    
    			=   τ/2(I+τ/2_h)^-1( L- _h )w.

	Applying Lemmas <ref> and <ref>, for all 1≤ℓ≤ k+1, one can get
	
    ( - ) w≤  Cτ( L- _h )w≤    C τw_H^min(ℓ,♭)h^min(ℓ,♭)-1
    ≤    C w_H^min(ℓ,♭)h^min(ℓ,♭).

	Using the triangle inequality and (<ref>), we have
	
    (I -)w≤(I -)w + 	( -)w≤ C w_H^min(ℓ,♭)h^min(ℓ,♭).





 §.§ Proof of Lemma <ref>


 

[Generalization of Lemma <ref>]
	Let α_i = (i!)^-1 for i = 0,1, ⋯,r. The operator inversion (∑_i = 1^s α_i (τ_h)^i-1)^-1 is well-defined if λ≪ 1. Moreover, 
	
    (∑_i = 1^s α_i (τ_h)^i-1)^-1 w_h≤ Cw_h.




	Consider (∑_i = 1^s α_i (τ_h)^i-1)v_h = w_h. Note that ∑_i = 1^s α_i (τ_h)^i-1 = I+∑_i = 2^s α_i (τ_h)^i-1. Then by the triangle inequality, we have 
	
    (1-∑_i = 2^s |α_i|( τ_h)^i-1)v_h≤w_h.
 
	Recall that τ_h≤ Cλ. By taking λ≪1, we can make 1-∑_i = 2^s |α_i| ( τ_h)^i-1>0 and hence v_h≤ Cw_h. As a result, w_h = 0 implies v_h = 0. Therefore, ∑_i = 1^s α_i (τ_h)^i-1 is an injection on a finite-dimensional space, and hence invertible. The norm estimate also follows from (<ref>) by taking v_h = (∑_i = 1^s α_i (τ_h)^i-1)^-1 w_h. 




[Generalization of Lemma <ref>]
	If w ∈ C^♭-1(Ω)∩ H^♭(Ω_h) and ♭≥ 2, then 
	
    ( L^i- _h^i ) w≤ Cw_H^min(ℓ,♭)h^min(ℓ ,♭)-i    ∀ 1≤ℓ≤ k+1,

	where 0≤ i≤min(ℓ,♭-1). 


	We prove by induction with respect to i. The case i = 0 is trivial, and the case i = 1 is proved in Lemma <ref>. Now we assume (<ref>) holds for i-1, and we prove the inequality holds for i with i≥ 2. Indeed, 
	
    ( L^i- _h^i) w = ( L^i-1 - _h^i-1) Lw + _h^i-1 ( Lw - _h  w).

	Noting that 2≤ i≤♭ -1, we have Lw∈ C^♭-2(Ω)∩ H^♭-1(Ω_h)⊂ C^i-1(Ω)∩ H^i(Ω_h)⊂ C^1(Ω)∩ H^2(Ω_h). Using the induction hypothesis with ℓ replaced by ℓ-1, we have 
	
    ( L^i-1 - _h^i-1) Lw≤    CL w_H^min(ℓ-1,♭-1)h^min(ℓ-1,♭-1)-(i-1)
    ≤     Cw_H^min(ℓ,♭)h^min(ℓ,♭)-i .

	With the fact _h≤L_h≤ Ch^-1 and Lemma <ref>, one can get
	
    _h^i-1 ( Lw - _h  w)≤    Ch^-(i-1)( L - _h ) w
    ≤     Cw_H^min(ℓ,♭)h^min(ℓ,♭)-i.

	Finally, by applying the triangle inequality to (<ref>), and invoking (<ref>) and (<ref>), we prove (<ref>) for i. By induction, (<ref>) holds for all 0≤ i≤min(ℓ,♭-1). 


Now we prove the approximation estimate of . 

	
    w -  w 
    
    			=   (∑_i = 1^s  α_i (τ_h)^i-1)^-1((∑_i = 1^q (i!)^-1 (τ L)^i-1)w + τ^q∑_i = q+1^s α_i (τ_h)^i-1-q L^q w) 
       -  w
    
    			=   (∑_i = 1^s  α_i (τ_h)^i-1)^-1((∑_i = 1^q (i!)^-1 (τ L)^i-1w) + τ^q∑_i = q+1^s α_i (τ_h)^i-1-q L^q w)
       -(∑_i = 1^s  α_i (τ_h)^i-1)^-1( (∑_i = 1^s  α_i (τ_h)^i-1) w).

			Combining the like terms, it yields
				
    w -  w 
    			=   (∑_i = 1^s  α_i (τ_h)^i-1)^-1(∑_i = 1^q (i!)^-1τ^i-1( L^i-1-_h^i-1) w)
       + (∑_i = 1^s  α_i (τ_h)^i-1)^-1(τ^q∑_i = q+1^s α_i (τ_h)^i-1-q(   L^q - _h^q)w ).

	Thus applying the estimate (∑_i = 1^s  α_i (τ_h)^i-1)^-1≤ C in Lemma <ref>, the triangle inequality, and τ_h≤λ, we have 
	
    (- ) w
    ≤   
    			C∑_i = 1^q  τ^i-1( L^i-1-_h^i-1) w
    			+Cτ^q∑_i = q+1^s   L^q w- _h^q w
    ≤    
    			C∑_i = 0^qτ^i( L^i-_h^i) w.

	We apply Lemma <ref> with ℓ = q≤min(r,k+1,♭-1). With the fact τ = λ h, one can get
	
    τ^i ( L^i-_h^i) w≤    C w_H^qτ^ih^q-i≤ Cw_H^qh^q  ∀ 0≤ i≤ q.

	Therefore, substituting (<ref>) into (<ref>), it gives
	
    ( - ) w≤ Cw_H^qh^q.

	Finally, recall that (I-) w≤ Cw_H^qh^q and apply the triangle inequality. We have
	
    (I - ) w≤(I - ) w + ( - ) w≤ Cw_H^qh^q.





§ CONCLUSIONS


In this paper, we perform stability analysis and optimal error estimates of a variant version of the widely used RKDG method. We focus on the linear advection equation in both one and two spatial dimensions and consider RKDG schemes of arbitrarily high order. We prove that even if we drop the highest-order polynomial modes at all inner RK stages, the resulting method will still maintain the same type of stability and optimal order of accuracy as those of the original RKDG method.

Through this paper, we establish theoretical groundwork to demonstrate the following fact: Within a multi-stage RK solver for time-dependent PDEs, it is feasible to employ cost-effective and low-order spatial discretizations at specific stages. Doing so allows for the creation of a more efficient fully discrete numerical solver that maintains the same stability and accuracy as the original method. The mathematical tools developed in this paper could be useful for the analysis of other numerical schemes beyond the classical method-of-lines framework and benefit the future development of numerical solvers for time-dependent PDEs.



§ ACKNOWLEDGMENT


The work of the author was partially supported by the NSF grant DMS-2208391. This work was also partially supported by the NSF grant DMS-1929284 while the author was in residence at the ICERM in Providence, RI, during the semester program “Numerical PDEs: Analysis, Algorithms, and Data Challenges". The author would like to thank the organizers, participants, and local faculty and staff for the wonderful workshops and many inspiring discussions. The author would also like to thank College of Arts and Sciences at The University of Alabama for the teaching release through the ASPIRE program, which made the visit to ICERM possible.





§ REMARKS ON ERROR ESTIMATES IN MULTI-STAGE FORM



In the proof of Theorems <ref> and <ref>, we formally rewrite the sdA-RKDG scheme into a one-step one-stage scheme and construct the desirable projection accordingly. However, in previous works for proving error estimates of the RKDG method, when constructing the reference function for deriving the error equation, one usually adheres to the multi-stage form of the RK method. Here, we provide remarks on the connections between the two approaches. To avoid unnecessary technicality, we will focus on the second-order case. We assume the exact solution is smooth and consider the 1D case only, for which we have L = L_h. We assume _h = L_h in this section so that the sdA-RKDG scheme formally coincides with the standard RKDG scheme. 

Analysis in multi-stage form. 
For the error estimate of RK2DG1 scheme, one starts with writing the exact solution u as

    u^(2) = u^n + τ/2 L u^n
    u^n+1 = u^n + τ Lu^(2) + τρ^n.

Next, we multiply a test function in V_h^k and then integrate by parts, which yields a similar equation as that of the RK2DG1 scheme. While in the strong form, it can be written as

    u^(2) =  u^n + τ/2 L u^n u^n+1 =  u^n + τ Lu^(2) + τρ^n.

After applying the relationship L= L_h, we have

    u^(2) =  u^n + τ/2 L_h u^n u^n+1 =  u^n + τ L_h  u^(2) + τρ^n.

However, if one directly works on u^n - u_h^n, it will lead to a suboptimal error estimate. The way to get around is to  consider ξ_h^n =  u^n - u_h^n. To this end, we add and subtract terms on both sides to get

	
    u^(2) =    u^n + τ/2 L_h u^n + τν_h^n,1 ,
     u^n+1 =    u^n + τ L_h  u^(2) + τν_h^n,2 + τρ^n.


Here ν_h^n,1=τ^-1(-Π)(u^(2)-u^n) and ν_h^n,2=τ^-1(-Π)(u^n+1-u^n). The ν_h^n,ℓ terms here are similar to the η̃ terms in <cit.>.

Then, one can subtract (<ref>) from (<ref>) to get the error equation

	
    ξ_h^(2) =   ξ_h^n + τ/2 L_h ξ_h + τν_h^n,1 ,
    ξ_h^n+1 =   ξ_h^n + τ L_h ξ_h^(2) + τν_h^n,2 + τρ^n.


Note that the error term ξ_h^n+1 satisfies an RK scheme for the linear advection equation together with a source term. Hence to derive optimal error estimates, one typically first derives a stability estimate of the RK scheme for the non-homogeneous equation u_t + β u_x = f(x,t). Then, one applies such stability estimate to (<ref>). Together with the fact ν_h^n,1, ν_h^n,2, and ρ^n are 𝒪(τ^2+h^2) terms, one can complete the optimal error estimate. 

Difference between  and . The first difference we notice is that in (<ref>), we use  to define the reference function, and in (<ref>), we use  to define the reference function. When _h = L_h, it yields

    w= (I+τ/2L_h)^-1(I + τ/2 L)w.

By simple algebraic manipulation and the relationship L_h  =  L, we have

    ( - )w =   (I+τ/2L_h)^-1( (I + τ/2 L)w -  (I + τ/2 L_h) w)
    
    =   τ/2(I+τ/2 L_h)^-1( L- L_h )w
    
    =   τ/2(I+τ/2 L_h)^-1( - )Lw.

One can show that, when τL_h≤ C λ≪ 1, (I+τL_h/2)^-1≤ C and thus

    ( - )w≤ Cτ( - )Lw≤ CLw_H^k+1τ h^k+1≤ Cw_H^k+2h^k+2.

Moreover, we have

    ( L - L_h)w = (L_h-L_h)w≤ CL_h( - )Lw≤ Cw_H^k+2h^k+1.

It can be seen that w is a 𝒪(h^k+2) order approximation of w and the operator  also satisfies the superconvergence property.  Although  and  are mathematically different, they are indeed very close and both satisfy the essential properties in Proposition <ref>. 

Compact form of (<ref>). For further comparison, we rewrite (<ref>) into a compact form

    u^n+1 =  u^n + τ L_h  u^n + τ^2/2L_h^2 u^n + τ(τ L_h ν_h^n,1 + ν_h^n,2 + ρ^n).

Note, if one replaces  in (<ref>) by , then this is of a similar form as that in (<ref>). In (<ref>), one needs to show that τ L_hν_h^n,1, ν_h^n,2, and ρ_h^n are 𝒪(τ^2+h^2) terms. While in (<ref>), one needs to show that ν_h^n, ζ_h^n, and ρ^n are 𝒪(τ^2+h^2) terms. In both cases, we write the reference function as u^n+1 = R_h,2 u^n + τ𝒪(τ^2+h^2), where = ,.  

Multi-stage form of (<ref>). When _h = L_h, one can rewrite (<ref>) into a multi-stage form.

	
    u^(2) =    u^n + τ/2L_h u^n,
     u^n+1 =    u^n + τ L_h  u^(2) + τ(ρ^n+ν_h^n + ζ_h^n).


Note, if one replaces  in (<ref>) by , this is of a similar form as that in (<ref>). Only the source terms are distributed differently across stages. 


abbrv

