ANY ORDER SPECTRAL VOLUME METHODS FOR DIFFUSION EQUATIONS USING THE LOCAL DISCONTINUOUS GALERKIN FORMULATION

. In this paper, we present and study two spectral volume (SV) schemes of arbitrary order for diffusion equations by using the local discontinuous Galerkin formulation to discretize the viscous flux. The basic idea of the scheme is to rewrite the diffusion equation into an equivalent first-order system first, and then use the SV method to solve the system. The SV scheme is designed with control volumes constructed by using the Gauss points or Radau points in subintervals of the underlying meshes, which leads to two SV schemes referred to as LSV and RSV schemes, respectively. The stability analysis for the linear diffusion equations based on alternating fluxes are provided, and optimal error estimates are established for both the exact solution and the auxiliary variable. Furthermore, a rigorous mathematical proof are given to demonstrate that the proposed RSV method is identical to the standard LDG method when applied to constant diffusion problems. Numerical experiments are presented to demonstrate the stability, accuracy and performance of the two SV schemes for both linear and nonlinear diffusion equations


Introduction
The spectral volume (SV) method is a class of high order Godunov-type finite volume method [10], which has been under development for several decades and is considered to be the current state-of-the-art for the numerical solution of hyperbolic conservation laws.The SV method was first proposed by Wang, Liu et al. and their collaborators for hyperbolic conservation laws on unstructured grids (see, e.g., [29][30][31]).Since then it attracted a lot of attention and was successfully implemented for solving various PDEs such as Euler equation [32], Navier-Stokes equations [22,23], and 3D Maxwell equations [17].Similar to the discontinuous Galerkin (DG) methods [5,6,8,11], SV methods adopt completely discontinuous functions as solution space and have many desired advantages such as the allowance of hanging nodes, compact stencils, easy ℎ adaptivity, high parallel efficiency, and so on.Comparisons between the SV method and other numerical methods (e.g., DG, spectral difference) in terms of accuracy and stability have also been conducted in the literature, see [21,25,34].
For equations containing higher order spatial derivatives, such as diffusion equations, however, SV methods can not be directly applied due to the discontinuous solution space at the element interfaces, which is not regular
While at the interface of each element (i.e.,  − 1 2 ,  ∈ Z  ), we choose the alternating fluxes.That is, B and qℎ are taken from opposite sides.To be more precise, We end with this subsection the discussion on specific SV schemes.By choosing different points   , s ,  ∈ Z  , we can get different SV schemes.Note that the stability and accuracy of the SV scheme are heavily dependent on the construction of control volume, i.e., the choice of   , s ,  ∈ Z  .In this paper, we restrict our discussion to the following two SV schemes.
GSV: both {  }  =1 and {s  }  =1 are chosen as Gauss-Legendre points, i.e.,   = s ,  ∈ Z  are  zeros of the Legendre polynomial   of degree .Remark 2.1.The purpose of introducing the two sets {  } +1 =0 and {s  } +1 =0 (separately corresponding to  , and x, ) is to ensure the stability of our numerical scheme (2.4).Actually, due to the fact that the numerical fluxes B and qℎ in (2.4) are taken from opposite sides, the associated points   and s (which is not necessary to be the same) should be carefully taken to match the choice of the fluxes for stability.The choice of {  }  =1 and {s  }  =1 in RSV scheme and our later theoretical analysis will demonstrate this point.
Remark 2.2.For convection diffusion equations, i.e.,   + ()  = (  )  , we can also design the corresponding GSV and RSV schemes by using the LDG formulation with the alternating flux for diffusion term, and the upwind flux for the convection term.For simplicity and clarity, we focus our attention on the diffusion equations in this paper.The same argument can also be applied to convection-diffusion equations.

