A sparse grid discontinuous Galerkin method for the high-dimensional Helmholtz equation with variable coefficients

TL;DR

This paper develops an efficient sparse grid discontinuous Galerkin method for high-dimensional Helmholtz equations with variable coefficients, achieving significant computational savings while maintaining accuracy.

Contribution

It introduces a semi-orthogonality property to enhance the sparse grid DG method's efficiency for high-dimensional Helmholtz problems with variable coefficients.

Findings

Achieves same order accuracy as original method

Produces a much sparser stiffness matrix

Verifies effectiveness with numerical tests up to six dimensions

Abstract

The simulation of high-dimensional problems with manageable computational resource represents a long-standing challenge. In a series of our recent work [25, 17, 18, 24], a class of sparse grid DG methods has been formulated for solving various types of partial differential equations in high dimensions. By making use of the multiwavelet tensor-product bases on sparse grids in conjunction with the standard DG weak formulation, such a novel method is able to significantly reduce the computation and storage cost compared with full grid DG counterpart, while not compromising accuracy much for sufficiently smooth solutions. In this paper, we consider the high-dimensional Helmholtz equation with variable coefficients and demonstrate that for such a problem the efficiency of the sparse grid DG method can be further enhanced by exploring a semi-orthogonality property associated with the multiwavelet bases, motivated by the work [21, 22, 19]. The detailed convergence analysis shows that the modified sparse grid DG method attains the same order accuracy, but the resulting stiffness matrix is much sparser than that by the original method, leading to extra computational savings. Numerical tests up to six dimensions are provided to verify the analysis.

Tables

Table 1. Table 1: Numerical errors and orders of accuracy for Example 4.1 computed by the modified and the original IPDG methods with k = 1 ,  2 ,  3 𝑘 1 2 3 k=1,\,2,\,3 . d = 2 𝑑 2 d=2 .

	Modified IPDG method				Original IPDG method
$N$	$L^2$ error	order	$H^1$ error	order	$L^2$ error	order	$H^1$ error	order
	$k=1$
2	6.12E-02		6.88E-01		6.12E-02		6.88E-01
3	1.62E-02	1.92	3.37E-01	1.03	1.62E-02	1.92	3.37E-01	1.03
4	4.02E-03	2.01	1.66E-01	1.02	4.02E-03	2.01	1.66E-01	1.02
5	9.99E-04	2.01	8.26E-02	1.01	9.99E-04	2.01	8.26E-02	1.01
6	2.52E-04	1.99	4.11E-02	1.01	2.52E-04	1.99	4.11E-02	1.01
	$k=2$
2	1.52E-03		5.27E-02		1.52E-03		5.27E-02
3	2.27E-04	2.74	1.34E-02	1.98	2.27E-04	2.74	1.34E-02	1.98
4	3.51E-05	2.70	3.35E-03	2.00	3.51E-05	2.70	3.35E-03	2.00
5	5.27E-06	2.74	8.37E-04	2.00	5.27E-06	2.74	8.37E-04	2.00
6	7.61E-07	2.79	2.09E-04	2.00	7.61E-07	2.79	2.09E-04	2.00
	$k=3$
2	9.17E-05		3.53E-03		9.17E-05		3.53E-03
3	6.03E-06	3.93	4.37E-04	3.01	6.03E-06	3.93	4.37E-04	3.01
4	3.84E-07	3.97	5.43E-05	3.01	3.84E-07	3.97	5.43E-05	3.01
5	2.43E-08	3.98	6.78E-06	3.00	2.43E-08	3.98	6.78E-06	3.00
6	1.54E-09	3.98	8.47E-07	3.00	1.54E-09	3.98	8.47E-07	3.00

Table 2. Table 2: Sparsity and condition number of the stiffness matrices for the modified and original IPDG methods. k = 1 ,  2 ,  3 𝑘 1 2 3 k=1,\,2,\,3 . d = 2 𝑑 2 d=2 . DOF is the number of degrees of freedom used for the sparse grid IPDG methods. NNZ is the number of nonzero elements in the stiffness matrix. O s subscript 𝑂 𝑠 O_{s} =log(NNZ)/ log(DOF). κ 2 subscript 𝜅 2 \kappa_{2} is the condition number.

$d=2$		Modified IPDG method			Original IPDG method
$N$	DOF	NNZ	$O_s$	${\kappa }_2$	NNZ	$O_s$	${\kappa }_2$
		$k=1$
2	32	592	1.84	8.23E+01	976	1.99	8.23E+01
3	80	2608	1.80	3.49E+02	5424	1.96	3.49E+02
4	448	9808	1.51	1.40E+03	26192	1.67	1.40E+03
5	1024	33160	1.50	5.54E+03	116122	1.68	5.54E+03
6	2304	103968	1.49	2.20E+04	493154	1.69	2.20E+04
		$k=2$
2	180	2976	1.54	3.51E+02	4920	1.64	3.51E+02
3	432	12864	1.56	1.37E+03	27120	1.68	1.37E+03
4	1008	48132	1.56	5.35E+03	131076	1.70	5.35E+03
5	2304	162324	1.55	2.11E+04	580352	1.71	2.11E+04
6	5184	509616	1.54	8.37E+04	2460726	1.72	8.37E+04
		$k=3$
2	128	9216	1.88	8.53E+02	15616	1.99	8.53E+02
3	768	39952	1.59	3.29E+03	85760	1.71	3.29E+03
4	1792	149264	1.59	1.28E+04	413952	1.73	1.28E+04
5	4096	505072	1.58	5.04E+04	1848064	1.73	5.04E+04
6	9216	1587200	1.56	2.00E+05	7869696	1.74	2.00E+05

Table 3. Table 3: Numerical errors and orders of accuracy for Example 4.1 computed by the modified and the original IPDG methods with k = 1 ,  2 ,  3 𝑘 1 2 3 k=1,\,2,\,3 . d = 3 𝑑 3 d=3 .

	Modified IPDG method				Original IPDG method
$N$	$L^2$ error	order	$H^1$ error	order	$L^2$ error	order	$H^1$ error	order
	$k=1$
2	1.30E-01		9.45E-01		1.30E-01		9.45E-01
3	3.54E-02	1.87	4.41E-01	1.10	3.54E-02	1.87	4.41E-01	1.10
4	8.88E-03	1.99	2.07E-01	1.09	8.88E-03	1.99	2.07E-01	1.09
5	2.16E-03	2.04	9.90E-02	1.07	2.16E-03	2.04	9.90E-02	1.07
6	5.42E-04	1.99	4.82E-02	1.04	5.42E-04	1.99	4.82E-02	1.04
	$k=2$
