A sparse grid discontinuous Galerkin method for the high-dimensional Helmholtz equation with variable coefficients
aa r X i v : . [ m a t h . NA ] J un A SPARSE GRID DISCONTINUOUS GALERKIN METHOD FOR THEHIGH-DIMENSIONAL HELMHOLTZ EQUATION WITH VARIABLE COEFFICIENTS
WEI GUO ∗ Abstract.
The simulation of high-dimensional problems with manageable computational resource represents a long stand-ing challenge. In a series of our recent work [25, 17, 18, 24], a class of sparse grid DG methods has been formulated for solvingvarious types of partial differential equations in high dimensions. By making use of the multiwavelet tensor-product bases onsparse grids in conjunction with the standard DG weak formulation, such a novel method is able to significantly reduce thecomputation and storage cost compared with full grid DG counterpart, while not compromising accuracy much for sufficientlysmooth solutions. In this paper, we consider the high-dimensional Helmholtz equation with variable coefficients and demonstratethat for such a problem the efficiency of the sparse grid DG method can be further enhanced by exploring a semi-orthogonalityproperty associated with the multiwavelet bases, motivated by the work [21, 22, 19]. The detailed convergence analysis showsthat the modified sparse grid DG method attains the same order accuracy, but the resulting stiffness matrix is much sparserthan that by the original method, leading to extra computational savings. Numerical tests up to six dimensions are providedto verify the analysis.
Keywords: discontinuous Galerkin method; Helmholtz equation; variable coefficients; semi-orthogonality;sparse grid; high dimensions.
1. Introduction.
Discontinuous Galerkin (DG) methods that make use of discontinuous functions asapproximations have been extensively studied for solving partial differential equations (PDEs) over the lastfew decades and become rather mature for various applications [13]. It is generally understood that DGmethods are more flexible compared with continuous finite element methods due to the lack of continuityrequirement, and thus enjoying many attractive properties. However, the situation changes when the di-mension of the underlying problem becomes large. Traditional grid-based methods including DG methodssuffer from the curse of dimensionality [6], which describes the scenario that the complexity and memorystorage of an algorithm grows exponentially with the dimension of the underly problem for a given level ofaccuracy. Moreover, compared with other high order schemes, DG methods often require more degrees offreedom (DOF), and hence the effect of the curse of dimensionality is more significant, making such methodsuncompetitive for high-dimensional calculations.Driven by the need for a method that is able to retain the attractive properties of DG methods and yetrequires feasible computational and storage cost for high-dimensional simulations, based on the sparse gridapproach [27, 12, 15], a class of novel DG methods has been developed for solving various types of high-dimensional PDEs in a series of work [25, 17, 18, 24, 23, 20]. For such a method, the underlying DG solutionis represented by a set of orthonormal tensor-product multiwavelet bases [1, 25, 17]. Meanwhile, followingthe standard sparse grid philosophy, rather than including all the anisotropic tensor-product bases, only theones with significant contribution to the approximation accuracy are chosen, leading to immense reduction incost complexity when the dimension d is large. In particular, the number of DOF scales as O ( h − | log h | d − )instead of O ( h − d ), where h denotes the mesh size in each direction; while the approximation property ofthe proposed sparse grid DG method is proven only slightly deteriorated for smooth solutions through boththeoretical and numerical verifications, see [25, 17].In this paper, we are concerned with the Helmholtz equation with variable coefficients, which arisesin many applications in science and engineering such as the representation of the solution of Schr¨odingerequation [22]. In [25], we developed an efficient high order sparse grid DG method based on the interiorpenalty DG (IPDG) weak formulation [5, 14, 26, 2, 3] for solving general elliptic equations, and such a methodcan be directly applied to simulate the problem. Due to the reduction of DOF, the resulting algebraic linearsystem enjoys a much smaller dimension than that by the traditional DG method, leading to computationalsavings in high dimensions. However, the linear system also becomes denser especially in the variable-coefficient case because of the hierarchical nature of the multiwavelet bases, and thus impeding the efficiencyadvantage of the sparse grid approach. To address the issue, we would like to explore the semi-orthogonalityproperty associated with the orthonormal multiwavelet bases, which can be thought of as a special type oforthogonality of the bases related to the underlying bilinear formulation. The idea of making use of semi-orthogonality on sparse grids was proposed in [21, 22, 19] in the continuous Ritz-Galerkin discretization ∗ Department of Mathematics and Statistics, Texas Tech University, Lubbock, TX, 70409. Research is supported by NSFgrants DMS-1620047, DMS-1830838. [email protected] ramework. In particular, by using the prewavelet as the bases, the standard bilinear formulation is carefullymodified according to the semi-orthogonality property, aiming to sparsify the stiffness matrix for the variable-coefficient equations. The analysis in [22] further shows that such a modification preserves the accuracy of theoriginal method under some extra smoothness requirement of the coefficient. In this work, we demonstratethat the modification of the bilinear form based on semi-orthogonality is also effective for our sparse gridIPDG method with multiwavelet bases. In the theoretical aspect, we refine the analysis in [22] and show thatthe error incurred by modifying the bilinear form is one order higher than the projection error for sufficientlysmooth problems, making the modified method produce virtually the same numerical result as the originalsparse grid IPDG method; while the associated linear system is much sparser, offering extra computationaland storage savings. Furthermore, under the sparse grid DG framework, the proposed methodology has thepotential to be extended to many other variable-coefficient PDEs.The rest of this paper is organized as follows. In Section 2, we review the fundamentals of the IPDGon sparse grids for solving the Helmholtz equation developed in [25, 17]. In Section 3, we formulate themodified scheme based on the semi-orthogonality property and perform an error analysis. In Section 4,numerical results in multi-dimensions (up to d = 6) are provided to validate the accuracy and performanceof the method. Conclusions and future work are given in Section 5.
2. IPDG Method on Sparse Grids for the Helmholtz Equation.
In this section, we review theconstruction together with several key theoretical results of the IPDG method on sparse grids for solvingthe following Helmholtz equation, − ∆ u + cu = f in Ω = [0 , d , (2.1) u =0 on ∂ Ω , where c is a non-negative function and the domain is the d -dimensional unit box. In this subsection, we briefly review the recentdevelopment on the DG finite element approximation space on sparse grids in [25, 17], which plays a keyrole in dimension reduction for the proposed DG method.We first introduce the hierarchical decomposition of DG approximation space in one dimension. Withoutloss of generality, we consider the domain Ω = [0 ,
1] and define the l -th level grid Ω l , composed of non-overlapping intervals I jl = (2 − l j, − l ( j + 1)], j = 0 , . . . , l − h l = 2 − l . Denoteby V kl = { v : v ∈ P k ( I jl ) , ∀ j = 0 , . . . , l − } the collection of piecewise polynomials of degree at most k defined on grid Ω l . Note that there exists anested structure with respect to different values of l : V k ⊂ V k ⊂ V k ⊂ V k ⊂ · · · , which allows us to define the increment space W kl , l = 1 , , . . . as the orthogonal complement of V kl − in V kl with respect to the standard L inner product on Ω, i.e. V kl − ⊕ W kl = V kl , W kl ⊥ V kl − . By letting W k := V k , which consists of all polynomials of up to degree k on the entire domain [0 , V kN , N ≥ N as V kN = M ≤ l ≤ N W kl . (2.2)In light of (2.2), we are able to represent a function in V kN by making use of the orthonormal multiwaveletbases of space W kl developed in [1]. We denote such orthonormal bases of W kl as v ji,l ( x ) , i = 1 , . . . , k + 1 , j = 0 , . . . , max(0 , l − − . he details about the construction of the multiwavelet bases can be found in [1, 25].It is worth noting that space W kl can be constructed by means of projection operators as well, whichwill be helpful in the analysis below. Denote P kl as the standard L projection operator from L (0 ,
1) onto V kl , i.e., V kl = P kl ( L (0 , Q kl := ( P kl − P kl − , if l ≥ ,P k , if l = 0 . Evidently, W kl = Q kl ( L (0 , . To summarize, we have obtained the hierarchical decomposition of the piecewise polynomial space V kN aswell as the corresponding projection operator P kN : V kN = M ≤ l ≤ N W kl , P kN = X ≤ l ≤ N Q kl . Then, we review the construction of approximation spaces in multi-dimensions. Consider the domain x = ( x , . . . , x d ) ∈ Ω = [0 , d . We first recall some notation on the norms and operations of multi-indices in N d , where N denotes the set of nonnegative integers. For α = ( α , . . . , α d ) ∈ N d , we define the l and l ∞ norms | α | := d X m =1 α m , | α | ∞ := max ≤ m ≤ d α m , the component-wise arithmetic operations α ± β := ( α ± β , . . . , α d ± β d ) , α := (2 α , . . . , α d ) , the relational operator α ≤ β ⇔ α m ≤ β m , ∀ m, and max( α , β ) := (max( α , β ) , . . . , max( α d , β d )) , min( α , β ) := (min( α , β ) , . . . , min( α d , β d )) . Denote := (0 , . . . , ∈ N d and := (1 , . . . , ∈ N d .Given a multi-index l = ( l , . . . , l d ) ∈ N d indicating the levels of the mesh in a multivariate sense, wedefine the grid Ω l consisting of non-overlapping elementary cells I jl = { x : x m ∈ (2 − l m j m , − l m ( j m +1)) , m =1 , . . . , d } , ≤ j ≤ l − . On grid Ω l , we define the tensor-product increment space W k l := W kl ,x × W kl ,x · · · × W kl d ,x d , where W kl m ,x m corresponds to the one-dimensional (1D) increment space in the m -th direction. Based on the1D hierarchical decomposition (2.2), we yield V kN := V kN,x × V kN,x · · · × V kN,x d = M | l | ∞ ≤ N W k l , where V kN is indeed the traditional tensor-product DG approximation space on grid Ω N := Ω ( N,...,N ) , i.e.,having 2 N cells in each dimension. V kN will be called the full grid space hereafter.Define the index set Θ l = { j ∈ N : ≤ j ≤ max( , l − − ) } . The basis functions for W k l areconstructed by a tensor product of the 1D multiwavelet bases v ji , l ( x ) := d Y m =1 v j m i m ,l m ( x m ) , j ∈ Θ l , ≤ i ≤ k + , (2.3) here k := ( k, . . . , k ) ∈ N d . Note that they form a set of orthonormal bases due to the property of the 1Dbases. The sparse grid DG finite element approximation space is defined as b V kN := M | l | ≤ N W k l . Evidently, b V kN is a subset of V kN . More importantly, its number of DOF scales as O (( k + 1) d N N d − ),as opposed to exponential dependence on N d for the full grid space V kN , see [25]. This is the key forcomputational savings in high dimensions.The standard L projection operator onto space W kl can be naturally written as Q k l := Q kl ,x ⊗ · · · ⊗ Q kl d ,x d , where Q kl m ,x m denotes the 1D projection operator Q kl m in the m -th dimension. Consequently, we write the L projection operator onto the sparse approximation space b V kN as b P kN := X | l | ≤ N Q k l . (2.4) To formulate the IPDG method on sparse grids, we start withintroducing some standard notation about jumps and averages for piecewise functions defined on grid Ω N .Let Γ := S T ∈ Ω N ∂ T be the union of the boundaries for all the elements in Ω N and S (Γ) := Π T ∈ Ω N L ( ∂T ) bethe set of L functions defined on Γ. For any q ∈ S (Γ) and q ∈ [ S (Γ)] d , their averages { q } , { q } and jumps[ q ] , [ q ] on the interior edges are defined as follows,[ q ] = q − n − + q + n + , { q } = 12 ( q − + q + ) , [ q ] = q − · n − + q + · n + , { q } = 12 ( q − + q + ) , where n denotes the unit normal, and ‘ − ’ and ‘+’ represent the directions of the underlying vector pointingto interior and exterior at an element edge e , respectively. If e is a boundary edge, then we let[ q ] = q n , { q } = q , where n is the outward unit normal.The sparse grid IPDG scheme for (2.1) developed in [25] is defined as follows. We look for u h ∈ b V kN ,such that B ( u h , v ) = L ( v ) , ∀ v ∈ b V kN (2.5)where B ( w, z ) = Z Ω ∇ w · ∇ z d x + Z Ω cwz d x − X e ∈ Γ Z e {∇ w } · [ z ] ds − X e ∈ Γ Z e {∇ z } · [ w ] ds + X e ∈ Γ σh Z e [ w ] · [ z ] ds, (2.6)and L ( v ) = Z Ω f v d x ds, (2.7)where σ is a positive penalty parameter, h = 2 − N is the uniform mesh size in each dimension. Note that thebilinear form B ( · , · ) is the same as in [26, 2], while instead of the expensive full grid piecewise polynomial space V kN , the more efficient sparse grid approximation space b V kN is employed, leading to a significant reductionin the size of the linear algebraic system especially when d is large, see [25] for a detailed discussion. elow, we summarize the theoretical results for the sparse grid DG methods (2.5). Define the energynorm of a function v ∈ H (Ω N ) by ||| v ||| := X T ∈ Ω N Z T |∇ v | d x + Z Ω c | u | d x + X e ∈ Γ h Z e (cid:26) ∂v∂ n (cid:27) ds + X e ∈ Γ h Z e [ v ] ds, where H (Ω N ) denotes the standard broken Sobolev space on grid Ω N . We also need the following semi-norm. For any set L = { i , . . . , i r } ⊂ { , . . . , d } , define L c to be the complement set of L in { , . . . , d } . Fora non-negative integer α and set L , we define the mixed derivative semi-norm for a function v | v | H α,L (Ω) := (cid:13)(cid:13)(cid:13)(cid:13)(cid:18) ∂ α ∂x αi · · · ∂ α ∂x αi r (cid:19) v (cid:13)(cid:13)(cid:13)(cid:13) L (Ω) , and | v | H q +1 (Ω) := max ≤ r ≤ d (cid:18) max L ⊂{ , ,...,d } | v | H q +1 ,L (Ω) (cid:19) . The space H q +1 (Ω) denotes the closure of C ∞ (Ω) in the semi-norm | · | H q +1 .We first review several properties related to the projection operators introduced in the previous subsec-tion. We refer the reader to [25, 17] for the proofs. To avoid unnecessary clutter of constants, the notation A . B is used henceforth to represent A ≤ C · B , where the generic constant C is independent of N and themesh level considered. Lemma 2.1.
Let Q k l be the L projection operator onto the increment space W k l . For any v ∈ H p +1 (Ω) , k ≥ , ≤ q ≤ max( p, k ) , N ≥ , d ≥ , we have k Q k l ( v ) k L (Ω) . − ( q +1) | l | | v | H q +1 (Ω) . Lemma 2.2 (
Projection error estimate. ). Let b P kN be the L projection operator onto the space b V kN introduced in (2.4) . For any v ∈ H p +1 (Ω) , k ≥ , ≤ q ≤ max( p, k ) , N ≥ , d ≥ , we have (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) b P kN ( v ) − v (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . N d − Nq | v | H q +1 (Ω) . The bilinear operator B ( · , · ) is known to enjoy the following properties. Lemma 2.3 (
Orthogonality [25]. ). Let u be the exact solution to (2.1) , and u h be the numerical solutionto (2.5) , then B ( u − u h , v ) = 0 , ∀ v ∈ b V kN . Lemma 2.4 (
Boundedness and stability [2, 4]. ). When σ is taken large enough, B ( w, z ) . ||| w ||| · ||| z ||| , ∀ w, z ∈ H (Ω N ); B ( w, w ) & ||| w ||| , ∀ v ∈ b V kN . In light of the Lemmas 2.2, 2.3 and 2.4, we can prove the following error estimate for the sparse gridIPDG method (2.5).
Theorem 2.5 (
Convergence [25]. ). Let u be the exact solution to (2.1) , and u h be the numerical solutionto (2.5) . For u ∈ H p +1 (Ω) , k ≥ , ≤ q ≤ min { p, k } , N ≥ , d ≥ , we have ||| u − u h ||| . N d − Nq | u | H q +1 (Ω) . This theorem implies a convergence rate of O ( h k ) up to the polylogarithmic term | log h | d in the energynorm when u is smooth enough. efore we proceed, we provide the following estimate, which plays a crucial role in the analysis later.Also note that such an estimate is closely related to the multilevel IPDG method, see e.g. [10, 11, 9]. Lemma 2.6.
Let P kn be the L projection operator onto the full grid space V kn , n = 0 , . . . , N . Then,for any v ∈ V kN , k ≥ , N ≥ , d ≥ , we have N X n =0 h − n k ( P kn − P kn − ) v k L (Ω) . ||| v ||| , where we let P k − = 0 , and h n = 2 − n denote the mesh size of grid Ω n . The proof is given in the Appendix A. By noting that P kn − P kn − = X | l | ∞ = n Q k l , we further have, for any v ∈ b V kN , N X n =0 n X | l | ∞ = n | l | ≤ N k Q k l v k L (Ω) . ||| v ||| . (2.8)
3. Modified IPDG Method on Sparse Grids.
In this section, we formulate a novel sparse gridIPDG method for solving the Helmholtz equation with variable coefficients. The key idea is to explore thesemi-orthogonality property associated with the orthonormal multiwavelet bases.We start with the following theorem for the constant-coefficient case.
Theorem 3.1 (
Semi-orthogonality. ). If the coefficient c ≥ is constant, and two multiwavelet basisfunctions v ji , l and v j ′ i ′ , l ′ ∈ b V kN satisfy | max( l , l ′ ) | > N , then B ( v ji , l , v j ′ i ′ , l ′ ) = 0 . (3.1) Proof : Since v ji , l , v j ′ i ′ , l ′ ∈ b V kN , we have | l | ≤ N, | l ′ | ≤ N. Together with | max( l , l ′ ) | > N , we claim that there are two different dimension directions m = m sothat l m = l ′ m and l m = l ′ m . Due to the orthonormal property of the multiwavelet bases, we follow theargument for the case of the prewavelet bases in [22] and prove that (3.1) holds. (cid:4) The semi-orthogonality property actually renders a highly desired sparse structure for the resultingstiffness matrix. On the other hand, when the method (2.5) is applied to (2.1) with variable coefficients,semi-orthogonality does not hold anymore, making the stiffness matrix much denser than that in the constant-coefficient case. To recover the sparse structure for such a variable-coefficient problem, we propose to modifythe bilinear formulation in light of semi-orthogonality, motivated by the work [21, 22, 19]. In particular,based on B ( · , · ) introduced in (2.6), we define the following modified bilinear form B so ( · , · ) : b V kN × b V kN → R B so ( v ji , l , v j ′ i ′ , l ′ ) = ( B ( v ji , l , v j ′ i ′ , l ′ ) | max( l , l ′ ) | ≤ N, | max( l , l ′ ) | > N. (3.2)The sparse grid IPDG weak formulation is modified accordingly as follows. We seek u soh ∈ b V kN such that B so ( u soh , v ) = L ( v ) , ∀ v ∈ b V kN . (3.3)Note that, by construction, the resulting stiffness matrix by the modified method enjoys the same sparsestructure as the constant-coefficient case, leading to computational savings. Meanwhile, when modifying thebilinear form by means of (3.2) we in fact commit a special type of variational crimes. Below, we show that he modified method will generate a numerical solution with the same order accuracy as the original method(2.5) under an extra smoothness assumption of c .We need the following mixed derivative norm to measure the smoothness for the coefficient c . For afunction w , define the norm k w k W ∞ ,k +1 (Ω) := max ≤ α ≤ k + (cid:13)(cid:13)(cid:13)(cid:13)(cid:18) ∂ α ∂x α · · · ∂ α d ∂x α d d (cid:19) w (cid:13)(cid:13)(cid:13)(cid:13) L ∞ (Ω) , and the space W ∞ ,k +1 (Ω) := { w ∈ L ∞ (Ω) : k w k W ∞ ,k +1 (Ω) < ∞} . For convenience of illustration, wefurther introduce several shorthand notation. Whenever given two multi-indexes l and l ′ , we denote by l max = max( l , l ′ ), l min = min( l , l ′ ), and Θ l , l ′ = { s ∈ { , , . . . , d } : l s = l ′ s } .We start with the following lemma, which will play a key role in the analysis. Lemma 3.2.
Let v ji , l and v j ′ i ′ , l ′ be the two multiwavelet basis functions of space W k l and W k l ′ , respectively.Denote Q j , j ′ l , l ′ := supp ( v ji , l ) ∩ supp ( v j ′ i ′ , l ′ ) . Assume c ( x ) ∈ W ∞ ,k +1 (Ω) , then (cid:12)(cid:12)(cid:12)(cid:12)Z Ω c ( x ) v ji , l ( x ) v j ′ i ′ , l ′ ( x ) d x (cid:12)(cid:12)(cid:12)(cid:12) . − P s ∈ Θ l , l ′ l min s +( k +1)( l max s − l min s ) k c k W ∞ ,k +1 ( Q j , j ′ l , l ′ ) k v ji , l k L ( Q j , j ′ l , l ′ ) k v j ′ i ′ , l ′ k L ( Q j , j ′ l , l ′ ) . (3.4) Proof : If Q j , j ′ l , l ′ is empty, then (3.4) is trivial. Below, we assume Q j , j ′ l , l ′ is nonempty. Since the multi-dimensionalmultiwavelet bases introduced in (2.3) are constructed by tensoring the 1D bases, we can rearrange v ji , l ( x )and v j ′ i ′ , l ′ ( x ) in the integrand according to l max and l min , yielding Z Ω c ( x ) v ji , l ( x ) v j ′ i ′ , l ′ ( x ) d x = Z Q j , j ′ l , l ′ c ( x ) v j max i max , l max ( x ) v j min i min , l min ( x ) d x . (3.5)Assume Θ l , l ′ = { m , m , . . . , m r } . Note that l max s ≥ l min s + 1 iff s ∈ Θ l , l ′ . In [17], we showed that (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)Z Q j , j ′ l , l ′ c ( x ) v j max i max , l max ( x ) v j min i min , l min ( x ) d x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . − ( k +1) P s ∈ Θ l , l ′ l max s (cid:13)(cid:13)(cid:13)(cid:13)(cid:18) ∂ k +1 , ∂x k +1 m · · · ∂ k +1 ∂x k +1 m r (cid:19) (cid:16) c · v j min i min , l min (cid:17)(cid:13)(cid:13)(cid:13)(cid:13) L ( Q j , j ′ l , l ′ ) . (3.6)Together with the elementary product rule and the fact that v j min i min , l min is a piecewise polynomial of degree k in each dimension, we obtain that (cid:18) ∂ k +1 , ∂x k +1 m · · · ∂ k +1 ∂x k +1 m r (cid:19) (cid:16) c · v j min i min , l min (cid:17) = ( k + 1) r · (cid:18) ∂∂x m · · · ∂∂x m r (cid:19) c · (cid:18) ∂ k ∂x km · · · ∂ k ∂x km r (cid:19) v j min i min , l min + · · · + (cid:18) ∂ k +1 ∂x k +1 m · · · ∂ k +1 ∂x k +1 m r (cid:19) c · v j min i min , l min . A direct application of the inverse inequality leads to (cid:13)(cid:13)(cid:13)(cid:13)(cid:18) ∂ k +1 , ∂x k +1 m · · · ∂ k +1 ∂x k +1 m r (cid:19) (cid:16) c · v j min i min , l min (cid:17)(cid:13)(cid:13)(cid:13)(cid:13) L ( Q j , j ′ l , l ′ ) . k P s ∈ Θ l , l ′ l min s k c k W ∞ ,k +1 ( Q j , j ′ l , l ′ ) k v j min i min , l min k L ( Q j , j ′ l , l ′ ) . (3.7)Also note that k v j max i max , l max k L ( Q j , j ′ l , l ′ ) = 1, and by definition, k v j max i max , l max k L ( Q j , j ′ l , l ′ ) k v j min i min , l min k L ( Q j , j ′ l , l ′ ) = k v ji , l k L ( Q j , j ′ l , l ′ ) k v j ′ i ′ , l ′ k L ( Q j , j ′ l , l ′ ) . This, together with (3.5), (3.6) and (3.7), completes the proof. Lemma 3.3.
Assume c ( x ) ∈ W ∞ ,k +1 (Ω) . If w ( x ) ∈ W k l , and z ( x ) ∈ W k l ′ , then (cid:12)(cid:12)(cid:12)(cid:12)Z Ω c ( x ) w ( x ) z ( x ) d x (cid:12)(cid:12)(cid:12)(cid:12) . − P s ∈ Θ l , l ′ l min s +( k +1)( l max s − l min s ) k c k W ∞ ,k +1 (Ω) k w k L (Ω) k z k L (Ω) . Proof : Note that, for w ( x ) ∈ W k l and z ( x ) ∈ W k l ′ , we have w ( x ) = X j ∈ Θ l X ≤ i ≤ k + w ji , l v ji , l ( x ) ,z ( x ) = X j ∈ Θ l ′ X ≤ i ≤ k + z ji , l ′ v ji , l ′ ( x ) . Then, based on Lemma 3.2, the proof immediately follows the same procedure as in Lemma 4.6 in [22] andhence is omitted for brevity. (cid:4)
The above lemma directly leads to the following estimate, which is useful when proving the boundednessand stability of the modified bilinear form B so ( · , · ). Lemma 3.4.
Assume c ( x ) ∈ W ∞ ,k +1 (Ω) , k ≥ . If w ( x ) ∈ W k l , and z ( x ) ∈ W k l ′ , then (cid:12)(cid:12)(cid:12)(cid:12)Z Ω c ( x ) w ( x ) z ( x ) d x (cid:12)(cid:12)(cid:12)(cid:12) . − P s ∈ Θ l , l ′ l max s k c k W ∞ ,k +1 (Ω) k w k L (Ω) k z k L (Ω) . Proof : Since k ≥ − X s ∈ Θ l , l ′ l min s + ( k + 1)( l max s − l min s ) ≤ − X s ∈ Θ l , l ′ l min s + ( l max s − l min s ) = − X s ∈ Θ l , l ′ l max s . (cid:4) Now, we are ready to establish the boundedness and stability of the modified bilinear form B so ( · , · ). Theorem 3.5 (
Boundedness and stability with semi-orthogonality. ). Assume c ( x ) ∈ W ∞ ,k +1 (Ω) , k ≥ , d ≥ . There exists an integer N , such that B so ( w, z ) . ||| w ||| · ||| z ||| ,B so ( w, w ) & ||| w ||| , where w, z ∈ b V kN with N ≥ N .Proof : We first show that, for w, z ∈ b V kN , | B ( w, z ) − B so ( w, z ) | . − N d N d ||| w ||| · ||| z ||| , (3.8)where N d := ( N d ≤ , d − N d > . ote that | B ( w, z ) − B so ( w, z ) | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) X | l | ≤ N, | l ′ | ≤ N | l max | >N Z Ω c ( x ) Q k l ( w ) Q k l ′ ( z ) d x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ X | l | ≤ N, | l ′ | ≤ N | l max | >N (cid:12)(cid:12)(cid:12)(cid:12)Z Ω c ( x ) Q k l ( w ) Q k l ′ ( z ) d x (cid:12)(cid:12)(cid:12)(cid:12) . X | l | ≤ N, | l ′ | ≤ N | l max | >N − P s ∈ Θ l , l ′ l max s k c k W ∞ ,k +1 (Ω) k L (Ω) k Q k l ( w ) k L (Ω) k Q k l ′ ( z ) k L (Ω) = k c k W ∞ ,k +1 (Ω) X | l | ≤ N, | l ′ | ≤ N | l max | >N −| l | ∞ −| l ′ | ∞ − P s ∈ Θ l , l ′ l max s | l | ∞ k Q k l ( w ) k L (Ω) | l ′ | ∞ k Q k l ′ ( z ) k L (Ω) (3.9)where we have used Lemma 3.4. The rest of the proof follows the same procedure as in Theorem 5.4 in[22]. The only difference is that we need the Lemma 3.3 and estimate (2.8) to account for the discontinuousapproximation space. For the completeness of the paper, we choose to provide the proof.First, in [22], it was shown that if | l | ≤ N , | l ′ | ≤ N , and | l max | > N , then | l | ∞ + | l ′ | ∞ + X s ∈ Θ l , l ′ l max s ≥ | l max | − N + N d . Then, X | l | ≤ N, | l ′ | ≤ N | l max | >N −| l | ∞ −| l ′ | ∞ − P s ∈ Θ l , l ′ l max s | l | ∞ k Q k l ( w ) k L (Ω) | l ′ | ∞ k Q k l ′ ( z ) k L (Ω) ≤ − N d X | l | ≤ N, | l ′ | ≤ N | l max | >N − ( | l max | − N ) | l | ∞ k Q k l ( w ) k L (Ω) | l ′ | ∞ k Q k l ′ ( z ) k L (Ω) = 2 − N d N X γ =1 − γ X | l | ≤ N, | l ′ | ≤ N | l max | = N + γ | l | ∞ k Q k l ( w ) k L (Ω) | l ′ | ∞ k Q k l ′ ( z ) k L (Ω) ≤ − N d vuuuut N X γ =1 − γ X | l | ≤ N, | l ′ | ≤ N | l max | = N + γ | l | ∞ k Q k l ( w ) k L (Ω) vuuuut N X γ =1 − γ X | l | ≤ N, | l ′ | ≤ N | l max | = N + γ | l ′ | ∞ k Q k l ′ ( z ) k L (Ω) , (3.10)where we have used the Cauchy-Schwarz inequality. Note that N X γ =1 − γ X | l | ≤ N, | l ′ | ≤ N | l max | = N + γ | l | ∞ k Q k l ( w ) k L (Ω) = N X n =0 X | l | ≤ N | l | ∞ = n n k Q k l ( w ) k L (Ω) N X γ =1 − γ X | l ′ | ≤ N | l max | = N + γ . (3.11)We bound the right-hand side of above equation as follows. First, we need the inequality proved in [22] that,for any l with | l | ≤ N , N X γ =1 − γ X | l ′ | ≤ N | l max | = N + γ ≤ CN d , here C is a constant independent of N and l . Together with the estimate (2.8), we derive that N X n =0 X | l | ≤ N | l | ∞ = n n k Q k l ( w ) k L (Ω) N X γ =1 − γ X | l ′ | ≤ N | l max | = N + γ . N d ||| w ||| . (3.12)Combining (3.9), (3.10), (3.11) and (3.12), we prove (3.8). Due to the boundedness and stability of B ( · , · ),we complete the proof of the theorem. (cid:4) Now, we are ready to establish the main result of the paper.
Theorem 3.6 (
Convergence with semi-orthogonality. ). Let u be the exact solution to the Helmholtzequation (2.1) , and u soh be the numerical solution to the modified IPDG formulation with semi-orthogonality (3.3) . For u ∈ H p +1 (Ω) , k ≥ , c ∈ W ∞ ,k +1 (Ω) , q = max( p, k ) , d ≥ , N ≥ N , where N is given inTheorem 3.5, we have ||| u − u soh ||| . N d − qN | u | H q +1 (Ω) . Proof : Following the theory of variational crime, see [16, 8], we consider the decomposition of the error ||| u − u soh ||| ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u − b P kN ( u ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u soh − b P kN ( u ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (3.13)For the second term on the right-hand side, due to stability of B so ( · , · ) in b V kN from Theorem 3.5, we have (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u soh − b P kN u (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . sup w ∈ b V kN B so ( u soh − b P kn ( u ) , w ) ||| w ||| = sup w ∈ b V kN B so ( u soh , w ) − B so ( b P kn ( u ) , w ) ||| w ||| = sup w ∈ b V kN B ( u h , w ) − B so ( b P kn ( u ) , w ) ||| w ||| = sup w ∈ b V kN B ( u, w ) − B so ( b P kn ( u ) , w ) ||| w ||| = sup w ∈ b V kN B ( u − b P kn ( u ) , w ) + B ( b P kn ( u ) , w ) − B so ( b P kn ( u ) , w ) ||| w ||| . (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u − b P kN (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + sup w ∈ b V kN B ( b P kn ( u ) , w ) − B so ( b P kn ( u ) , w ) ||| w ||| , (3.14)where u h denotes the numerical solution by the original IPDG method (2.6). In the derivation, we havealso used orthogonality and boundedness of B ( · , · ). Note that the first term on the right-hand side is theprojection error and has been estimated in Lemma 2.2; while the second term measures the effect of modifyingthe bilinear form. For the rest of the proof, we show that | B ( b P kN ( u ) , v ) − B so ( b P kN ( u ) , v ) | . N d − ( q +1) N | u | H q +1 (Ω) ||| v ||| . (3.15) ince Q k l (cid:16) b P kN ( u ) (cid:17) = Q k l ( u ) for | l | ≤ N , we derive that | B ( b P kN ( u ) , v ) − B so ( b P kN ( u ) , v ) |≤ X | l | ≤ N, | l ′ | ≤ N | l max | >N (cid:12)(cid:12)(cid:12)(cid:12)Z Ω c ( x ) Q k l ( u ) Q k l ′ ( v ) d x (cid:12)(cid:12)(cid:12)(cid:12) . X | l | ≤ N, | l ′ | ≤ N | l max | >N − P s ∈ Θ l , l ′ l min s +( k +1)( l max s − l min s ) k c k W ∞ ,k +1 (Ω) k Q k l ( u ) k L (Ω) k Q k l ′ ( v ) k L (Ω) . k c k W ∞ ,k +1 (Ω) | u | H q +1 (Ω) X | l | ≤ N, | l ′ | ≤ N | l max | >N − P s ∈ Θ l , l ′ l min s +( k +1)( l max s − l min s ) − ( q +1) | l | k Q k l ′ ( v ) k L (Ω) , (3.16)where we have used Lemmas 3.3 and 2.1. Notice that X s ∈ Θ l , l ′ l min s + ( k + 1)( l max s − l min s ) + ( q + 1) | l | ≥ X s ∈ Θ l , l ′ ( k + 1)( l max s − l min s ) + ( q + 1) | l | ≥ ( q + 1) | l max | . Then, we have | B ( b P kN ( u ) , v ) − B so ( b P kN ( u ) , v ) | . k c k W ∞ ,k +1 (Ω) | u | H q +1 (Ω) X | l | ≤ N, | l ′ | ≤ N | l max | >N − ( q +1) | l max | k Q k l ′ ( v ) k L (Ω) = 2 − ( q +1) N k c k W ∞ ,k +1 (Ω) | u | H q +1 (Ω) N X γ =1 − ( q +1) γ X | l | ≤ N, | l ′ | ≤ N | l max | = N + γ k Q k l ′ ( v ) k L (Ω) . Similar to (3.9)-(3.12), we make use of the Cauchy-Schwarz inequality and bound the sum on the right-handside as N X γ =1 − ( q +1) γ X | l | ≤ N, | l ′ | ≤ N | l max | = N + γ k Q k l ′ ( v ) k L (Ω) ≤ N X γ =1 − γ X | l | ≤ N, | l ′ | ≤ N | l max | = N + γ k Q k l ′ ( v ) k L (Ω) . N d ||| v ||| . The proof of (3.15) is complete.By combining (3.13), (3.14), (3.15), and Lemma 2.2 about the projection error estimate, we completethe proof. (cid:4)
Remark 3.7. sup w ∈ b V kN B ( b P kN u,w ) − B so ( b P kN u,w ) ||| w ||| quantifies the variational crime from modifying the bi-linear form and is indeed one order higher than the projection error. From the numerical results presentedin the next section, we will see that if the coefficient c and the solution u are sufficiently smooth, then themodified sparse grid IPDG method will generate almost the same numerical results as the original method,while the resulting linear system is much sparser, leading to additional computational savings when solving (2.1) with variable coefficients. Remark 3.8.
The proposed framework can be extended to other variable-coefficient problems if a similarestimate as in Lemma 3.2 is derived. For instance, for Poisson’s equation with variable coefficients −∇ · ( K ( x ) ∇ u ) = f , under some condition of K , it is possible to estimate (cid:12)(cid:12)(cid:12)(cid:12)Z Ω ∇ v ji , l ( x ) · (cid:16) K ( x ) ∇ v j ′ i ′ , l ′ ( x ) (cid:17) d x (cid:12)(cid:12)(cid:12)(cid:12) and the boundary terms in the IPDG formulation with | l | ≤ N , | l | ≤ N and | l max | > N , and devise amodified sparse grid method with semi-orthogonality in an attempt to sparsify the stiffness matrix and savecost. We leave the investigation to the future work. . Numerical Results. In this section, we provide numerical results to demonstrate the performance ofthe modified sparse grid IPDG method as well as verify the analysis established for simulating the Helmholtzequations with variable coefficients. The penalty parameter σ is set to be empirical values σ = 5 · k · d . Example 4.1.
We solve the Helmholtz equation (2.1) with the smooth coefficient c ( x ) = d Y m =1 (1 − x m ) , up to d = 6 . The right-hand f is chosen such that u ( x ) = d Y m =1 sin( πx m ) . This problem was solved in [19] by a continuous Ritz-Galerkin discretization on sparse grids usingprewavelets in conjunction with semi-orthogonality. We first consider the case of d = 2 and compare andcontrast the performance of both modified and original sparse grid IPDG methods in terms of accuracyand efficiency. In Table 4.1, we summarize the convergence study for both methods with k = 1 , ,
3. Itis observed that, for this problem, the L and H errors given by the two methods are virtually the same,and the associated convergence rates are close to k + 1-th order for the L error and k -th order for the H error. Such an observation validates the convergence analysis in Theorem 3.6 that the variational crimeincurred by the proposed modification of the bilinear form is one order higher than the projection error andhence becoming negligible in the simulation when the smoothness requirements are fulfilled. To demonstratethe efficiency advantage of the modified method, we provide the numbers of nonzero entries of the stiffnessmatrices as well as the conditional numbers for both methods with various mesh levels N and polynomialdegrees k in Table 4.2. Further, in Figure 4.1, we plot the sparsity patterns of the matrices with N = 5 and k = 1. One can see that the modified method enjoys much sparser stiffness matrices than the original one,while the condition numbers remain almost unchanged. More specifically, we compute the order of sparsitydefined by O s = log(NNZ) / log(DOF), where NNZ is referred to as the number of nonzero elements. Notethat NNZ scales as O (DOF . ) for the modified method; while it scales as O (DOF . ) for the original method.Such reduction leads to great computational savings when assembling as well as solving the algebraic linearsystem. nz = 33160 (a) (b) Fig. 4.1 . Example 4.1. The sparsity patterns of the stiffness matrices by (a) the modified sparse grid IPDG method and(b) the original sparse grid IPDG method. N = 5 , k = 1 , d = 2 . We then consider the case of d = 3 and perform a similar comparison between the two IPDG methods.We provide the convergence study results in Table 4.3 with k = 1 , , , including the L and H errorsand the associated orders of accuracy for both methods. As in the above d = 2 case, both methods give able 4.1 Numerical errors and orders of accuracy for Example 4.1 computed by the modified and the original IPDG methods with k = 1 , , . d = 2 . Modified IPDG method Original IPDG method
N L error order H error order L error order H error order k = 12 6.12E-02 6.88E-01 6.12E-02 6.88E-013 1.62E-02 1.92 3.37E-01 1.03 1.62E-02 1.92 3.37E-01 1.034 4.02E-03 2.01 1.66E-01 1.02 4.02E-03 2.01 1.66E-01 1.025 9.99E-04 2.01 8.26E-02 1.01 9.99E-04 2.01 8.26E-02 1.016 2.52E-04 1.99 4.11E-02 1.01 2.52E-04 1.99 4.11E-02 1.01 k = 22 1.52E-03 5.27E-02 1.52E-03 5.27E-023 2.27E-04 2.74 1.34E-02 1.98 2.27E-04 2.74 1.34E-02 1.984 3.51E-05 2.70 3.35E-03 2.00 3.51E-05 2.70 3.35E-03 2.005 5.27E-06 2.74 8.37E-04 2.00 5.27E-06 2.74 8.37E-04 2.006 7.61E-07 2.79 2.09E-04 2.00 7.61E-07 2.79 2.09E-04 2.00 k = 32 9.17E-05 3.53E-03 9.17E-05 3.53E-033 6.03E-06 3.93 4.37E-04 3.01 6.03E-06 3.93 4.37E-04 3.014 3.84E-07 3.97 5.43E-05 3.01 3.84E-07 3.97 5.43E-05 3.015 2.43E-08 3.98 6.78E-06 3.00 2.43E-08 3.98 6.78E-06 3.006 1.54E-09 3.98 8.47E-07 3.00 1.54E-09 3.98 8.47E-07 3.00 Table 4.2
Sparsity and condition number of the stiffness matrices for the modified and original IPDG methods. k = 1 , , . d = 2 .DOF is the number of degrees of freedom used for the sparse grid IPDG methods. NNZ is the number of nonzero elements inthe stiffness matrix. O s =log(NNZ)/ log(DOF). κ is the condition number. d = 2 Modified IPDG method Original IPDG method N DOF NNZ O s κ NNZ O s κ k = 12 32 592 1.84 8.23E+01 976 1.99 8.23E+013 80 2608 1.80 3.49E+02 5424 1.96 3.49E+024 448 9808 1.51 1.40E+03 26192 1.67 1.40E+035 1024 33160 1.50 5.54E+03 116122 1.68 5.54E+036 2304 103968 1.49 2.20E+04 493154 1.69 2.20E+04 k = 22 180 2976 1.54 3.51E+02 4920 1.64 3.51E+023 432 12864 1.56 1.37E+03 27120 1.68 1.37E+034 1008 48132 1.56 5.35E+03 131076 1.70 5.35E+035 2304 162324 1.55 2.11E+04 580352 1.71 2.11E+046 5184 509616 1.54 8.37E+04 2460726 1.72 8.37E+04 k = 32 128 9216 1.88 8.53E+02 15616 1.99 8.53E+023 768 39952 1.59 3.29E+03 85760 1.71 3.29E+034 1792 149264 1.59 1.28E+04 413952 1.73 1.28E+045 4096 505072 1.58 5.04E+04 1848064 1.73 5.04E+046 9216 1587200 1.56 2.00E+05 7869696 1.74 2.00E+05 able 4.3 Numerical errors and orders of accuracy for Example 4.1 computed by the modified and the original IPDG methods with k = 1 , , . d = 3 . Modified IPDG method Original IPDG method
N L error order H error order L error order H error order k = 12 1.30E-01 9.45E-01 1.30E-01 9.45E-013 3.54E-02 1.87 4.41E-01 1.10 3.54E-02 1.87 4.41E-01 1.104 8.88E-03 1.99 2.07E-01 1.09 8.88E-03 1.99 2.07E-01 1.095 2.16E-03 2.04 9.90E-02 1.07 2.16E-03 2.04 9.90E-02 1.076 5.42E-04 1.99 4.82E-02 1.04 5.42E-04 1.99 4.82E-02 1.04 k = 22 1.13E-03 4.54E-02 1.13E-03 4.54E-023 2.10E-04 2.43 1.19E-02 1.94 2.10E-04 2.43 1.19E-02 1.944 3.73E-05 2.49 3.01E-03 1.98 3.73E-05 2.49 3.01E-03 1.985 6.15E-06 2.60 7.54E-04 2.00 6.15E-06 2.60 7.54E-04 2.006 9.71E-07 2.66 1.88E-04 2.00 9.71E-07 2.66 1.88E-04 2.00 k = 32 1.12E-04 1.30E-03 1.12E-04 1.30E-033 7.40E-06 3.92 1.39E-04 3.23 7.40E-06 3.92 1.39E-04 3.234 4.72E-07 3.97 1.58E-05 3.14 4.72E-07 3.97 1.58E-05 3.145 2.99E-08 3.98 1.86E-06 3.08 2.99E-08 3.98 1.86E-06 3.086 1.89E-09 3.98 2.26E-07 3.05 1.89E-09 3.98 2.26E-07 3.05almost the same numerical errors. k + 1-th and k -th order of accuracy is observed for the L and H errors,respectively, which agrees with the error estimates established in Theorem 3.6. In Figure 4.2, we plot thesparsity patterns of the stiffness matrices by two methods with N = 4 and k = 1. As expected, the stiffnessmatrix by the modified method is much sparser than that by the original method. Furthermore, we provideDOF, NNZ and the associated order of sparsity O s for both methods in Table 4.4. It is observed that NNZscales as O (DOF . ) for the modified method; while it scales as O (DOF . ) for the original method. (a) (b) Fig. 4.2 . Example 4.1. The sparsity patterns of the stiffness matrices by (a) the modified sparse grid IPDG method and(b) the original sparse grid IPDG method. N = 4 , k = 1 , d = 3 . Last, we summarize the simulation results for d = 4 , , L errors, the associated orders of accuracy, NNZ, and the orders of sparsity O s for themodified method only. The observation is similar to that in the cases of d = 2 ,
3. The modified method able 4.4 Sparsity of the stiffness matrices for the modified and original IPDG methods. k = 1 , , . d = 3 . DOF is the numberof degrees of freedom used for the sparse grid IPDG methods. NNZ is the number of nonzero elements in the stiffness matrix. O s =log(NNZ)/ log(DOF). d = 3 Modified IPDG method Original IPDG method N DOF NNZ O s NNZ O s k = 12 104 4.31E+03 1.80 1.05E+04 1.993 304 2.34E+04 1.76 8.26E+04 1.984 832 1.08E+05 1.72 5.51E+05 1.975 2176 4.47E+05 1.69 3.27E+06 1.956 5504 1.69E+06 1.67 1.80E+07 1.94 k = 22 351 4.93E+04 1.84 1.19E+05 1.993 1026 2.65E+05 1.80 9.38E+05 1.984 2808 1.22E+06 1.76 6.26E+06 1.975 7344 5.03E+06 1.73 3.73E+07 1.966 18576 1.91E+07 1.71 2.06E+08 1.95 k = 32 832 2.76E+05 1.86 6.69E+05 1.993 2432 1.48E+06 1.82 5.26E+06 1.994 6656 6.81E+06 1.79 3.51E+07 1.975 17408 2.81E+07 1.76 2.10E+08 1.966 44032 1.07E+08 1.73 1.16E+09 1.95is able to achieve the expected orders of convergence, while slight order reduction is observed due to thepolylogarithmic term appearing in the error estimate. Note that unlike the modified method, the originalmethod suffers the almost fully populated stiffness matrix, making its simulation exceed our computationresource limit. Currently, we are allowed to handle matrices with NNZ no greater than 1 E + 10. Based onthe comparison drawn in d = 2 ,
5. Conclusion.
In this paper, we developed a modified sparse grid IPDG method using semi-orthogonalityfor solving the Helmholtz equation with variable coefficients. The original IPDG features a sparse finite el-ement space which scales as O ( h − | log h | d − ) for d -dimensional problems, translating into a significantcost reduction when d is large. On the other hand, when applied to the variable-coefficient problem, themethod suffers a dense stiffness matrix, impeding its efficiency advantage to some extent. Based on thesemi-orthogonality property associated with the orthonormal multiwavelet bases, the IPDG bilinear form ismodified aiming to sparsify the resulting linear algebraic matrix. A numerical analysis demonstrates that theerror incurred by the modification is one order higher than the projection error, making it negligible in thesimulations. Numerical results in up to six dimensions are shown to validate the analysis and demonstratethat the modified method enjoys sparser stiffness matrix than the original one, leading to extra computa-tional savings. Future work includes the study of adaptive algorithm for less regular solutions as well asextension of the method to other types of high-dimensional variable-coefficient equations. A. Proof of Lemma 2.6.
In the proof, we need the averaging projection operator A proposed in[10, 11], which decomposes the space V kN into the conforming subspace V cN and nonconforming part V ncN ,i.e, V cN := A V kN ⊂ V kN ∩ H (Ω) , V ncN := V kN − V cN . Furthermore, A satisfies the Jackson estimate |||A v ||| a . ||| v ||| a , ||| ( I − A ) v ||| a . h ||| v ||| a , v ∈ V kN , (A.1) able 4.5 The L errors, orders of accuracy and sparsity of the stiffness matrices for the modified IPDG methods. k = 1 , , . d = 4 . DOF is the number of degrees of freedom used for the sparse grid IPDG methods. NNZ is the number of nonzeroelements in the stiffness matrix. O s =log(NNZ)/ log(DOF). d = 4 Modified IPDG method N DOF L error order NNZ O s k = 12 304 1.49E-01 2.76E+04 1.793 1008 7.35E-02 1.02 1.78E+05 1.754 3072 2.15E-02 1.77 9.71E+05 1.725 8832 5.62E-03 1.93 4.69E+06 1.696 24320 1.44E-03 1.97 2.06E+07 1.67 k = 22 1539 1.28E-03 7.04E+05 1.833 5103 2.25E-04 2.51 4.52E+06 1.804 15552 3.94E-05 2.51 2.47E+07 1.765 44712 6.63E-06 2.57 1.19E+08 1.746 123120 1.09E-06 2.61 5.25E+08 1.71 k = 32 4864 7.83E-05 7.02E+06 1.863 16128 5.30E-06 3.88 4.51E+07 1.824 49152 3.44E-07 3.95 2.46E+08 1.79 Table 4.6
The L errors, orders of accuracy and sparsity of the stiffness matrices for the modified IPDG methods. k = 1 , , . d = 5 . DOF is the number of degrees of freedom used for the sparse grid IPDG methods. NNZ is the number of nonzeroelements in the stiffness matrix. O s =log(NNZ)/ log(DOF). d = 5 Modified IPDG method N DOF L error order NNZ O s k = 12 832 1.30E-01 1.60E+05 1.783 3072 9.37E-02 0.47 1.20E+06 1.744 10272 4.10E-02 1.19 7.52E+06 1.715 32064 1.31E-02 1.65 4.15E+07 1.696 95104 3.71E-03 1.81 2.07E+08 1.67 k = 22 6318 1.06E-03 9.22E+06 1.833 23328 1.97E-04 2.42 6.89E+07 1.794 78003 3.60E-05 2.45 4.32E+08 1.775 243486 6.29E-06 2.52 2.38E+09 1.746 722196 1.08E-06 2.54 1.19E+10 1.72 k = 32 26624 6.79E-05 1.64E+08 1.863 98304 4.69E-06 3.85 1.22E+09 1.824 328704 3.10E-07 3.92 7.67E+09 1.79where the energy norm ||| v ||| a is defined by ||| v ||| a := X T ∈ Ω N Z T |∇ v | d x + X e ∈ Γ h Z e [ v ] ds, v ∈ V kN , able 4.7 The L errors, orders of accuracy and sparsity of the stiffness matrices for the modified IPDG methods. k = 1 , . d = 6 .DOF is the number of degrees of freedom used for the sparse grid IPDG methods. NNZ is the number of nonzero elements inthe stiffness matrix. O s =log(NNZ)/ log(DOF). d = 6 Modified IPDG method N DOF L error order NNZ O s k = 12 2176 1.03E-01 8.78E+05 1.783 8832 8.71E-02 0.24 7.48E+06 1.744 32064 5.70E-02 0.61 5.29E+07 1.715 107712 2.33E-02 1.29 3.27E+08 1.696 341504 7.85E-03 1.57 1.81E+09 1.67 k = 22 24786 8.47E-04 1.14E+08 1.833 100602 1.67E-04 2.35 9.69E+08 1.804 365229 3.13E-05 2.41 6.85E+09 1.77and I denotes the identity operator. The construction and detailed analysis of A were established in[10, 11, 9]. Clearly, ||| v ||| a ≤ ||| v ||| .For any v ∈ V kN , we denote by w := A ( v ) and z := ( I − A ) v . We first prove that N X n =0 h − n k ( P kn − P kn − ) w k L (Ω) . | w | H (Ω) . (A.2)Note that there is a standard nested relation for conforming finite element spaces: V c ⊂ V c ⊂ · · · V cN ⊂ H (Ω) , where the space V cn defined on mesh Ω n is a subspace of V kn , n = 0 , . . . , N . Following the idea of the BPXmultilevel precondtioner [7], we define the Ritz projection operator R n : V cN → V cn a ( R n ( w ) , r ) = a ( w, r ) ∀ r ∈ V cn , n = 1 , . . . , N, where a ( w, r ) = Z Ω ∇ w · ∇ rd x , and we let R − = 0. Denote w n := ( R n − R n − ) w ∈ V cn . Then, w = P Nn =0 w n . Since V cn ⊂ V kn ,( P kl − P kl − ) w n = 0 , l > n. (A.3)In addition, the Ritz operator enjoys the property R n − ( R n − R n − ) = 0 . Together with the standard estimate of R n − by the conforming finite element analysis, we have k w n k L (Ω) = k ( I − R n − ) w n k L (Ω) . h n | w n | H (Ω) . (A.4) hen, we derive N X l =0 h − l k ( P kl − P kl − ) w k L (Ω) = N X l =0 h − l k ( P kl − P kl − ) N X n =0 w n k L (Ω) = N X l =0 h − l k N X n =0 ( P kl − P kl − ) w n k L (Ω) = N X l =0 h − l N X n,m =0 (( P kl − P kl − ) w n , ( P kl − P kl − ) w m )= N X l =0 h − l N X n,m = l (( P kl − P kl − ) w n , ( P kl − P kl − ) w m ) ( due to (A.3))= N X n,m =0 min( n,m ) X l =0 h − l (( P kl − P kl − ) w n , ( P kl − P kl − ) w m ) ≤ N X n,m =0 min( n,m ) X l =0 h − l k ( P kl − P kl − ) w n k L (Ω) k ( P kl − P kl − ) w m k L (Ω) . N X n,m =0 min( n,m ) X l =0 h − l k w n k L (Ω) k w m k L (Ω) . N X n,m =0 min( n,m ) X l =0 h − l h n h m | w n | H (Ω) | w m | H (Ω) ( due to (A.4)) . N X n,m =0 (cid:18) (cid:19) | n − m | | w n | H (Ω) | w m | H (Ω) . By an elementary linear algebra result, see [8], together with the identity a ( w n , w m ) = 0, n = m , we finallyget N X n,m =0 (cid:18) (cid:19) | n − m | | w n | H (Ω) | w m | H (Ω) . N X n =0 | w n | H (Ω) = | w | H (Ω) , and (A.2) is proved. Together with (A.2) and the Jackson estimate (A.1), we immediately obtain N X n =0 h − n k ( P kn − P kn − ) w k L (Ω) . ||| w ||| a . ||| v ||| a . (A.5)Now we derive the estimate for z . A direct application of the Poincar´e-Friedrichs type inequality forspace V kN , see [2], yields k z k L (Ω) . ||| z ||| a . This, together with the Jackson estimate (A.1), gives N X n =0 h − n k ( P kn − P kn − ) z k L (Ω) . N X n =0 h − n k z k L (Ω) . N X n =0 h − n ||| z ||| a . h N X n =0 h − n ! ||| v ||| a . ||| v ||| a . (A.6)Combining the estimates (A.5) and (A.6), we complete the proof as N X n =0 h − n k ( P kn − P kn − ) v k L (Ω) . N X n =0 h − n k ( P kn − P kn − ) w k L (Ω) + N X n =0 h − n k ( P kn − P kn − ) z k L (Ω) . ||| v ||| a ≤ ||| v ||| . (cid:4) SIAM J. Math. Anal. , 24(1):246–262,1993.[2] D. Arnold. An interior penalty finite element method with discontinuous elements.
SIAM J. Numer. Anal. , 19(4):742–760,1982.[3] D. Arnold, F. Brezzi, B. Cockburn, and L. Marini. Unified analysis of discontinuous Galerkin methods for elliptic problems.
SIAM J. Numer. Anal. , 39(5):1749–1779, 2002.[4] D. Arnold, F. Brezzi, B. Cockburn, and L. Marini. Unified analysis of discontinuous Galerkin methods for elliptic problems.
SIAM J. Numer. Anal. , 39:1749–1779, 2002.[5] G. Baker. Finite element methods for elliptic equations using nonconforming elements.
Math. Comp. , 31(137):45–59,1977.[6] R. Bellman.
Adaptive control processes: a guided tour , volume 4. Princeton University Press Princeton, 1961.[7] J. H. Bramble, J. E. Pasciak, and J. Xu. Parallel multilevel preconditioners.
Math. Comput. , 55(191):1–22, 1990.[8] S. Brenner and R. Scott.
The mathematical theory of finite element methods , volume 15. Springer Science & BusinessMedia, 2007.[9] K. Brix, M. Campos Pinto, C. Canuto, and W. Dahmen. Multilevel preconditioning of discontinuous Galerkin spectralelement methods. Part I: geometrically conforming meshes.
IMA J. Numer. Anal. , 35(4):1487–1532, 2014.[10] K. Brix, M. C. Pinto, and W. Dahmen. A multilevel preconditioner for the interior penalty discontinuous Galerkin method.
SIAM J. Numer. Anal. , 46(5):2742–2768, 2008.[11] K. Brix, M. C. Pinto, W. Dahmen, and R. Massjung. Multilevel preconditioners for the interior penalty discontinuousGalerkin method II - quantitative studies.
Commun. Comput. Phys. , 5(2-4):296–325, 2009.[12] H.-J. Bungartz and M. Griebel. Sparse grids.
Acta numer. , 13:147–269, 2004.[13] B. Cockburn, G. E. Karniadakis, and C.-W. Shu. The development of discontinuous Galerkin methods. In
DiscontinuousGalerkin Methods , pages 3–50. Springer, 2000.[14] J. Douglas and T. Dupont. Interior penalty procedures for elliptic and parabolic galerkin methods. In
Computing methodsin applied sciences , pages 207–216. Springer, 1976.[15] J. Garcke and M. Griebel.
Sparse grids and applications . Springer, 2013.[16] C. Großmann and H.-G. Roos.
Numerik partieller Differentialgleichungen . Springer-Verlag, 1992.[17] W. Guo and Y. Cheng. A sparse grid discontinuous Galerkin method for high-dimensional transport equations and itsapplication to kinetic simulations.
SIAM J. Sci. Comput. , 38(6):A3381–A3409, 2016.[18] W. Guo and Y. Cheng. An adaptive multiresolution discontinuous Galerkin method for time-dependent transport equationsin multidimensions.
SIAM J. Sci. Comput. , 39(6):A2962–A2992, 2017.[19] R. Hartmann and C. Pflaum. A prewavelet-based algorithm for the solution of second-order elliptic differential equationswith variable coefficients on sparse grids.
Numer. Algorithms , pages 1–28, 2018.[20] Y. Liu, Y. Cheng, S. Chen, and Y.-T. Zhang. Krylov implicit integration factor discontinuous galerkin methods on sparsegrids for high dimensional reaction-diffusion equations.
J. Comput. Phys. , 388:90 – 102, 2019.[21] C. Pflaum. A multilevel algorithm for the solution of second order elliptic differential equations on sparse grids.
Numer.Math. , 79(1):141–155, 1998.[22] C. Pflaum and R. Hartmann. A sparse grid discretization of the helmholtz equation with variable coefficients in highdimensions.
SIAM J. Numer. Anal. , 54(4):2707–2727, 2016.[23] Z. Tao, A. Chen, M. Zhang, and Y. Cheng. Sparse grid central discontinuous Galerkin method for linear hyperbolicsystems in high dimensions. arXiv preprint arXiv:1807.04137 , 2018.[24] Z. Tao, W. Guo, and Y. Cheng. Sparse grid discontinuous Galerkin methods for the Vlasov-Maxwell system.
J. Comput.Phys. , accepted, 2019.[25] Z. Wang, Q. Tang, W. Guo, and Y. Cheng. Sparse grid discontinuous Galerkin methods for high-dimensional ellipticequations.
J. Comput. Phys. , 314:244–263, 2016.[26] M. Wheeler. An elliptic collocation-finite element method with interior penalties.
SIAM J. Numer. Anal. , 15(1):152–161,1978.[27] C. Zenger. Sparse grids. In