SV methods as a Petrov-Galerkin method
In this subsection, we reformulate the SV scheme (2.4) into its equivalent Petrov-Galerkin form, and then establish the stability of SV scheme under the the framework of Glerkin method.
We begin with the construction of two types of control volumes, which are defined by With these control volumes V , , V, , we define the piecewise constant function space associated with V , , V, as Obviously, any function  * ∈  ℎ , w* ∈ Vℎ can be represented as where  * , , w* , , (, ) ∈ Z  × Z 0  are coefficients as functions of ,   ,  ⊂ [, ] is the characteristic function defined as   = 1 in  and   = 0 otherwise.Define the so-called broken Sobolev space as follows: For all  ∈ Z  and any (, ,  * , w* ) Then it is easy to check that the SV solution ( ℎ ,  ℎ ) of (2.4) satisfies or equivalently, Conversely, if  ℎ ,  ℎ ∈  ℎ satisfy (2.10) or (2.11), then by choosing  * =  V, , w* =  V, , we find that  ℎ ,  ℎ satisfy (2.4).In other words, the SV method (2.4) is equivalent to the Petrov-Galerkin method (2.10) or (2.11).

Stability analysis for linear diffusion problems
In this subsection, we investigate the stability of the two SV schemes for the linear diffusion problem (2.2), i.e.,  = ().In this case,  = √︀ () is continuous and thus the numerical fluxes β =  in (2.5) or (2.6).We begin with the introduction of the special transformation from the trial space  ℎ to the test spaces  ℎ , Vℎ .

Transformation from the trial space to the test space
Since the transformation we defined in this subsection is closely related to some quadratures points and their associated weights.We introduce some numerical quadratures and quadrature errors first.
Given any  ∈  1 ([−1, 1]), suppose   , s ,  ∈ Z 0 +1 are points of some numerical quadrature to calculate ∫︀ 1 −1  () d and   , Ā are associated weights.Define If   are taken as the Gauss points, then the above quadrature is referred as the standard Gauss-Legendre quadrature with ,  = 1, . . .,  and  0 =  +1 = 0. Similarly, if   are taken as the right Radau points, we have  0 = 0 and call the quadrature right Radau numerical quadrature.If   are taken as the left Radau points, we have  +1 = 0. Note that the -point Gauss quadrature and ( + 1)-point right/left Radau quadrature is exact for polynomials of degree not more than 2 − 1 and 2, respectively.
Define the residual of the numerical quadrature by By a scaling from [−1, 1] to   , we get the residual of the numerical quadrature in each interval   .That is, where  , = ℎ 2   , Ā, = ℎ 2 Ā for all (, ) ∈ Z  × Z 0 +1 .Now we are ready to present the transformation from the trial space to the test space.for all  ∈  ℎ , we define a transformation ℱ :  ℎ →  ℎ by where the  + 1 coefficients  * , ,  ∈ Z 0  are given as By a direct calculation, we have Consequently, we use the inverse inequality to get For any function ,  ∈ ℋ ℎ , we define Due to the special transformation ℱ, F, we have the following relationship between inner product and the discrete inner product: )︁ .
By using (2.14) and (2.15) and the integration by parts, we get Similarly, there holds Here for any  ,   ( ), R ( ) denote the numerical quadrature errors defined in (2.12).