2	1.13E-03		4.54E-02		1.13E-03		4.54E-02
3	2.10E-04	2.43	1.19E-02	1.94	2.10E-04	2.43	1.19E-02	1.94
4	3.73E-05	2.49	3.01E-03	1.98	3.73E-05	2.49	3.01E-03	1.98
5	6.15E-06	2.60	7.54E-04	2.00	6.15E-06	2.60	7.54E-04	2.00
6	9.71E-07	2.66	1.88E-04	2.00	9.71E-07	2.66	1.88E-04	2.00
	$k=3$
2	1.12E-04		1.30E-03		1.12E-04		1.30E-03
3	7.40E-06	3.92	1.39E-04	3.23	7.40E-06	3.92	1.39E-04	3.23
4	4.72E-07	3.97	1.58E-05	3.14	4.72E-07	3.97	1.58E-05	3.14
5	2.99E-08	3.98	1.86E-06	3.08	2.99E-08	3.98	1.86E-06	3.08
6	1.89E-09	3.98	2.26E-07	3.05	1.89E-09	3.98	2.26E-07	3.05

Table 4. Table 4: Sparsity of the stiffness matrices for the modified and original IPDG methods. k = 1 ,  2 ,  3 𝑘 1 2 3 k=1,\,2,\,3 . d = 3 𝑑 3 d=3 . DOF is the number of degrees of freedom used for the sparse grid IPDG methods. NNZ is the number of nonzero elements in the stiffness matrix. O s subscript 𝑂 𝑠 O_{s} =log(NNZ)/ log(DOF).

$d=3$		Modified IPDG method		Original IPDG method
$N$	DOF	NNZ	$O_s$	NNZ	$O_s$
		$k=1$
2	104	4.31E+03	1.80	1.05E+04	1.99
3	304	2.34E+04	1.76	8.26E+04	1.98
4	832	1.08E+05	1.72	5.51E+05	1.97
5	2176	4.47E+05	1.69	3.27E+06	1.95
6	5504	1.69E+06	1.67	1.80E+07	1.94
		$k=2$
2	351	4.93E+04	1.84	1.19E+05	1.99
3	1026	2.65E+05	1.80	9.38E+05	1.98
4	2808	1.22E+06	1.76	6.26E+06	1.97
5	7344	5.03E+06	1.73	3.73E+07	1.96
6	18576	1.91E+07	1.71	2.06E+08	1.95
		$k=3$
2	832	2.76E+05	1.86	6.69E+05	1.99
3	2432	1.48E+06	1.82	5.26E+06	1.99
4	6656	6.81E+06	1.79	3.51E+07	1.97
5	17408	2.81E+07	1.76	2.10E+08	1.96
6	44032	1.07E+08	1.73	1.16E+09	1.95

Table 5. Table 5: The L 2 superscript 𝐿 2 L^{2} errors, orders of accuracy and sparsity of the stiffness matrices for the modified IPDG methods. k = 1 ,  2 ,  3 𝑘 1 2 3 k=1,\,2,\,3 . d = 4 𝑑 4 d=4 . DOF is the number of degrees of freedom used for the sparse grid IPDG methods. NNZ is the number of nonzero elements in the stiffness matrix. O s subscript 𝑂 𝑠 O_{s} =log(NNZ)/ log(DOF).

$d=4$		Modified IPDG method
$N$	DOF	$L^2$ error	order	NNZ	$O_s$
$k=1$
2	304	1.49E-01		2.76E+04	1.79
3	1008	7.35E-02	1.02	1.78E+05	1.75
4	3072	2.15E-02	1.77	9.71E+05	1.72
5	8832	5.62E-03	1.93	4.69E+06	1.69
6	24320	1.44E-03	1.97	2.06E+07	1.67
$k=2$
2	1539	1.28E-03		7.04E+05	1.83
3	5103	2.25E-04	2.51	4.52E+06	1.80
4	15552	3.94E-05	2.51	2.47E+07	1.76
5	44712	6.63E-06	2.57	1.19E+08	1.74
6	123120	1.09E-06	2.61	5.25E+08	1.71
$k=3$
2	4864	7.83E-05		7.02E+06	1.86
3	16128	5.30E-06	3.88	4.51E+07	1.82
4	49152	3.44E-07	3.95	2.46E+08	1.79

Table 6. Table 6: The L 2 superscript 𝐿 2 L^{2} errors, orders of accuracy and sparsity of the stiffness matrices for the modified IPDG methods. k = 1 ,  2 ,  3 𝑘 1 2 3 k=1,\,2,\,3 . d = 5 𝑑 5 d=5 . DOF is the number of degrees of freedom used for the sparse grid IPDG methods. NNZ is the number of nonzero elements in the stiffness matrix. O s subscript 𝑂 𝑠 O_{s} =log(NNZ)/ log(DOF).

$d=5$		Modified IPDG method
$N$	DOF	$L^2$ error	order	NNZ	$O_s$
$k=1$
2	832	1.30E-01		1.60E+05	1.78
3	3072	9.37E-02	0.47	1.20E+06	1.74
4	10272	4.10E-02	1.19	7.52E+06	1.71
5	32064	1.31E-02	1.65	4.15E+07	1.69
6	95104	3.71E-03	1.81	2.07E+08	1.67
$k=2$
2	6318	1.06E-03		9.22E+06	1.83
3	23328	1.97E-04	2.42	6.89E+07	1.79
4	78003	3.60E-05	2.45	4.32E+08	1.77
5	243486	6.29E-06	2.52	2.38E+09	1.74
6	722196	1.08E-06	2.54	1.19E+10	1.72
$k=3$
2	26624	6.79E-05		1.64E+08	1.86
3	98304	4.69E-06	3.85	1.22E+09	1.82
4	328704	3.10E-07	3.92	7.67E+09	1.79

Table 7. Table 7: The L 2 superscript 𝐿 2 L^{2} errors, orders of accuracy and sparsity of the stiffness matrices for the modified IPDG methods. k = 1 ,  2 𝑘 1 2 k=1,\,2 . d = 6 𝑑 6 d=6 . DOF is the number of degrees of freedom used for the sparse grid IPDG methods. NNZ is the number of nonzero elements in the stiffness matrix. O s subscript 𝑂 𝑠 O_{s} =log(NNZ)/ log(DOF).

$d=6$		Modified IPDG method
$N$	DOF	$L^2$ error	order	NNZ	$O_s$
$k=1$
2	2176	1.03E-01		8.78E+05	1.78
3	8832	8.71E-02	0.24	7.48E+06	1.74
4	32064	5.70E-02	0.61	5.29E+07	1.71
5	107712	2.33E-02	1.29	3.27E+08	1.69
6	341504	7.85E-03	1.57	1.81E+09	1.67
$k=2$
2	24786	8.47E-04		1.14E+08	1.83
3	100602	1.67E-04	2.35	9.69E+08	1.80
4	365229	3.13E-05	2.41	6.85E+09	1.77