𝐿 2 stability
We first study the bilinear terms  Then for all (, ) 2 ) and  = (, ) is defined by (2.3) with  replaced by .
Proof.In each element   , if not otherwise stated, we use the notation | + 1 2 , | − 1 2 to denote the value of the function  at the boundary points  + 1  2 ,  − 1 2 , respectively.That is, Recalling the definition of  1  in (2.7) and using (2.14) and (2.15), we get where in the last step, we have used (2.20) and We next estimate .Noticing that  0 =  +1 = 0 for LSV schemes, we have  = 0.As for RSV schemes, if the fluxes are chosen as (2.5) (i.e., p =  − ), then  , ,  ∈ Z  are taken as  interior right Radau points and thus, Similarly, if the fluxes are chosen as (2.6) (i.e., p =  + ), then  , ,  ∈ Z  are chosen as interior left Radau points, which yields Then for both the LSV and RSV, we have  = 0. Consequently, where in the last step, we have used the integration by parts.Then (2.22)In other words, (2.25) holds true for both GSV and RRSV, which indicates (, ℱ) and (, F) are both equivalent to the standard  2 -norm.Now we are ready to present the  2 stability for both GSV and RSV schemes.
Theorem 2.4.Let ( ℎ ,  ℎ ) be the solution of (2.11).Then both GSV and RSV are  2 stable, i.e., Proof.By choosing  = , w =  in (2.22), (2.23) and using the numerical quadrature error, we get for all ,  ∈  ℎ and  * = ℱ,  * = F that Recalling the definition of  in (2.3) and using the continuity of the linear variable function  = (), we have Then summing up all  from 1 to  yields We next estimate the numerical quadrature errors   (  ) and R (  ).Given any function  , we define    |  ∈ P  the Lagrange interpolation function of  with the interpolation point  , , 0 ≤  ≤  + 1, here  =  for GSV and  =  + 1 for the RSV.Since   ( ) = 0,  ∈ P 2−1 for GSV and   ( ) = 0,  ∈ P 2 for RSV, then we have for both GSV and RSV that The above identity is also valid for R (   ).Consequently, Here   denotes the cell average of  in the interval   .Similarly, there holds Using the inverse inequality and the fact that (, )  =   , we have for all ,  ∈  ℎ that Substituting the above two inequalities into (2.27),we get for both LSV and RSV that Especially, we choose (, ) = ( ℎ ,  ℎ ) in the above inequality and use (2.11), (2.25) and Cauchy-Schwarz inequality and then obtain d d Then (2.26) follows from the Gronwall inequality and the equivalence (2.25).
We end with this subsection the discussion on the relationship between the proposed SV schemes and the LDG schemes.Recall the bilinear form of the LDG schemes for (2.2), which is defined by The LDG method for (2.2) is: In light of (2.22), (2.23) and (2.20), we easily obtain for all , w ∈  ℎ (, ;  * , w* ) =  LDG (, ; , w) where  * = ℱ, w* = F w.As we may observe, the bilinear forms of the proposed LSV and RSV schemes are equivalent to those of the LDG schemes, up to some Gauss or Radau numerical quadrature errors.Especially, for constant coefficients problem (i.e.,  = ), we have Note that the left/right Radau numerical quadrature is exact for polynomial of degree not more than 2.Then for RSV, R Consequently, (, ;  * , w* ) =  LDG (, ; , w), ∀,  ∈  ℎ , which indicates that the RSV method is identical to the LDG method for constant diffusion problems with the source function (, ) ∈ P 2 .

Error estimates for linear diffusion problems
This subsection is dedicated to the error estimates for linear diffusion problems, i.e.,  = ().To study the convergence of SV methods, we first introduce some interpolation functions.For any function  ∈ ℋ ℎ , we denote by  + ℎ ,  − ℎ  ∈  ℎ the standard Lagrange interpolation of , which satisfy the following  + 1 conditions in each element   : Similarly, we can define the interpolation functions Ī+ ℎ , Ī− ℎ  ∈  ℎ of , with the interpolation points x, substituting of  , in the above definition.By the standard approximation theory of the interpolation function, we have Due to the special construction of the interpolation function, we have for both fluxes (2.5) and (2.6) that Theorem 2.5.Let  be the solution of (2.1) satisfying  ∈  +2 ,   ∈  +1 , and ( ℎ ,  ℎ ) be the solution of (2.11) with the initial solution chosen as  ℎ (, 0) =  ℎ  0 .Then for both the flux choices (2.5) and (2.6) and both GSV and RSV, By taking (, ) = (  ,   ) in (2.29) and the orthogonality ( −  ℎ ,  −  ℎ ;  * , w* ) = 0 for all , w ∈  ℎ , we get Recalling the definitions of  1  in (2.7) and using (2.32), we easily get for all  * ∈  ℎ , w* ∈ Vℎ that On the other hand, by denoting | , = (x , ) and using the continuity of () and (2.3), we have In light of (2.32), we easily obtain Consequently, we have from (2.16) and the second inequality of (2.19) that Substituting the above inequality into (2.34) and using the Cauchy-Schwarz inequality gives where  is a constant independent of the mesh size ℎ.Integrating the above inequality from 0 to  and using (2.25) and the Gronwall inequality, we get Then the desired result follows from the approximation properties of ( ℎ ,  ℎ ).
Remark 2.6.The stability result and optimal error estimates established in Theorems 2.4 and 2.5 are mainly based on linear diffusion equations.Although our numerical experiments show that the stability result and optimal error estimtates hold true for nonlinear equations, its theoretical analysis is still lacking.Just as demonstrated in Theorem 2.4, the numerical quadrature errors appeared in the right hand side of (2.27) can not be bounded by the  2 norms of  and  for nonlinear equations, which leads to the difficulty in stability.Choosing suitable SV numerical fluxes may provide a new approach to study the stability of SV method for nonlinear equations.The theoretical analysis of SV methods for nonlinear problems is interesting yet difficult, and it deserves a separate study.

SV method for multidimensional linear diffusion problems
In this section, we present the SV scheme for multidimensional linear diffusion problems and investigate its stability and convergence.The analysis for multidimensional nonlinear diffusion problems deserves a separate study and thus we skip it here.To simplify our analysis and make the idea clear, we will use the two dimensional linear diffusion problems as a model.The methodology we adopt can also be applied to linear diffusion equation in higher dimensions, e.g., three dimension.
We consider the SV method for the following two dimensional linear diffusion problems where both  0 and  := (, ) ≥ 0 are smooth.For simplicity, we still consider the periodic boundary condition.We rewrite the above equation into its equivalent system:

SV schemes
Denote by  ℎ the rectangular partition of Ω.That is, For any  ∈  ℎ , we denote by ℎ   , ℎ   the lengths of -and -directional edges of  , respectively.ℎ is the maximal length of all edges, and . We assume that the mesh  ℎ is quasi-uniform in the sense that there exist constants  1 ,  2 > 0 such that Denote by Q  (, ) = P  () × P  () the tensor product bi- polynomial space, and define the finite element space Let −1 =  0 <  1 < . . .<   <  +1 = 1, and −1 = s0 < s1 < . . .< s < s+1 = 1, and define Assume that   is the affine mapping from τ to  , and Two types of control volumes can be constructed by using the above ( + 2) 2 points, i.e., The SV schemes for (3.2) is to find a  ℎ ∈  ℎ , q ℎ ∈ U ℎ such that for all (, ) Here qℎ , B denote the numerical fluxes, n 0 = (1, 1), and n is the outward normal vector.We still take the alternating flux as our numerical fluxes, i.e., Similar to the 1D case, we consider two classes of SV schemes here, i.e., GSV schemes in which   , s are chosen as Gauss points, and RSV schemes in which   , s are chosen as right or left Radau points dependent on the choice of numerical fluxes.
We next reformulate the SV scheme (3.3) into its equivalent Petorv-Glerkin form.To this end, we first define two piecewise test spaces as follows.