Equations

$$-\Delta u+cu=$$

$$u=$$

$$V_l^k=\{v:v\in P^k(I_l^j),\forall j=0,\dots ,2^l-1\}$$

$$V_0^k\subset V_1^k\subset V_2^k\subset V_3^k\subset \cdots ,$$

$$V_{l-1}^k\oplus W_l^k=V_l^k,W_l^k\perp V_{l-1}^k.$$

$$V_N^k=\underset{0\le l\le N}{⨁}{W}_l^k.$$

$$v_{i,l}^j(x),i=1,\dots ,k+1,j=0,\dots ,\max (0,2^{l-1}-1).$$

$$Q_l^k:=\{\begin{array}{ll}{P}_l^k-P_{l-1}^k,& \text{if}l\ge 1,\\ P_0^k,& \text{if}l=0.\end{array}$$

$$W_l^k=Q_l^k(L^2(0,1)).$$

$$V_N^k=\underset{\begin{array}{c}0\le l\le N\end{array}}{⨁}{W}_l^k,P_N^k=\sum _{\begin{array}{c}0\le l\le N\end{array}}{Q}_l^k.$$

$$|\mathbf{\alpha }{|}_1:=\sum _{m=1}^d{\alpha }_m,|\mathbf{\alpha }{|}_{\infty }:=\underset{1\le m\le d}{\max }{\alpha }_m,$$

$$\mathbf{\alpha }\pm \mathbf{\beta }:=({\alpha }_1\pm {\beta }_1,\dots ,{\alpha }_d\pm {\beta }_d),2^{\mathbf{\alpha }}:=(2^{{\alpha }_1},\dots ,2^{{\alpha }_d}),$$

$$\mathbf{\alpha }\le \mathbf{\beta }\iff {\alpha }_m\le {\beta }_m,\forall m,$$

$$\max (\mathbf{\alpha },\mathbf{\beta }):=(\max ({\alpha }_1,{\beta }_1),\dots ,\max ({\alpha }_d,{\beta }_d)),\min (\mathbf{\alpha },\mathbf{\beta }):=(\min ({\alpha }_1,{\beta }_1),\dots ,\min ({\alpha }_d,{\beta }_d)).$$

$${\mathbf{W}}_{\mathbf{l}}^k:=W_{l_1,x_1}^k\times W_{l_2,x_2}^k\cdots \times W_{l_d,x_d}^k,$$

$${\mathbf{V}}_N^k:=V_{N,x_1}^k\times V_{N,x_2}^k\cdots \times V_{N,x_d}^k=\underset{|\mathbf{l}{|}_{\infty }\le N}{⨁}{\mathbf{W}}_{\mathbf{l}}^k,$$

$$v_{\mathbf{i},\mathbf{l}}^{\mathbf{j}}(\mathbf{x}):=\prod _{m=1}^d{v}_{i_m,l_m}^{j_m}(x_m),\mathbf{j}\in {\Theta }_{\mathbf{l}},1\le \mathbf{i}\le \mathbf{k}+1,$$

$${\hat{\mathbf{V}}}_N^k:=\underset{\begin{array}{c}|\mathbf{l}{|}_1\le N\end{array}}{⨁}{\mathbf{W}}_{\mathbf{l}}^k.$$

$${\mathbf{Q}}_{\mathbf{l}}^k:=Q_{l_1,x_1}^k\otimes \cdots \otimes Q_{l_d,x_d}^k,$$

$${\hat{\mathbf{P}}}_N^k:=\sum _{\begin{array}{c}|\mathbf{l}{|}_1\le N\end{array}}{\mathbf{Q}}_{\mathbf{l}}^k.$$

$$[q]=q^{-}{\mathbf{n}}^{-}+q^{+}{\mathbf{n}}^{+},$$

$$[\mathbf{q}]={\mathbf{q}}^{-}\cdot {\mathbf{n}}^{-}+{\mathbf{q}}^{+}\cdot {\mathbf{n}}^{+},$$

$$[q]=q\mathbf{n},\{\mathbf{q}\}=\mathbf{q},$$

$$B(u_h,v)=L(v),\forall v\in {\hat{\mathbf{V}}}_N^k$$

$$B(w,z)=$$

$$-$$

$$L(v)={\int }_{\Omega }fvd\mathbf{x}ds,$$

$${\left|\left|\left|v\right|\right|\right|}^2:=\sum _{\begin{array}{c}T\in {\Omega }_N\end{array}}{\int }_T|\nabla v{|}^{2}d\mathbf{x}+{\int }_{\Omega }c|u{|}^{2}d\mathbf{x}+\sum _{\begin{array}{c}e\in \Gamma \end{array}}h{\int }_e{\left\{\frac{\partial v}{\partial \mathbf{n}}\right\}}^{2}ds+\sum _{\begin{array}{c}e\in \Gamma \end{array}}\frac{1}{h}{\int }_e[v{]}^{2}ds,$$

$$|v{|}_{H^{\alpha ,L}\left(\Omega \right)}:={\Vert \left(\frac{{\partial }^{\alpha }}{\partial x_{i_1}^{\alpha }}\cdots \frac{{\partial }^{\alpha }}{\partial x_{i_r}^{\alpha }}\right)v\Vert }_{L^2\left(\Omega \right)},$$

$$|v{|}_{{\mathcal{H}}^{q+1}\left(\Omega \right)}:=\underset{1\le r\le d}{\max }\left(\underset{L\subset \{1,2,\dots ,d\}}{\max }|v{|}_{H^{q+1,L}\left(\Omega \right)}\right).$$

$$\Vert {\mathbf{Q}}_{\mathbf{l}}^k(v){\Vert }_{L^2(\Omega )}\lesssim 2^{-(q+1)|\mathbf{l}{|}_1}|v{|}_{{\mathcal{H}}^{q+1}(\Omega )}.$$

$$\left|\left|\left|{\hat{\mathbf{P}}}_N^k(v)-v\right|\right|\right|\lesssim N^d{2}^{-Nq}|v{|}_{{\mathcal{H}}^{q+1}(\Omega )}.$$

Full text

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 as approximations have been extensively studied for solving partial differential equations (PDEs) over the last few decades and become rather mature for various applications [13]. It is generally understood that DG methods are more flexible compared with continuous finite element methods due to the lack of continuity requirement, and thus enjoying many attractive properties. However, the situation changes when the dimension of the underlying problem becomes large. Traditional grid-based methods including DG methods suffer from the curse of dimensionality [6], which describes the scenario that the complexity and memory storage of an algorithm grows exponentially with the dimension of the underly problem for a given level of accuracy. Moreover, compared with other high order schemes, DG methods often require more degrees of freedom (DOF), and hence the effect of the curse of dimensionality is more significant, making such methods uncompetitive for high-dimensional calculations.

Driven by the need for a method that is able to retain the attractive properties of DG methods and yet requires feasible computational and storage cost for high-dimensional simulations, based on the sparse grid approach [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 solution is represented by a set of orthonormal tensor-product multiwavelet bases [1, 25, 17]. Meanwhile, following the standard sparse grid philosophy, rather than including all the anisotropic tensor-product bases, only the ones with significant contribution to the approximation accuracy are chosen, leading to immense reduction in cost complexity when the dimension $d$ is large. In particular, the number of DOF scales as $\mathcal{O}(h^{-1}|\log h|^{d-1})$ instead of $\mathcal{O}(h^{-d})$, where $h$ denotes the mesh size in each direction; while the approximation property of the proposed sparse grid DG method is proven only slightly deteriorated for smooth solutions through both theoretical and numerical verifications, see [25, 17].

In this paper, we are concerned with the Helmholtz equation with variable coefficients, which arises in many applications in science and engineering such as the representation of the solution of Schrödinger equation [22]. In [25], we developed an efficient high order sparse grid DG method based on the interior penalty DG (IPDG) weak formulation [5, 14, 26, 2, 3] for solving general elliptic equations, and such a method can be directly applied to simulate the problem. Due to the reduction of DOF, the resulting algebraic linear system enjoys a much smaller dimension than that by the traditional DG method, leading to computational savings 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 efficiency advantage of the sparse grid approach. To address the issue, we would like to explore the semi-orthogonality property associated with the orthonormal multiwavelet bases, which can be thought of as a special type of orthogonality 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 framework. In particular, by using the prewavelet as the bases, the standard bilinear formulation is carefully modified 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 the original method under some extra smoothness requirement of the coefficient. In this work, we demonstrate that the modification of the bilinear form based on semi-orthogonality is also effective for our sparse grid IPDG method with multiwavelet bases. In the theoretical aspect, we refine the analysis in [22] and show that the error incurred by modifying the bilinear form is one order higher than the projection error for sufficiently smooth problems, making the modified method produce virtually the same numerical result as the original sparse grid IPDG method; while the associated linear system is much sparser, offering extra computational and storage savings. Furthermore, under the sparse grid DG framework, the proposed methodology has the potential 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 IPDG on sparse grids for solving the Helmholtz equation developed in [25, 17]. In Section 3, we formulate the modified 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 performance of 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 the construction together with several key theoretical results of the IPDG method on sparse grids for solving the following Helmholtz equation,

[TABLE]

where $c$ is a non-negative function and the domain is the $d$-dimensional unit box.

2.1 DG Finite Element Spaces on Sparse Grids

In this subsection, we briefly review the recent development on the DG finite element approximation space on sparse grids in [25, 17], which plays a key role in dimension reduction for the proposed DG method.

We first introduce the hierarchical decomposition of DG approximation space in one dimension. Without loss of generality, we consider the domain $\Omega=[0,1]$ and define the $l$-th level grid $\Omega_{l}$, composed of non-overlapping intervals $I_{l}^{j}=(2^{-l}j,2^{-l}(j+1)]$, $j=0,\ldots,2^{l}-1$ with uniform cell size $h_{l}=2^{-l}.$ Denote by

[TABLE]

the collection of piecewise polynomials of degree at most $k$ defined on grid $\Omega_{l}$. Note that there exists a nested structure with respect to different values of $l$:

[TABLE]

which allows us to define the increment space $W_{l}^{k}$, $l=1,2,\ldots$ as the orthogonal complement of $V_{l-1}^{k}$ in $V_{l}^{k}$ with respect to the standard $L^{2}$ inner product on $\Omega$, i.e.

[TABLE]

By letting $W_{0}^{k}:=V_{0}^{k}$, which consists of all polynomials of up to degree $k$ on the entire domain $[0,1]$, we arrive at the hierarchical decomposition of the standard piecewise polynomial space $V_{N}^{k}$, $N\geq 1$ on grid $\Omega_{N}$ as

[TABLE]

In light of (2), we are able to represent a function in $V_{N}^{k}$ by making use of the orthonormal multiwavelet bases of space $W_{l}^{k}$ developed in [1]. We denote such orthonormal bases of $W_{l}^{k}$ as

[TABLE]

The details about the construction of the multiwavelet bases can be found in [1, 25].

It is worth noting that space $W_{l}^{k}$ can be constructed by means of projection operators as well, which will be helpful in the analysis below. Denote $P_{l}^{k}$ as the standard $L^{2}$ projection operator from $L^{2}(0,1)$ onto $V_{l}^{k}$, i.e., $V_{l}^{k}=P_{l}^{k}(L^{2}(0,1))$. Then, we define the increment projector,

[TABLE]

Evidently,

[TABLE]

To summarize, we have obtained the hierarchical decomposition of the piecewise polynomial space $V_{N}^{k}$ as well as the corresponding projection operator $P_{N}^{k}$:

[TABLE]

Then, we review the construction of approximation spaces in multi-dimensions. Consider the domain $\mathbf{x}=(x_{1},\ldots,x_{d})\in\Omega=[0,1]^{d}$. We first recall some notation on the norms and operations of multi-indices in $\mathbb{N}_{0}^{d}$, where $\mathbb{N}_{0}$ denotes the set of nonnegative integers. For ${\bm{\alpha}}=(\alpha_{1},\ldots,\alpha_{d})\in\mathbb{N}_{0}^{d},$ we define the $l^{1}$ and $l^{\infty}$ norms

[TABLE]

the component-wise arithmetic operations

[TABLE]

the relational operator

[TABLE]

and

[TABLE]

Denote $\mathbf{0}:=(0,\ldots,0)\in\mathbb{N}_{0}^{d}$ and $\mathbf{1}:=(1,\ldots,1)\in\mathbb{N}_{0}^{d}$.

Given a multi-index $\mathbf{l}=(l_{1},\ldots,l_{d})\in\mathbb{N}_{0}^{d}$ indicating the levels of the mesh in a multivariate sense, we define the grid $\Omega_{\mathbf{l}}$ consisting of non-overlapping elementary cells $I_{\mathbf{l}}^{\mathbf{j}}=\{\mathbf{x}:x_{m}\in(2^{-l_{m}}j_{m},2^{-l_{m}}(j_{m}+1)),\,m=1,\ldots,d\}$, $\mathbf{0}\leq\mathbf{j}\leq 2^{\mathbf{l}}-\mathbf{1}$. On grid $\Omega_{\mathbf{l}}$, we define the tensor-product increment space

[TABLE]

where $W_{l_{m},x_{m}}^{k}$ corresponds to the one-dimensional (1D) increment space in the $m$-th direction. Based on the 1D hierarchical decomposition (2), we yield

[TABLE]

where ${\bf V}_{N}^{k}$ is indeed the traditional tensor-product DG approximation space on grid $\Omega_{N}:=\Omega_{(N,\ldots,N)}$, i.e., having 2N cells in each dimension. ${\bf V}_{N}^{k}$ will be called the full grid space hereafter.

Define the index set $\Theta_{\mathbf{l}}=\{\mathbf{j}\in\mathbb{N}_{0}:\mathbf{0}\leq\mathbf{j}\leq\max(\mathbf{0},2^{\mathbf{l}-\mathbf{1}}-\mathbf{1})\}$. The basis functions for $\mathbf{W}_{\mathbf{l}}^{k}$ are constructed by a tensor product of the 1D multiwavelet bases

[TABLE]

where $\mathbf{k}:=(k,\ldots,k)\in\mathbb{N}^{d}_{0}$. Note that they form a set of orthonormal bases due to the property of the 1D bases. The sparse grid DG finite element approximation space is defined as

[TABLE]

Evidently, $\widehat{{\bf V}}_{N}^{k}$ is a subset of ${\bf V}_{N}^{k}$. More importantly, its number of DOF scales as $O((k+1)^{d}2^{N}N^{d-1})$, as opposed to exponential dependence on $Nd$ for the full grid space ${\bf V}_{N}^{k}$, see [25]. This is the key for computational savings in high dimensions.

The standard $L^{2}$ projection operator onto space $\mathbf{W}_{l}^{k}$ can be naturally written as

[TABLE]

where ${Q}^{k}_{l_{m},x_{m}}$ denotes the 1D projection operator ${Q}^{k}_{l_{m}}$ in the $m$-th dimension. Consequently, we write the $L^{2}$ projection operator onto the sparse approximation space $\widehat{{\bf V}}_{N}^{k}$ as

[TABLE]

2.2 IPDG Method on Sparse Grids

To formulate the IPDG method on sparse grids, we start with introducing some standard notation about jumps and averages for piecewise functions defined on grid $\Omega_{N}$. Let $\Gamma:=\bigcup_{T\in\Omega_{N}}\partial_{T}$ be the union of the boundaries for all the elements in $\Omega_{N}$ and $S(\Gamma):=\Pi_{T\in\Omega_{N}}L^{2}(\partial T)$ be the set of $L^{2}$ functions defined on $\Gamma$. For any $q\in S(\Gamma)$ and $\mathbf{q}\in[S(\Gamma)]^{d}$, their averages $\{q\},\{\mathbf{q}\}$ and jumps $[q],[\mathbf{q}]$ on the interior edges are defined as follows,

[TABLE]

where ${\bf n}$ denotes the unit normal, and ‘$-$’ and ‘+’ represent the directions of the underlying vector pointing to interior and exterior at an element edge $e$, respectively. If $e$ is a boundary edge, then we let

[TABLE]

where $\bm{n}$ is the outward unit normal.

The sparse grid IPDG scheme for (1) developed in [25] is defined as follows. We look for $u_{h}\in\widehat{{\bf V}}_{N}^{k}$, such that

[TABLE]

where

[TABLE]

and

[TABLE]

where $\sigma$ is a positive penalty parameter, $h=2^{-N}$ is the uniform mesh size in each dimension. Note that the bilinear form $B(\cdot,\cdot)$ is the same as in [26, 2], while instead of the expensive full grid piecewise polynomial space ${\bf V}_{N}^{k}$, the more efficient sparse grid approximation space $\widehat{\mathbf{V}}_{N}^{k}$ is employed, leading to a significant reduction in the size of the linear algebraic system especially when $d$ is large, see [25] for a detailed discussion.

Below, we summarize the theoretical results for the sparse grid DG methods (5). Define the energy norm of a function $v\in H^{2}(\Omega_{N})$ by

[TABLE]

where $H^{2}(\Omega_{N})$ denotes the standard broken Sobolev space on grid $\Omega_{N}$. We also need the following semi-norm. For any set $L=\left\{i_{1},\ldots,i_{r}\right\}\subset\{1,\ldots,d\}$, define $L^{c}$ to be the complement set of $L$ in $\{1,\ldots,d\}$. For a non-negative integer $\alpha$ and set $L$, we define the mixed derivative semi-norm for a function $v$

[TABLE]

and

[TABLE]

The space $\mathcal{H}^{q+1}(\Omega)$ denotes the closure of $C^{\infty}(\overline{\Omega})$ in the semi-norm $|\cdot|_{\mathcal{H}^{q+1}}$.

We first review several properties related to the projection operators introduced in the previous subsection. We refer the reader to [25, 17] for the proofs. To avoid unnecessary clutter of constants, the notation $A\lesssim B$ is used henceforth to represent $A\leq C\cdot B$, where the generic constant $C$ is independent of $N$ and the mesh level considered.

Lemma 2.1**.**

Let ${\bf Q}_{\mathbf{l}}^{k}$ be the $L^{2}$ projection operator onto the increment space $\mathbf{W}_{\mathbf{l}}^{k}$. For any $v\in\mathcal{H}^{p+1}(\Omega)$, $k\geq 1$, $1\leq q\leq\max(p,k)$, $N\geq 1$, $d\geq 2$, we have

[TABLE]

Lemma 2.2** (Projection error estimate.).**

Let $\widehat{{\bf P}}_{N}^{k}$ be the $L^{2}$ projection operator onto the space $\widehat{{\bf V}}_{N}^{k}$ introduced in (4). For any $v\in\mathcal{H}^{p+1}(\Omega)$, $k\geq 1$, $1\leq q\leq\max(p,k)$, $N\geq 1$, $d\geq 2$, we have

[TABLE]

The bilinear operator $B(\cdot,\cdot)$ is known to enjoy the following properties.

Lemma 2.3** (Orthogonality [25].).**

Let $u$ be the exact solution to (1), and $u_{h}$ be the numerical solution to (5), then

[TABLE]

Lemma 2.4** (Boundedness and stability [2, 4].).**

When $\sigma$ is taken large enough,

[TABLE]

In light of the Lemmas 2.2, 2.3 and 2.4, we can prove the following error estimate for the sparse grid IPDG method (5).

Theorem 2.5** (Convergence [25].).**

Let $u$ be the exact solution to (1), and $u_{h}$ be the numerical solution to (5). For $u\in\mathcal{H}^{p+1}(\Omega)$, $k\geq 1$, $1\leq q\leq\min\{p,k\}$, $N\geq 1$, $d\geq 2$, we have

[TABLE]

This theorem implies a convergence rate of $O(h^{k})$ up to the polylogarithmic term $|\log_{2}h|^{d}$ in the energy norm when $u$ is smooth enough.

Before 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 ${\bf P}^{k}_{n}$ be the $L^{2}$ projection operator onto the full grid space ${\bf V}_{n}^{k}$, $n=0,1\ldots,N$. Then, for any $v\in{\bf V}^{k}_{N}$, $k\geq 1$, $N\geq 1$, $d\geq 2$, we have

[TABLE]

where we let ${\bf P}^{k}_{-1}=0$, and $h_{n}=2^{-n}$ denote the mesh size of grid $\Omega_{n}$.

The proof is given in the Appendix A. By noting that

[TABLE]

we further have, for any $v\in\widehat{\mathbf{V}}_{N}^{k}$,

[TABLE]

3 Modified IPDG Method on Sparse Grids

In this section, we formulate a novel sparse grid IPDG method for solving the Helmholtz equation with variable coefficients. The key idea is to explore the semi-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\geq 0$ is constant, and two multiwavelet basis functions $v_{\mathbf{i},\mathbf{l}}^{\mathbf{j}}$ and $v_{\mathbf{i}^{\prime},\mathbf{l}^{\prime}}^{\mathbf{j}^{\prime}}\in\widehat{\mathbf{V}}_{N}^{k}$ satisfy $|\max(\mathbf{l},\mathbf{l}^{\prime})|_{1}>N$, then

[TABLE]

Proof: Since $v_{\mathbf{i},\mathbf{l}}^{\mathbf{j}}$, $v_{\mathbf{i}^{\prime},\mathbf{l}^{\prime}}^{\mathbf{j}^{\prime}}\in\widehat{\mathbf{V}}_{N}^{k}$, we have

[TABLE]

Together with $|\max(\mathbf{l},\mathbf{l}^{\prime})|_{1}>N$, we claim that there are two different dimension directions $m_{1}\neq m_{2}$ so that $l_{m_{1}}\neq l^{\prime}_{m_{1}}$ and $l_{m_{2}}\neq l^{\prime}_{m_{2}}$. Due to the orthonormal property of the multiwavelet bases, we follow the argument for the case of the prewavelet bases in [22] and prove that (9) holds.

$\blacksquare$

The semi-orthogonality property actually renders a highly desired sparse structure for the resulting stiffness matrix. On the other hand, when the method (5) is applied to (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 modify the bilinear formulation in light of semi-orthogonality, motivated by the work [21, 22, 19]. In particular, based on $B(\cdot,\cdot)$ introduced in (2.2), we define the following modified bilinear form $B^{so}(\cdot,\cdot):\widehat{\mathbf{V}}_{N}^{k}\times\widehat{\mathbf{V}}_{N}^{k}\rightarrow\mathbb{R}$

[TABLE]

The sparse grid IPDG weak formulation is modified accordingly as follows. We seek $u_{h}^{so}\in\widehat{\mathbf{V}}_{N}^{k}$ such that

[TABLE]

Note that, by construction, the resulting stiffness matrix by the modified method enjoys the same sparse structure as the constant-coefficient case, leading to computational savings. Meanwhile, when modifying the bilinear form by means of (10) we in fact commit a special type of variational crimes. Below, we show that the modified method will generate a numerical solution with the same order accuracy as the original method (5) under an extra smoothness assumption of $c$.

We need the following mixed derivative norm to measure the smoothness for the coefficient $c$. For a function $w$, define the norm

[TABLE]

and the space $\mathcal{W}^{\infty,k+1}(\Omega):=\{w\in L^{\infty}(\Omega):\|w\|_{\mathcal{W}^{\infty,k+1}(\Omega)}<\infty\}$. For convenience of illustration, we further introduce several shorthand notation. Whenever given two multi-indexes $\mathbf{l}$ and $\mathbf{l}^{\prime}$, we denote by $\mathbf{l}^{\max}=\max(\mathbf{l},\mathbf{l}^{\prime})$, $\mathbf{l}^{\min}=\min(\mathbf{l},\mathbf{l}^{\prime})$, and $\Theta_{\mathbf{l},\mathbf{l}^{\prime}}=\{s\in\{1,2,\ldots,d\}:l_{s}\neq l_{s}^{\prime}\}$.

We start with the following lemma, which will play a key role in the analysis.

Lemma 3.2**.**

Let $v_{\mathbf{i},\mathbf{l}}^{\mathbf{j}}$ and $v_{\mathbf{i}^{\prime},\mathbf{l}^{\prime}}^{\mathbf{j}^{\prime}}$ be the two multiwavelet basis functions of space $\mathbf{W}_{\mathbf{l}}^{k}$ and $\mathbf{W}_{\mathbf{l}^{\prime}}^{k}$, respectively. Denote

[TABLE]

Assume $c(\mathbf{x})\in\mathcal{W}^{\infty,k+1}(\Omega)$, then

[TABLE]

Proof: If $\mathcal{Q}_{\mathbf{l},\mathbf{l}^{\prime}}^{\mathbf{j},\mathbf{j}^{\prime}}$ is empty, then (12) is trivial. Below, we assume $\mathcal{Q}_{\mathbf{l},\mathbf{l}^{\prime}}^{\mathbf{j},\mathbf{j}^{\prime}}$ is nonempty. Since the multi-dimensional multiwavelet bases introduced in (3) are constructed by tensoring the 1D bases, we can rearrange $v_{\mathbf{i},\mathbf{l}}^{\mathbf{j}}(\mathbf{x})$ and $v_{\mathbf{i}^{\prime},\mathbf{l}^{\prime}}^{\mathbf{j}^{\prime}}(\mathbf{x})$ in the integrand according to $\mathbf{l}^{\max}$ and $\mathbf{l}^{\min}$, yielding

[TABLE]

Assume $\Theta_{\mathbf{l},\mathbf{l}^{\prime}}=\{m_{1},m_{2},\ldots,m_{r}\}$. Note that $l^{\max}_{s}\geq l^{\min}_{s}+1$ iff $s\in\Theta_{\mathbf{l},\mathbf{l}^{\prime}}$. In [17], we showed that

[TABLE]

Together with the elementary product rule and the fact that $v_{\mathbf{i}^{\min},\mathbf{l}^{\min}}^{\mathbf{j}^{\min}}$ is a piecewise polynomial of degree $k$ in each dimension, we obtain that

[TABLE]

A direct application of the inverse inequality leads to

[TABLE]

Also note that $\|v_{\mathbf{i}^{\max},\mathbf{l}^{\max}}^{\mathbf{j}^{\max}}\|_{L^{2}(\mathcal{Q}_{\mathbf{l},\mathbf{l}^{\prime}}^{\mathbf{j},\mathbf{j}^{\prime}})}=1$, and by definition,

[TABLE]

This, together with (13), (14) and (15), completes the proof.

$\blacksquare$

Lemma 3.3**.**

Assume $c(\mathbf{x})\in\mathcal{W}^{\infty,k+1}(\Omega)$. If $w(\mathbf{x})\in\mathbf{W}^{k}_{\mathbf{l}}$, and $z(\mathbf{x})\in\mathbf{W}^{k}_{\mathbf{l}^{\prime}}$, then

[TABLE]

Proof: Note that, for $w(\mathbf{x})\in\mathbf{W}^{k}_{\mathbf{l}}$ and $z(\mathbf{x})\in\mathbf{W}^{k}_{\mathbf{l}^{\prime}}$, we have

[TABLE]

Then, based on Lemma 3.2, the proof immediately follows the same procedure as in Lemma 4.6 in [22] and hence is omitted for brevity.

$\blacksquare$

The above lemma directly leads to the following estimate, which is useful when proving the boundedness and stability of the modified bilinear form $B^{so}(\cdot,\cdot)$.

Lemma 3.4**.**

Assume $c(\mathbf{x})\in\mathcal{W}^{\infty,k+1}(\Omega)$, $k\geq 1$. If $w(\mathbf{x})\in\mathbf{W}^{k}_{\mathbf{l}}$, and $z(\mathbf{x})\in\mathbf{W}^{k}_{\mathbf{l}^{\prime}}$, then

[TABLE]

Proof: Since $k\geq 1$,

[TABLE]

$\blacksquare$

Now, we are ready to establish the boundedness and stability of the modified bilinear form $B^{so}(\cdot,\cdot)$.

Theorem 3.5** (Boundedness and stability with semi-orthogonality.).**

Assume $c(\mathbf{x})\in\mathcal{W}^{\infty,k+1}(\Omega)$, $k\geq 1$, $d\geq 2$. There exists an integer $N_{0}$, such that

[TABLE]

where $w,z\in\widehat{\mathbf{V}}_{N}^{k}$ with $N\geq N_{0}$.

Proof: We first show that, for $w,z\in\widehat{\mathbf{V}}_{N}^{k}$,

[TABLE]

where

[TABLE]

Note that

[TABLE]

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 (8) to account for the discontinuous approximation space. For the completeness of the paper, we choose to provide the proof.

First, in [22], it was shown that if $|\mathbf{l}|_{1}\leq N$, $|\mathbf{l}^{\prime}|_{1}\leq N$, and $|\mathbf{l}^{\max}|_{1}>N$, then

[TABLE]

Then,

[TABLE]

where we have used the Cauchy-Schwarz inequality. Note that

[TABLE]

We bound the right-hand side of above equation as follows. First, we need the inequality proved in [22] that, for any $\mathbf{l}$ with $|\mathbf{l}|_{1}\leq N$,

[TABLE]

where $C$ is a constant independent of $N$ and $\mathbf{l}$. Together with the estimate (8), we derive that

[TABLE]

Combining (LABEL:eq:bilinear_diff), (LABEL:eq:bilinear_diff2), (19) and (20), we prove (16). Due to the boundedness and stability of $B(\cdot,\cdot)$, we complete the proof of the theorem.

$\blacksquare$

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 Helmholtz equation (1), and $u^{so}_{h}$ be the numerical solution to the modified IPDG formulation with semi-orthogonality (11). For $u\in\mathcal{H}^{p+1}(\Omega)$, $k\geq 1$, $c\in\mathcal{W}^{\infty,k+1}(\Omega)$, $q=\max(p,k)$, $d\geq 2$, $N\geq N_{0}$, where $N_{0}$ is given in Theorem 3.5, we have

[TABLE]

Proof: Following the theory of variational crime, see [16, 8], we consider the decomposition of the error

[TABLE]

For the second term on the right-hand side, due to stability of $B^{so}(\cdot,\cdot)$ in $\widehat{\mathbf{V}}_{N}^{k}$ from Theorem 3.5, we have

[TABLE]

where $u_{h}$ denotes the numerical solution by the original IPDG method (2.2). In the derivation, we have also used orthogonality and boundedness of $B(\cdot,\cdot)$. Note that the first term on the right-hand side is the projection error and has been estimated in Lemma 2.2; while the second term measures the effect of modifying the bilinear form. For the rest of the proof, we show that

[TABLE]

Since ${\bf Q}_{\mathbf{l}}^{k}\left(\widehat{{\bf P}}^{k}_{N}(u)\right)={\bf Q}_{\mathbf{l}}^{k}(u)$ for $|\mathbf{l}|_{1}\leq N$, we derive that

[TABLE]

where we have used Lemmas 3.3 and 2.1. Notice that

[TABLE]

Then, we have

[TABLE]

Similar to (LABEL:eq:bilinear_diff)-(20), we make use of the Cauchy-Schwarz inequality and bound the sum on the right-hand side as

[TABLE]

The proof of (23) is complete.

By combining (21), (22), (23), and Lemma 2.2 about the projection error estimate, we complete the proof.

$\blacksquare$

Remark 3.7**.**

$\sup_{w\in\widehat{\mathbf{V}}_{N}^{k}}\frac{B(\widehat{{\bf P}}^{k}_{N}u,w)-B^{so}(\widehat{{\bf P}}^{k}_{N}u,w)}{{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|w\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}}$* quantifies the variational crime from modifying the bilinear form and is indeed one order higher than the projection error. From the numerical results presented in the next section, we will see that if the coefficient $c$ and the solution $u$ are sufficiently smooth, then the modified 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 (1) with variable coefficients.*

Remark 3.8**.**

The proposed framework can be extended to other variable-coefficient problems if a similar estimate as in Lemma 3.2 is derived. For instance, for Poisson’s equation with variable coefficients $-\nabla\cdot(\mathbf{K}(\mathbf{x})\nabla u)=f$, under some condition of $\mathbf{K}$, it is possible to estimate

[TABLE]

and the boundary terms in the IPDG formulation with $|\mathbf{l}|_{1}\leq N$, $|\mathbf{l}|_{1}\leq N$ and $|\mathbf{l}^{\max}|_{1}>N$, and devise a modified sparse grid method with semi-orthogonality in an attempt to sparsify the stiffness matrix and save cost. We leave the investigation to the future work.

4 Numerical Results

In this section, we provide numerical results to demonstrate the performance of the modified sparse grid IPDG method as well as verify the analysis established for simulating the Helmholtz equations with variable coefficients. The penalty parameter $\sigma$ is set to be empirical values $\sigma=5\cdot k\cdot d$.

Example 4.1**.**

We solve the Helmholtz equation (1) with the smooth coefficient

[TABLE]

up to $d=6$. The right-hand $f$ is chosen such that

[TABLE]

This problem was solved in [19] by a continuous Ritz-Galerkin discretization on sparse grids using prewavelets in conjunction with semi-orthogonality. We first consider the case of $d=2$ and compare and contrast the performance of both modified and original sparse grid IPDG methods in terms of accuracy and efficiency. In Table 1, we summarize the convergence study for both methods with $k=1,\,2,\,3$. It is observed that, for this problem, the $L^{2}$ and $H^{1}$ 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^{2}$ error and $k$-th order for the $H^{1}$ error. Such an observation validates the convergence analysis in Theorem 3.6 that the variational crime incurred by the proposed modification of the bilinear form is one order higher than the projection error and hence becoming negligible in the simulation when the smoothness requirements are fulfilled. To demonstrate the efficiency advantage of the modified method, we provide the numbers of nonzero entries of the stiffness matrices as well as the conditional numbers for both methods with various mesh levels $N$ and polynomial degrees $k$ in Table 2. Further, in Figure 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 sparsity defined by $O_{s}=\log(\text{NNZ})/\log(\text{DOF})$, where NNZ is referred to as the number of nonzero elements. Note that NNZ scales as $O(\text{DOF}^{1.5})$ for the modified method; while it scales as $O(\text{DOF}^{1.7})$ for the original method. Such reduction leads to great computational savings when assembling as well as solving the algebraic linear system.

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 3 with $k=1,\,2,\,3,$ including the $L^{2}$ and $H^{1}$ errors and the associated orders of accuracy for both methods. As in the above $d=2$ case, both methods give almost the same numerical errors. $k+1$-th and $k$-th order of accuracy is observed for the $L^{2}$ and $H^{1}$ errors, respectively, which agrees with the error estimates established in Theorem 3.6. In Figure 2, we plot the sparsity patterns of the stiffness matrices by two methods with $N=4$ and $k=1$. As expected, the stiffness matrix by the modified method is much sparser than that by the original method. Furthermore, we provide DOF, NNZ and the associated order of sparsity $O_{s}$ for both methods in Table 4. It is observed that NNZ scales as $O(\text{DOF}^{1.7})$ for the modified method; while it scales as $O(\text{DOF}^{1.95})$ for the original method.

Last, we summarize the simulation results for $d=4,\,5,\,6$ in Tables 5-7, respectively. In particular, we report the $L^{2}$ errors, the associated orders of accuracy, NNZ, and the orders of sparsity $O_{s}$ for the modified method only. The observation is similar to that in the cases of $d=2,\,3$. The modified method is able to achieve the expected orders of convergence, while slight order reduction is observed due to the polylogarithmic term appearing in the error estimate. Note that unlike the modified method, the original method suffers the almost fully populated stiffness matrix, making its simulation exceed our computation resource limit. Currently, we are allowed to handle matrices with NNZ no greater than $1E+10$. Based on the comparison drawn in $d=2,\,3$ together with the Theorems 2.5 and 3.6, we believe the errors and the associated convergence rates of the two methods are also comparable.

5 Conclusion

In this paper, we developed a modified sparse grid IPDG method using semi-orthogonality for solving the Helmholtz equation with variable coefficients. The original IPDG features a sparse finite element space which scales as $O(h^{-1}|\log_{2}h|^{d-1})$ for $d$-dimensional problems, translating into a significant cost reduction when $d$ is large. On the other hand, when applied to the variable-coefficient problem, the method suffers a dense stiffness matrix, impeding its efficiency advantage to some extent. Based on the semi-orthogonality property associated with the orthonormal multiwavelet bases, the IPDG bilinear form is modified aiming to sparsify the resulting linear algebraic matrix. A numerical analysis demonstrates that the error incurred by the modification is one order higher than the projection error, making it negligible in the simulations. Numerical results in up to six dimensions are shown to validate the analysis and demonstrate that the modified method enjoys sparser stiffness matrix than the original one, leading to extra computational savings. Future work includes the study of adaptive algorithm for less regular solutions as well as extension of the method to other types of high-dimensional variable-coefficient equations.

Appendix A Proof of Lemma 2.6

In the proof, we need the averaging projection operator $\mathcal{A}$ proposed in [10, 11], which decomposes the space ${\bf V}_{N}^{k}$ into the conforming subspace ${\bf V}_{N}^{c}$ and nonconforming part ${\bf V}_{N}^{nc}$, i.e,

[TABLE]

Furthermore, $\mathcal{A}$ satisfies the Jackson estimate

[TABLE]

where the energy norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{a}$ is defined by

[TABLE]

and $\mathcal{I}$ denotes the identity operator. The construction and detailed analysis of $\mathcal{A}$ were established in [10, 11, 9]. Clearly, ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{a}\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$.

For any $v\in{\bf V}_{N}^{k}$, we denote by $w:=\mathcal{A}(v)$ and $z:=(\mathcal{I}-\mathcal{A})v$. We first prove that

[TABLE]

Note that there is a standard nested relation for conforming finite element spaces:

[TABLE]

where the space ${\bf V}^{c}_{n}$ defined on mesh $\Omega_{n}$ is a subspace of ${\bf V}^{k}_{n}$, $n=0,\ldots,N$. Following the idea of the BPX multilevel precondtioner [7], we define the Ritz projection operator ${\bf R}_{n}:{\bf V}^{c}_{N}\rightarrow{\bf V}^{c}_{n}$

[TABLE]

where

[TABLE]

and we let ${\bf R}_{-1}=0$. Denote $w_{n}:=({\bf R}_{n}-{\bf R}_{n-1})w\in{\bf V}_{n}^{c}$. Then, $w=\sum_{n=0}^{N}w_{n}$. Since ${\bf V}_{n}^{c}\subset{\bf V}_{n}^{k}$,

[TABLE]

In addition, the Ritz operator enjoys the property

[TABLE]

Together with the standard estimate of ${\bf R}_{n-1}$ by the conforming finite element analysis, we have

[TABLE]

Then, we derive

[TABLE]

By an elementary linear algebra result, see [8], together with the identity $a(w_{n},w_{m})=0$, $n\neq m$, we finally get

[TABLE]

and (26) is proved. Together with (26) and the Jackson estimate (25), we immediately obtain

[TABLE]

Now we derive the estimate for $z$. A direct application of the Poincaré-Friedrichs type inequality for space $\mathbf{V}_{N}^{k}$, see [2], yields

[TABLE]

This, together with the Jackson estimate (25), gives

[TABLE]

Combining the estimates (29) and (30), we complete the proof as

[TABLE]

$\blacksquare$