Transformation from the trial space to the test space
Similar to the 1D problem, we need to define a special transformation from the trial space to the test space, which is of great importance in our stability analysis.
Similarly, we can define a transformation F from the vector trial space U ℎ to the vector test space V ℎ with the points   , replaced by ḡ , , and there also holds In each element  =    ×    , we define Since the transformations ℱ and F share the same properties, we next only consider the property of ℱ and the same argument can be applied to F.
The proof of Lemma 3.1 is given in Appendix A.
where ℰ  is define in (3.16) and Consequently, for all v = ( 1 ,  2 ) ∈ ℋ ℎ × ℋ ℎ , Proof.We only prove (3.17) as the same argument can be applied to (3.18).Let  * = ℱ.By (3.12), we have where  1 ,  2 are defined in (3.14) and (3.15).Using (3.13) and the integration by parts, we get Recalling the definition of  1 ,  2 in (3.Then the desire result (3.17Here  := (, , ) is defined in (3.2).For simplicity, we adopt the notation () instead of (, , ) if no confusion arises.Then the LDG method for (3.2) is: find In light of (3.5), we get By using (3.13) and (3.20), we have for all p = ( 1 ,  2 ) Similarly, there holds for all w = ( 1 ,  2 ) where Ē , Ē  ,  = 1, 2 are given in (3.16) and (3.19) with R  , R  substituting of    ,    , respectively.As indicated by (3.21) and (3.22), the proposed GSV and RSV schemes are equivalent to the counterpart LDG schemes up to some Gauss or Radau numerical quadrature errors.Especially, for constant coefficients problems, i.e.,  =  in (3.2), using the left or right Radau quadrature is exact for polynomials of degree not more than 2, we get for all  ∈  ℎ , p = ( which indicates that the RSV method is exactly the same as the LDG method when applied to constant diffusion problems.

Stability
To study the stability of the SV schemes, we first estimate the numerical quadrature errors appeared in (3.21) and (3.22).Lemma 3.3.For any ,  ∈  ℎ , let ℰ   ((), ),  = 1, 2 be defined in (3.19) with ( B, ) = ( + ,  − ) or ( B, ) = ( − ,  + ), and ℰ  ((), ) be defined in (3.16).Denote Then for both GSV and RSV, The proof of Lemma 3.3 is given in the Appendix A. Similar to the 1D case, we define (, ℱ) = ∑︁ In light of (3.13), we have (, ℱ) = (, ) + ∑︁ Recalling the definition of ℰ  in (3.16), we have ℰ(, ) = 0, ,  ∈  ℎ for the RSV scheme.As for the GSV scheme, by using the error of Gauss numerical quadrature in [−1, 1], there holds for some  ∈ (−1, 1) By a scaling from [−1, 1] to    and    and using the fact that    ,     are both constants for all  ∈  ℎ , we get for any  =    ×    that Here  is a positive constant only dependent on .Then a direct calculation from the inverse inequality yields Consequently, for both GSV and RSV, Now we are ready to present the stability result for both RSV and GSV schemes.

Error estimates
This subsection is dedicated to error estimates of the SV method for 2D equations in rectangular meshes.To this end, we first introduce a suitable projection similar to the one-dimensional case, which is defined by where the superscripts indicate the application of the one-dimensional operator   ℎ with respect to the variable  with (  ℎ ,   ℎ q) given in (2.31).The operator   ℎ  follows the same definition along the  direction.The definition of (  ℎ ,   ℎ q) ensures that for both alternating fluxes, For any function  ∈ ℋ ℎ , we define We see that    and    is continuous about  and , respectively.By the approximation theory, we have for all  ≤  ≤  + 1 that Furthermore, a direct calculation yields that Now we are ready to present the error estimates for the SV method.
Since the exact solution (, q) also satisfy (3.6), we take (, p) = (  ,  q ) in (3.28) and use the orthogonality to get On the other hand, recalling the definitions of  1  ,  2  in (3.5), we get for all Using the error decomposition (3.33), we get In light of the properties of   in (3.31) and the fact   is continuous about , we have Similarly, there holds Substituting the last three equations into  1  yields Following the same argument, we have for all which yields, together with (3.35) and (3.26) that Integrating the above inequality from 0 to , using the Gronwall inequality and the  2 norm equivalence (3.25), we have Then the desired result follows from standard approximation property of the interpolation function and the inequality ‖  ‖ +1 ‖‖ +3 due to   = −∇ • (∇).

Numerical results
In this section, we present some numerical examples to test the efficiency and accuracy of the SV method for solving the one and two dimensional diffusion equations.In our experiments, we implement the GSV and RSV numerical schemes with  = 1, 2, 3 and the fluxes choice (2.5).We shall compute the  2 error for both the exact solution  and the auxiliary variable  or q, which is denoted by ‖  ‖ 0 and ‖  ‖ 0 or (‖ q ‖ 0 ).To diminish the time discretization error, we use the fourth order Runge-Kutta method with time step △ = 0.005ℎ 2 .
Example 1.We solve the following problem The source function (, ) is chosen such that the exact solution is (, ) =  − sin().
We consider two cases: the linear equation with variable coefficient  = 1 + sin 2  and the nonlinear equation with  = sin 2 .We use non-uniform meshes of  elements in our experiment, which are obtained by randomly and independently perturbing each node of a uniform mesh by up to some percentage.That is, where () returns a uniformly distributed random number in (0, 1).Table 1 shows the  2 errors and the corresponding convergence rate calculated from the GSV and RSV schemes for  = 1 + sin 2  and  = 1, 2, 3. We observe that, for both GSV and RSV, the convergence rates for the errors ‖  ‖ 0 and ‖  ‖ 0 can achieve (ℎ +1 ).These numerical results verify the theoretical findings (2.33) in Theorem 2.5.
Listed in Table 2 are  2 errors and the rates of convergence of two SV schemes for the nonlinear problem with  = sin 2 .Similar to the linear equation, we again observe an optimal convergence order of (ℎ +1 ) for both the errors ‖  ‖ 0 and ‖  ‖ 0 , which indicates that the error estimate (2.33) also holds true for the nonlinear problems.Furthermore, as demonstrated from Tables 1 and 2, the GSV method seems to have larger  2 errors than the counterpart RSV method on the same meshes.with the periodic boundary condition (0, , ) = (2, , ) and (, 0, ) = (, 2, ).
We compute the numerical solution at  = 0.1.The computational results are given in Tables 3 and 4, from which we observe a convergence rate of (ℎ +1 ) for both ‖  ‖ 0 and ‖ q ‖ 0 in cases  = 1 and  = 1 + sin 2 .These results are consistent with the error estimates (3.34) in Theorem 3.5, which indicates that the error bound given in (3.34) is sharp.

Concluding remarks
In this work, we present and develop two classes of arbitrary high order SV schemes for solving diffusion equations by using the LDG formulation to discretize the viscous flux.We prove that both the proposed GSV and RSV schemes are energy stable when alternating fluxes are used.Optimal error estimates for both the exact solution and the auxiliary variable are also established.The comparison between the proposed SV method and the standard LDG method is given.A interesting discovery is that: the RSV is identical to the LDG method when applied to constant diffusion equations.Generally speaking, the SV method shares similar behaviors on stability, accuracy and convergence, from the perspective of theoretical analysis.While in terms of computation complexity, LDG requires both volume and surface integrations.In contrast, SV requires only surface integrations and thus requires less computational efforts or CPU time to reach the same accuracy than the counterpart LDG method.This advantage in computation efficiency makes SV method more appealing for solving some comparatively complicated nonlinear diffusion equations, i.e., the porous medium equation (PME).The current study is our first attempt to investigate the difference between the LDG method and SV method, more theoretical and numerical investigations are needed to find advantages of SV methods.Our ongoing works include the study of SV methods for nonlinear diffusion problems and their superconvergence properties.Appendix A.
We begin with some preliminaries.For any function  ∈  ℎ , let  * = ℱ be defined in (3.9) with the coefficients given by (3.10).Then a direct calculation from (3.10) yields )︀ .
By using (3.7), (3.8) and the integration by parts, Following the same argument, we can prove that Substituting the last two equations into (A.2) yields the desired result (3.12).We next prove (3.13). where )︁ .
By using (3.10) and (A.1), we have This finishes the proof the first inequality of (3.23).The second inequality of (3.23) follows from the same argument and thus we omit it here.

Table 2 .
2 errors and corresponding convergence rates of the SV method for Example 1 with  = sin 2  and  = 1, 2, 3.

Table 3 .
2 errors and corresponding convergence rates of the SV method for Example 2 with  = 1 and  = 1, 2, 3.Nonuniform meshes of  ×  rectangles are obtained by randomly and independently perturbing each node in the  and  axes of a uniform mesh by up to some percentage.That is, for all (, ) ∈ Z 0  × Z 0  ,