Dispersion optimized quadratures for isogeometric analysis

TL;DR

This paper introduces optimized quadrature blending schemes for isogeometric analysis that significantly reduce dispersion errors, enhance eigenvalue convergence, and improve accuracy in wave propagation and eigenvalue problems.

Contribution

The authors develop dispersion-optimized quadrature schemes for high-continuity isogeometric analysis, achieving superconvergence and improved robustness over existing methods.

Findings

Two extra orders of convergence in eigenvalue errors

Optimal convergence order in eigenfunction errors

Validated improved accuracy through numerical examples

Abstract

We develop and analyze quadrature blending schemes that minimize the dispersion error of isogeometric analysis up to polynomial order seven with maximum continuity in the span ($C^{p-1}$). The schemes yield two extra orders of convergence (superconvergence) on the eigenvalue errors, while the eigenfunction errors are of optimal convergence order. Both dispersion and spectrum analysis are unified in the form of a Taylor expansion for eigenvalue errors. As a consequence, the schemes increase the accuracy and robustness of isogeometric analysis for wave propagation as well as the differential eigenvalue problems. We analyze the methods' robustness and efficacy and utilize numerical examples to verify our analysis of the performance of the proposed schemes.

Tables

Table 1. Table 1: Leading order coefficients of the discrete dispersion relations for the different quadratures described and for different polynomial orders p = 4 , 5 , 6 , 7 𝑝 4 5 6 7 p=4,5,6,7 .

$p$	$G_{p+1}$	$G{L}_{p+1}$	$G_p$	$O_p$
4	$\frac{3}{20\cdot 9!}$	$\frac{79}{560\cdot 9!}$	$\frac{11}{70\cdot 9!}$	$\frac{317}{24\cdot 11!}$
5	$\frac{5}{12\cdot 11!}$	$\frac{29}{70\cdot 11!}$	$\frac{211}{504\cdot 11!}$	$\frac{35039}{420\cdot 13!}$
6	$\frac{691}{420\cdot 13!}$	$\frac{91177}{55440\cdot 13!}$	$\frac{5069}{3080\cdot 13!}$	$\frac{15479}{24\cdot 15!}$
7	$\frac{35}{4\cdot 15!}$	$\frac{105103}{12012\cdot 15!}$	$\frac{60061}{6864\cdot 15!}$	$\frac{91067}{15\cdot 17!}$

Table 2. Table 2: Eigenfunction H 1 superscript 𝐻 1 H^{1} -seminorm errors for u 3 subscript 𝑢 3 u_{3} when using C 1 superscript 𝐶 1 C^{1} quadratic isogeometric analysis with different quadrature rules.

$N$	$O_2$	$G_3$	$G{L}_3$	$G_2$
20	8.007409E-02	8.007620E-02	8.007355E-02	8.007966E-02
40	1.962886E-02	1.962889E-02	1.962885E-02	1.962894E-02
80	4.882983E-03	4.882983E-03	4.882983E-03	4.882984E-03
160	1.219233E-03	1.219233E-03	1.219233E-03	1.219233E-03
320	3.047138E-04	3.047138E-04	3.047138E-04	3.047138E-04
Order	2.01	2.01	2.01	2.01

Table 3. Table 3: Effectivity index (EI) of the error estimator R ​ ( λ ~ j h , u ~ j h ) 𝑅 superscript subscript ~ 𝜆 𝑗 ℎ subscript superscript ~ 𝑢 ℎ 𝑗 R(\tilde{\lambda}_{j}^{h},\tilde{u}^{h}_{j}) and the convergence order of EI to 1.

	$p=1$		$p=2$		$p=3$
$N$	${\lambda }_1$	${\lambda }_4$	${\lambda }_1$	${\lambda }_4$	${\lambda }_1$	${\lambda }_4$
5	1.00119	0.88632	0.96997	0.67227	0.92794	0.64385
10	1.00040	1.00127	0.99231	0.89264	0.97974	0.75699
20	1.00011	1.00192	0.99806	0.97046	0.99468	0.92242
40	1.00003	1.00057	0.99952	0.99242	1.00223	0.97899
Order	1.83	2.28	1.96	1.82	1.7	1.39

Table 4. Table 4: Wall-clock time (in seconds) required for the matrix assembly and the entire simulation of the eigenvalue problem. The test problem is ( 2.2 ) with Ω ¯ = [ 0 , 1 ] 2 ¯ Ω superscript 0 1 2 \overline{\Omega}=[0,1]^{2} , which is discretized uniformly with 100 × 100 100 100 100\times 100 elements.

	Gauss rule		Optimal rule
$p$	Assembly time	Total time	Assembly time	Total time
2	1.10	41.54	1.24	41.18
3	1.48	44.86	1.65	45.23
4	1.71	47.12	1.91	47.25
5	2.06	48.84	2.23	49.34
6	2.36	50.93	2.51	51.07
7	2.78	54.69	2.97	54.93

Equations

$$\Delta u+{\omega }^{2}u$$

$$u$$

$$-\Delta u=\lambda u\text{in}\Omega .$$

$$0<{\lambda }_1<{\lambda }_2\le {\lambda }_3\le \cdots$$

$$(u_j,u_k)={\delta }_{jk},$$

$$a(u,v)=\lambda b(u,v)\forall v\in H_0^1(\Omega ),$$

$$a(u_j,u_k)={\lambda }_{j}b(u_j,u_k)={\lambda }_j{\delta }_{jk}.$$

$$a(u^h,v_h)={\lambda }^{h}b(u^h,v_h)\forall v_h\in V_h.$$

$$B_a^0(\xi )$$

$$B_a^p(\xi )$$

$$|\lambda -{\lambda }^h|\le C{h}^{2p}\Vert u{\Vert }_{p+1,\Omega }^2\text{and}\Vert u-u^h{\Vert }_{1,\Omega }\le C{h}^p\Vert u{\Vert }_{p+1,\Omega },$$

$${\int }_{\hat{K}}\hat{f}(\hat{\mathbf{x}})\text{d}\hat{\mathbf{x}}\approx \sum _{l=1}^{N_q}{\hat{\varpi }}_l\hat{f}(\widehat{n_l}),$$

$${\int }_{K}f(\mathbf{x})\text{d}\mathbf{x}\approx \mathcal{Q}(f)=\sum _{l=1}^{N_q}{\varpi }_{l,K}f(n_{l,K}),$$

$${\int }_{K}f(\mathbf{x})\text{d}\mathbf{x}\approx {\mathcal{Q}}_{\tau }(f)=\tau {\mathcal{Q}}_1(f)+(1-\tau ){\mathcal{Q}}_2(f),$$

$${\tilde{a}}_h({\tilde{u}}^h,v_h)={\tilde{\lambda }}^h{\tilde{b}}_h({\tilde{u}}^h,v_h)\forall v_h\in V_h,$$

$${\tilde{a}}_h(w,v)=\sum _{K\in {\mathcal{T}}_h}\sum _{l=1}^{N_q}{\varpi }_{l,K}^{(1)}\nabla w(n_{l,K}^{(1)})\cdot \nabla v(n_{l,K}^{(1)})$$

$${\tilde{b}}_h(w,v)=\sum _{K\in {\mathcal{T}}_h}\sum _{l=1}^{N_q}{\varpi }_{l,K}^{(2)}w(n_{l,K}^{(2)})v(n_{l,K}^{(2)}),$$

$$\mathbf{K}\tilde{\mathbf{U}}={\tilde{\lambda }}^h\mathbf{M}\tilde{\mathbf{U}},$$

$$\Vert u_j-u_j^h{\Vert }_E^2={\lambda }_j\Vert u_j-u_j^h{\Vert }_{0,\Omega }^2+{\lambda }_j^h-{\lambda }_j.$$

$$\|u_{j}-\tilde{u}_{j}^{h}\|_{E}^{2}=\tilde{\lambda}_{j}^{h}-\lambda_{j}+\lambda_{j}\|u_{j}-\tilde{u}_{j}^{h}\|^{2}_{0,\Omega}+\|\tilde{u}^{h}_{j}\|_{E}^{2}-\|\tilde{u}^{h}_{j}\|_{E,h}^{2}+\lambda_{j}\Big{(}1-\|\tilde{u}^{h}_{j}\|^{2}_{0,\Omega}\Big{)}.$$

$$U(x)=\sum _{a=0}^N{U}^a{B}_a^p(x)$$

$${\tilde{B}}_h(U,v_h)=0,\forall v_h\in V_h,$$

$$(20+{\Lambda }^2)(U^{a-2}+U^{a+2})+(40+26{\Lambda }^2)(U^{a-1}+U^{a+1})-(120-66{\Lambda }^2)U^a=0,$$

$$U^a=e^{ia{\mu }_{G_3}^{(2)}h}$$

$$(20+{\Lambda }^2){\cos }^2({\mu }_{G_3}^{(2)}h)+(20+13{\Lambda }^2)\cos ({\mu }_{G_3}^{(2)}h)+(-40+16{\Lambda }^2)=0,$$

$${\mu }_{G_3}^{(2)}h=\Lambda -\frac{1}{1440}{\Lambda }^5-\frac{1}{6720}{\Lambda }^7+\mathcal{O}({\Lambda }^9).$$

$$(16+{\Lambda }^2)(U^{a-2}+U^{a+2})+(32+20{\Lambda }^2)(U^{a-1}+U^{a+1})-(96-54{\Lambda }^2)U^a=0,$$

$$(16+{\Lambda }^2){\cos }^2({\mu }_{G{L}_3}^{(2)}h)+(16+10{\Lambda }^2)\cos ({\mu }_{G{L}_3}^{(2)}h)+(-32+13{\Lambda }^2)=0.$$

$${\mu }_{G{L}_3}^{(2)}h=\Lambda +\frac{1}{2880}{\Lambda }^5-\frac{1}{16128}{\Lambda }^7+\mathcal{O}({\Lambda }^9).$$

$$(24+{\Lambda }^2)(U^{a-2}+U^{a+2})+(48+32{\Lambda }^2)(U^{a-1}+U^{a+1})-(144-78{\Lambda }^2)U^a=0,$$

$$(24+{\Lambda }^2){\cos }^2({\mu }_{G_2}^{(2)}h)+(24+16{\Lambda }^2)\cos ({\mu }_{G_2}^{(2)}h)+(-48+19{\Lambda }^2)=0.$$

Full text

††journal: arXiv

1 Introduction

Isogeometric analysis is a widely-used numerical method introduced by Hughes and his collaborators [1, 2] in 2005. Approximation, stability, and error estimates are established in [3]. Structural vibrations and wave propagation problems are investigated using isogeometric analysis in [4]. Spectrum analysis of the method shows that the method significantly improves accuracy in the spectral calculation over the classical finite element method [2]. Further advantages of the method on spectral approximation properties are investigated in [5].

A duality principle, which induces a bijective map from spectral analysis to dispersion analysis, is established [6]. The spectral analysis for structural vibrations (eigenvalue problems) and the dispersion analysis for wave propagation are then unified. Although the cost per degree of freedom of isogeometric analysis is higher than for finite elements [7, 8, 9, 10], the dispersion error is dramatically smaller than that of finite elements. In this paper, we propose a blending scheme for isogeometric analysis which increases the convergence by two orders with respect to the mesh size. With this motivation, under the framework of unified analysis, we study the dispersion optimization of the isogeometric analysis.

Dispersion analysis for the Galerkin finite element methods has been studied actively in literature; see for example Thomson and Pinsky [11, 12], Ihlenburg and Babuska [13], Ainsworth [14], and others [15, 16, 17]. In [11], Thomson and Pinsky study the dispersive effects of the Galerkin methods with different local approximation basis (Legendre, spectral, and Fourier elements) for the Helmholtz equation and it is found that the choice of the basis has a negligible effect on the dispersion errors. In [14], a $2p$ convergence rate of the dispersion error is shown for arbitrary $p$-th order finite element methods. For the linear case, a fourth order superconvergence result is obtained by a modified integration rule for finite elements in [18].

As early as 1984, Marfurt [19] conjectured that the most promising and cost-effective method for computational wave propagation is to employ a weighted average of the finite element and spectral element schemes. In 2010, Ainsworth and Wajid [20] introduced the optimal blending of these two schemes of arbitrary polynomial order. This optimal blending delivers two extra orders of convergence on the dispersion errors. In 2016, a dispersion optimized edge-based mimetic finite difference method for Maxwell’s equations in cold plasma was developed in [21]. The authors presented a generalized form of mass lumping and an optimization procedure to reduce the numerical dispersion error from second to fourth order accuracy. In [22], the authors described a similar result on dispersion by optimally blending finite element and spectral element methods for Maxwell’s equations.

To evaluate the bilinear form, in the case of spectral and finite elements, Ainsworth and Wajid pointed out that spectral elements use the Gauss-Lobatto quadrature rule while finite elements use the Gauss-Legendre quadrature rule [20] on the same polynomial space. The optimal blending of spectral and finite elements is essentially a blending of the Gauss-Legendre and Gauss-Lobatto quadrature rules. Thus, this new blended quadrature rule minimizes the dispersion errors.

In this paper, we minimize the dispersion error for isogeometric analysis by blending quadrature rules appropriately while the eigenfunction errors do not degenerate. We study optimally-blended quadratures for isogeometric analysis up to polynomial order seven. We minimize the dispersion error to obtain two additional orders of error convergence.

The rest of this paper is organized as follows. Section 2 presents the discretization of an eigenvalue problem and its generalized Pythagorean eigenvalue theorem. In Section 3, we present the discrete dispersion relation and derive the optimized dispersion error expansions for isogeometric analysis up to order seven. Dispersion and spectrum analysis are unified in the form of a Taylor expansion for the eigenvalue errors. Section 4 analyzes the error of the blending schemes and an a posteriori error estimator based on the superconvergence result of the eigenvalue error is derived. Section 5 presents numerical examples to demonstrate the performance of the proposed blending schemes. Section 6 describes our concluding remarks.

2 Problem setting

We consider stationary waves as described by the Helmholtz equation

[TABLE]

where $\Omega=(0,1)^{d},d\in\{1,2,3\},$ is a unit cube (unit interval when $d=1$ and unit square when $d=2$), $\Delta=\nabla^{2}$ is the Laplacian and $\omega=\omega_{f}/c$ with $\omega_{f}$ being the frequency of a particular sinusoidal wave and $c$ being the speed of sound of the medium. Denoting $\lambda=\omega^{2}$, in the view of duality and unified analysis of discrete approximations for wave propagation studied in [6], (2.1) is also posed as a second order elliptic eigenvalue problem: find real values $\lambda$ and nonzero functions $u$ such that

[TABLE]

The eigenvalue problem (2.2) has a countable set of eigenvalues $\lambda_{j}\in\mathbb{R}^{+}$ (c.f., [23, Sec. 9.8] and [24, 25, 26, 27])

[TABLE]

and an associated set of orthonormal eigenfunctions $u_{j}$

[TABLE]

where $(\cdot,\cdot)$ denotes the $L^{2}$-inner product on $\Omega$. Herein, the Kronecker delta is defined as $\delta_{lm}=1$ when $l=m$ while zero otherwise.

2.1 Discretization

For an open bounded set $S\subset\mathbb{R}^{d}$ with Lipschitz boundary, we denote by $H^{m}(S)$ the Sobolev spaces and $H^{m}_{0}(S)$ the Sobolev spaces with functions vanishing at the boundary, where $m>0$ specifies the order of the weak derivative. The variational formulation of (2.2) is to find $\lambda\in\mathbb{R}^{+}$ and $u\in H^{1}_{0}(\Omega)$ such that

[TABLE]

where $a(w,v)=(\nabla w,\nabla v)$ and $b(w,v)=(w,v)$. Let $(\lambda_{j},u_{j})$ be an eigenpair, then $a(u_{j},v)=\lambda_{j}b(u_{j},v)$. For orthonormal eigenfunctions $u_{j}$ (in the sense of $L^{2}$-inner product; see (2.4)), they are also orthogonal in the energy inner product

[TABLE]

Let $\prod$ be the product symbol and $N_{k}$, $k=1,\cdots,d,$ be positive integers associated with the space variable $x_{k}$. Assume that a uniform tensor product mesh with $\prod_{k=1}^{d}N_{k}$ elements is placed on $\overline{\Omega}=[0,1]^{d}$ with grid nodes located at $(h_{1}n_{1},\cdots,h_{d}n_{d})$, where $h_{k}=\frac{1}{N_{k}},k=1,\cdots,d,$ is the size of the $k^{\text{th}}$ dimension and $n_{k}=0,1,\cdots,N_{k}.$ In case of one dimension, we simplify the notation as $N,h,$ respectively. We denote each element as $K$ and their collection as $\mathcal{T}_{h}$ such that $\bar{\Omega}=\cup_{K\in\mathcal{T}_{h}}K$. Due to the tensor product structure of the discretization, the element size is $h=\sqrt{\sum_{k=1}^{d}h_{k}^{2}}$. The Galerkin-type numerical methods seek $\lambda^{h}\in\mathbb{R}^{+}$ and $u^{h}\in V_{h}\subset H_{0}^{1}(\Omega)$ such that

[TABLE]

Different solution/trial spaces $V_{h}$ lead to different numerical methods. Under this framework, we consider the following methods: finite element, spectral element, and isogeometric analysis. We utilize the B-spline basis functions for isogeometric analysis. Following [28, 2, 29], the definition of the $p$-th order B-spline basis functions in one dimension is as follows. Let $\Xi=\{\xi_{0},\xi_{1},\cdots,\xi_{N_{0}}\}$ be an ordered knot vector with $0=\xi_{0}\leq\xi_{1}\leq\cdots\leq\xi_{N_{0}}=1$, that is, a nondecreasing sequence of real numbers called knots. The B-spline basis function of degree $p$, denoted as $B_{a}^{p}(\xi)$ (here $a$ refers to an index with slightly abuse of notation $a(\cdot,\cdot)$), is defined as

[TABLE]

In this paper, for isogeometric analysis, we utilize B-splines on uniform tensor product meshes with non-repeating knots, that is, the B-splines with maximum continuity, while for finite element method, we utilize $C^{0}$ B-spline basis functions. For multiple dimensions, the B-spline basis functions are constructed by tensor products of these one-dimensional B-spline basis functions; we refer to [2, 29] for details. Let $B^{p}_{n_{k}}$ be the one-dimensional basis functions in dimension $k=1,\cdots,d$ for $p$-th order finite element and isogeometric analysis. Provided the tensor product structure of the discretization, a $d$-dimensional basis function can be written as $\prod_{k=1}^{d}B^{p}_{n_{k}}$. Then $V_{h}=\text{span}\{\prod_{k=1}^{d}B^{p}_{n_{k}}\}_{n_{k}=0,1,\cdots,N_{k}}$.

In the framework of finite elements, the eigenpairs $(\lambda^{h},u^{h})$ have the following properties (see for example, [25, 27, 30, 31])

[TABLE]

where $C$ is a constant independent of $h$.

In practice, the integrals involved in $a(u^{h},v_{h})$ and $b(u^{h},v_{h})$ are evaluated numerically, that is, approximated by quadrature rules. On a reference element $\hat{K}$, a quadrature rule is of the form

[TABLE]

where $\hat{\varpi}_{l}$ are the weights, $\hat{n_{l}}$ are the nodes, and $N_{q}$ is the number of quadrature points. For each element $K$, we assume that there is an invertible affine map $\sigma$ such that $K=\sigma(\hat{K})$, which leads to the correspondence between the functions on $K$ and $\hat{K}$. Assuming $J_{K}$ is the corresponding Jacobian of the mapping, (2.10) induces a quadrature rule, denoted as $\mathcal{Q}$, over the element $K$ given by

[TABLE]

where $\varpi_{l,K}=\text{det}(J_{K})\hat{\varpi}_{l}$ and $n_{l,K}=\sigma(\hat{n}_{l})$. Let $G_{N_{q}}$ and $GL_{N_{q}}$ denote the $N_{q}$-point Gauss-Legendre quadrature rule and Gauss-Lobatto quadrature rule, respectively. The detailed description of these rules is given in [32]. We define the blending quadrature rule $\mathcal{Q}_{\tau}$ as

[TABLE]

where $\tau$ is the blending parameter and $\mathcal{Q}_{1},\mathcal{Q}_{2}$ are different quadrature rules. In this paper, we seek the optimal blending parameter to minimize the dispersion errors of the isogeometric analysis. We denote by $O_{p}$ the optimal blending scheme for the $p$-th order isogeometric analysis.

Applying quadrature rules to (2.7), we have the approximated form

[TABLE]

where

[TABLE]

and

[TABLE]

where $\{\varpi_{l,K}^{(1)},n_{l,K}^{(1)}\}$ and $\{\varpi_{l,K}^{(2)},n_{l,K}^{(2)}\}$ specify two (possibly different) quadrature rules.

We represent the eigenfunctions as a linear combination of the B-spline (or Lagrange) basis functions and substitute all the basis functions for $V_{h}$ in (2.13). Applying quadrature rules, this leads to the generalized matrix eigenvalue problem

[TABLE]

where $\mathbf{K}_{ab}=\tilde{a}_{h}(B_{a},B_{b}),\mathbf{M}_{ab}=\tilde{b}_{h}(B_{a},B_{b})$ with $B_{a}$ and $B_{b}$ being generic basis functions, and $\tilde{\mathbf{U}}$ is the corresponding representation of the eigenvector as the coefficients of the basis functions. Once the matrix eigenvalue problem (2.16) is solved, we obtain eigenpairs. Throughout this paper, we refer to $(\lambda_{j},u_{j})$ as one exact eigenpair, $(\lambda_{j}^{h},u_{j}^{h})$ as one approximate eigenpair when the inner products are not modified by the numerical quadrature, and $(\tilde{\lambda}_{j}^{h},\tilde{u}_{j}^{h})$ one approximate eigenpair when modified by the numerical quadrature.

2.2 Pythagorean eigenvalue theorem and its generalization

Denoting the energy norm as $\|\cdot\|_{E}=\sqrt{a(\cdot,\cdot)}$, the Pythagorean eigenvalue theorem [27, 5] states

Lemma 1**.**

For each discrete mode, with the normalization $\|u_{j}\|_{0,\Omega}=1$ and $\|u_{j}^{h}\|_{0,\Omega}=1$, there holds

[TABLE]

In general, applying quadrature rules to the inner products $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$ results in quadrature errors. We denote by $\|\cdot\|_{E,h}=\sqrt{\tilde{a}_{h}(\cdot,\cdot)}$ the (approximate) energy norm evaluated by a quadrature rule (2.14), the generalized Pythagorean eigenvalue theorem, c.f., [33], is stated as follows

Theorem 1**.**

For each discrete mode, with the normalization $\|u_{j}\|_{0,\Omega}=1$ and $\tilde{b}_{h}(\tilde{u}_{j}^{h},\tilde{u}_{j}^{h})=1$, there holds

[TABLE]

Remark 1**.**

The Pythagorean eigenvalue theorem reveals the properties of the numerical approximation of the eigenvalue problem when the inner products are fully represented in the discrete setting, while the generalized theorem accounts for the effect of the numerical integration.

3 The discrete dispersion relation and dispersion optimization

In the linear element case, isogeometric analysis and finite elements coincide and result in the same discrete dispersion relation; see for example [14, 34, 20]. We start with quadratic elements and focus on the one-dimensional case with $\Omega=[0,1]$ to illustrate the idea. We seek an approximate solution of the form

[TABLE]

satisfying

[TABLE]

where $\tilde{B}_{h}(w,v)=\tilde{a}_{h}(w,v)-\lambda\tilde{b}_{h}(w,v).$

3.1 Quadratic elements

We consider the $C^{1}$ quadratic B-spline basis function of the isogeometric analysis (for quadratic finite elements, we refer to [14, 20]). Applying quadrature rule $G_{3}$, one obtains the following equation for the value $U^{a}$ of the approximation at node $x_{a}=ah$ where $a=3,4,\cdots,N-2$ (that is, a node away from boundary),

[TABLE]

where $\Lambda=\sqrt{\lambda}h$. Let $i^{2}=-1$. This equation admits nontrivial solutions of the form

[TABLE]

provided that $\mu^{(2)}_{G_{3}}$ satisfies

[TABLE]

which is known as the discrete dispersion relation (c.f., [14]) for the discrete method with a particular quadrature rule. Solving (3.5) for $\mu^{(2)}_{G_{3}}h$ and writing the expression as a series in $\Lambda$ (assuming $\Lambda=\sqrt{\lambda}h<1$), we obtain the discrete dispersion error

[TABLE]

Now applying $GL_{3}$, one obtains

[TABLE]

which leads to the discrete dispersion relation

[TABLE]

Solving (3.8) for $\mu^{(2)}_{GL_{3}}h$ and writing the expression as a series in $\Lambda$, we obtain

[TABLE]

Lastly, we apply $G_{2}$ to obtain

[TABLE]

which leads to the discrete dispersion relation

[TABLE]

Solving (3.11) for $\mu^{(2)}_{G_{2}}h$ and writing the expression as a series in $\Lambda$, we obtain

[TABLE]

Remark 2**.**

$G_{3}$* integrates both the stiffness and mass terms exactly while $GL_{3}$ and $G_{2}$ integrate the stiffness terms exactly but under-integrate the mass ones. However, all of them yield the optimal order $h^{2p}$ of convergence. Similar calculations also indicate that once the stiffness term is under-integrated, for example by $GL_{2}$ or $G_{1}$, the optimal convergence is lost, which is expected in the view of Strang’s second lemma. The impact of under-integration on convergence is discussed in [27] among others.*

The leading coefficient $-\frac{1}{1440}$ in (3.6) and the leading coefficient $\frac{1}{2880}$ in (3.9) allows us to conjecture that the blending rule $\frac{1}{3}G_{3}+\frac{2}{3}GL_{3}$ eliminates the leading term in the dispersion error since there holds $\frac{1}{3}\cdot(-\frac{1}{1440})+\frac{2}{3}\cdot\frac{1}{2880}=0.$ Hence, the signs and coefficients of (3.6), (3.9), and (3.12) allow us to propose the following blendings

[TABLE]

All of them lead to the discrete dispersion relation for the optimal scheme

[TABLE]

which gives the optimal dispersion error

[TABLE]

Remark 3**.**

The blending rules allow us to eliminate the fifth order error and give a seventh order error. Moreover, the coefficient of the seventh order error is the sum of the corresponding ones in (3.6), (3.9), and (3.12) with weights specified in (3.13). However, calculation shows that this is not the case for the coefficients of the ninth or higher orders.

To show that the blending is optimal, applying the general blending rule

[TABLE]

gives the dispersion error

[TABLE]

where $\tau=1/3$ eliminates the fifth order error term, hence it is the optimal blending parameter. Similarly, one can verify the optimal blending parameters in (3.13). Alternatively, one can optimize the dispersion error directly from the error expansion in the most general form for the mass

[TABLE]

where $\alpha$ and $\beta$ are parameters representing the approximated mass entries. Herein, partition of unity requires the coefficient of $U^{a}$ in the mass term in (3.18) to be $(1-2\alpha-2\beta)$. This leads to a dispersion error expansion

[TABLE]

To optimize the dispersion error, we set and solve

[TABLE]

to obtain

[TABLE]

which in return simplifies (3.19) to the optimized dispersion error expression

[TABLE]

Thus, this alternative way leads to the same optimized dispersion error. In fact, values of $\alpha$ and $\beta$ uniquely determine the coefficients in the blending schemes (3.13).

3.2 Cubic elements

Now we consider the cubic B-spline element case with maximum continuity, that is, $C^{2}$ basis functions. Similarly, applying $G_{4},GL_{4}$, and $G_{3}$, we obtain the following discrete dispersion relations

[TABLE]

They lead to the dispersion error expansions given below

[TABLE]

Following the procedure for the quadratic case, we obtain the corresponding optimal blending schemes

[TABLE]

which all result in the following optimized dispersion error expression

[TABLE]

Remark 4**.**

As in the quadratic case, one can assign $\alpha,\beta$, and $\gamma$ to the mass terms to derive the optimized dispersion error expression. There are other blending schemes which give the same optimized dispersion error expression. Here, we list a few

[TABLE]

Herein, blending more quadrature rules does not reduce the dispersion and eigenvalue errors further.

3.3 Higher-order elements

For higher order elements, the calculations and derivations become more complicated. In the following, we list our results up to order seven for isogeometric analysis with maximum continuity basis functions. As before, for orders $p=4,5,6,7$, we apply $G_{p+1},GL_{p+1}$, and $G_{p}$ to compute the integrals. Their corresponding dispersion error expressions are of the form

[TABLE]

for $Q=G_{p+1},GL_{p+1},G_{p}$, while the optimized dispersion expressions are of the form

[TABLE]

In Table 1, we list the coefficients $\epsilon^{Q}_{p}$ and $\epsilon_{p}^{o}$ for $p=4,5,6,7$. As before, the optimal blending schemes are not unique; we list a few of them below

[TABLE]

These blending schemes can be rewritten as

[TABLE]

where $C_{1,p}$ and $C_{2,p}$ depend on $p$ and satisfy

[TABLE]

Thus, once the general form of one of them is obtained, the other one can be derived from (3.32).

Remark 5**.**

The analysis of the possible generalization of the optimal blending schemes to arbitrary order $p$ with variable continuities $C^{k},k=1,\cdots,p-1,$ is an open question and is still under investigation.

3.4 Extension to multiple dimensions

The extension is similar to the extension done for finite elements in [20]. The conclusion is that the optimal blending for arbitrary dimension coincides with the one we derive above for the one-dimensional case and is independent of the number of spatial dimensions.

In multiple dimensions, a $d$-dimensional basis function is a tensor product of one-dimensional basis functions, that is, $\prod_{k=1}^{d}B_{n_{k}}$ (we drop $p$ here as the following derivation is independent of the polynomial order $p$). We seek a solution of the form

[TABLE]

where each $U_{k}(x_{k})$ has a one-dimensional representation (3.1), which is further written by using (3.4) as

[TABLE]

with $\lambda_{k}$ (only in this subsection) being the squared wave number in the dimension $k=1,\cdots,d$. This allows us to also write

[TABLE]

where $U_{n_{1},\cdots,n_{d}}$ with $n_{k}=0,1,\cdots,N_{k},k=1,\cdots,d$ are the coefficients of the linear combination of the multi-dimensional basis functions.

From (3.33), we have

[TABLE]

Suppose that in each dimension $k$, $U_{k}(x_{k})$ satisfies (3.2), then

[TABLE]

By applying the one-dimensional equations (3.37) with $k=1,\cdots,d$, we obtain

[TABLE]

The multi-dimensional problem is to find $U({\boldsymbol{x}})$ satisfying for all the basis functions

[TABLE]

Plugging (3.36) and (3.38) into (3.39), the multidimensional problem admits a nontrivial solution of the form $\eqref{eq:mdbsol}$ provided that

[TABLE]

Replacing the $\lambda$ in the expression (3.2) for one dimension and (3.39) for multiple dimensions with the numerical approximated $\tilde{\lambda}^{h}_{k}$ and $\lambda^{h}$, respectively, the same derivations lead to

[TABLE]

Thus, subtracting (3.40) from (3.41), the dispersion error for multidimensional problems consists of the dispersion errors for each dimension, that is,

[TABLE]

which means that the numerical schemes for the multi-dimensional problem (2.2) preserve the error estimations of the one-dimensional problem when using the tensor product structure for the mesh discretization.

3.5 Duality with spectrum analysis

There is a symmetry between (2.1) and (2.2). The duality between spectrum analysis and dispersion analysis has been established in [6]. We denote the approximate squared frequencies by $\tilde{\lambda}^{p,h}_{Q}$, where $Q$ specifies the quadrature rule. For spectrum analysis in the quadratic case, we obtain

[TABLE]

where the first equation is also given in [6].

Applying the Taylor expansion on these expressions gives

[TABLE]

Similarly, this spectrum analysis can be done for cubic and higher order methods. In general, for multiple dimensions, in the view of (3.42), we obtain error expressions for $p=2,\cdots,7$ in the forms

[TABLE]

for $Q=G_{p+1},GL_{p+1},G_{p}$ and

[TABLE]

for the optimal schemes. For the cubic case, the coefficients $\epsilon_{3}^{o}=1/145152$ and $\epsilon^{Q}_{3}=1/60480$, $1/100800,13/604800$ for $Q=G_{3},GL_{3}$, and $G_{2}$, respectively, while for $p=4$, $5,6,7$, coefficients $\epsilon^{Q}_{p}$ and $\epsilon_{p}^{o}$ are given in Table 1.

Remark 6**.**

Equations (3.45) with (3.28) and (3.46) with (3.29) reveal the duality principle of dispersion and spectrum analysis in the error expansion form. The duality principle remains valid for the optimal blending schemes. The different signs of the coefficients $\epsilon^{Q}_{p}$ and $\epsilon_{p}^{o}$ in the error expressions are consequences of duality.

4 Error analysis

In the framework of finite element analysis, the eigenvalue and eigenfunction errors for both $(\lambda_{j}^{h},u_{j}^{h})$ and $(\tilde{\lambda}_{j}^{h},\tilde{u}_{j}^{h})$ are of the forms (2.9), c.f., [27, 30]. We now estimate the errors of the approximated eigenpair $(\tilde{\lambda}_{j}^{h},\tilde{u}_{j}^{h})$ under the framework of isogeometric analysis.

4.1 Eigenvalue estimates

The eigenvalue errors when using the standard quadrature rules, such as Gauss and Lobatto rules, converge at rate $h^{2p}$, c.f., [31, 4, 5]. This is also confirmed by our theoretical finding (3.45) for orders $p=2,\cdots,7$ and we conjecture that this is true for arbitrary order. To see this, squaring both sides of (3.45) and using $\Lambda=\sqrt{\lambda}h$ gives the following estimate

[TABLE]

which reduces to

[TABLE]

where $C$ is a constant independent of $h$. Clearly, this yields the $2p$ order convergence for the eigenvalues. Now, for isogeometric analysis with order up to $p=7$, we have the following eigenvalue estimate

Theorem 2**.**

For a fixed $j\geq 1$ and $p=1,\cdots,7$, assume $\tilde{\lambda}_{j}^{h}=\tilde{\lambda}^{p,h}_{O_{p}}$. For $h$ such that $\max{\{\sqrt{\lambda_{j}}h,\sqrt{\tilde{\lambda}_{j}^{h}}h\}}<1$, we have

[TABLE]

where $C$ is a constant independent of $h$.

Proof.

For order $p=1$, this is shown in [20]. For $p=2,\cdots,7,$ we have the error representation (3.46). For a fixed $j$, squaring both sides of (3.46) and using $\Lambda=\sqrt{\lambda_{j}}h$ gives the following estimate

[TABLE]

which reduces to the desired result. ∎

This theorem shows that the optimally-blended schemes produce two extra orders of convergence for the eigenvalue errors.

4.2 Eigenfunction estimates

In this section, we establish the optimal convergence rates of the eigenfunction errors. First, we establish the coercivity of the bilinear forms (2.14) and (2.15).

Lemma 2**.**

Given the blending scheme of the form (3.31), the bilinear forms (2.14) and (2.15) are coercive, that is, there holds

[TABLE]

where $\tilde{\alpha},\tilde{\beta}>0$ are constants independent of $h$.

Proof.

Without loss of generality, we consider the blending scheme $G_{p+1}+C_{2,p}\cdot(G_{p+1}-G_{p})$ for $p$-th order isogeometric analysis. Since $G_{p+1}$ and $G_{p}$ are exact for polynomial spaces of order $2p+1$ and $2p-1$, respectively, the blending scheme is exact for polynomial space of order $2p-1$. Hence, by Theorem 4.1.2 in [35], the bilinear form $\tilde{a}_{h}(\cdot,\cdot)$ is coercive.

Now, for $\tilde{b}_{h}(v_{h},v_{h})$, if the quadrature rule applied to the integral is $G_{p+1}$, then $\tilde{b}_{h}(v_{h},v_{h})=\|v_{h}\|^{2}_{0,\Omega}.$ In the following, we treat the quadrature rules as operators on integrals and let $\mathcal{I}$ be the identity operator. We calculate

[TABLE]

By Lemma 3.1 in [30], we have

[TABLE]

where $C$ is a constant independent of $h$.

Thus, for sufficiently small $h$, we have a constant $\tilde{\beta}>0$ such that

[TABLE]

which completes the proof. ∎

The boundedness of $\tilde{b}_{h}(\cdot,\cdot)$ is a consequence of the proof. With coercivity, we immediately conclude that the eigenvalues of (2.13) are positive. Before we establish the eigenfunction error estimate, we present the following inequality, which can be obtained by applying the Aubin-Nitsche Lemma (duality argument) on the discrete solution operator $T_{h}:L^{2}(\Omega)\to V_{h}$ defined as $\tilde{a}_{h}(T_{h}(\tilde{u}^{h}),v_{h})=\tilde{\lambda}^{h}\tilde{b}_{h}(\tilde{u}^{h},v_{h}),\forall\ v_{h}\in V_{h}$ as in (2.13), c.f., [35, Theorem 3.2.4] or [36, Section 2.3.4].

Lemma 3**.**

Suppose $\Omega=[0,1]^{d}\subset\mathbb{R}^{d}.$ For $u_{j},\tilde{u}_{j}^{h}\in H^{1}_{0}(\Omega)$, there holds

[TABLE]

where $C$ is a constant independent of $h$.

Theorem 3**.**

For a fixed $j\geq 1$, assume that $u_{j}$ and $\tilde{u}_{j}^{h}$ are normalized, that is, $b(u_{j},u_{j})=1$ and $\tilde{b}_{h}(\tilde{u}_{j}^{h},\tilde{u}_{j}^{h})=1$, and the signs of eigenfunctions of $u_{j}$ and $\tilde{u}_{j}^{h}$ are chosen such that $b(u_{j},\tilde{u}_{j}^{h})>0$. Then for sufficiently small $h$, we have the estimate

[TABLE]

where $C$ is a constant independent of $h$.

Proof.

With the normalization $b(u_{j},u_{j})=1$ and $\tilde{b}_{h}(\tilde{u}_{j}^{h},\tilde{u}_{j}^{h})=1$, by (2.5) and (2.13) we have

[TABLE]

Analogously to the derivations described in [30, 37, 38], we estimate the modified bilinear forms, that is, for $v,w\in V_{h}$ there holds

[TABLE]

Thus, for a fixed $j$, by definition and using (2.5) and (2.13), we obtain

[TABLE]

Rearranging terms and taking absolute value yields

[TABLE]

By applying Lemma 3, Theorem 2, (4.6), and (4.8), for sufficiently small $h$, we have a constant $\tilde{C}>0$ such that

[TABLE]

Taking the square root on both sides completes the proof. ∎

5 Numerical experiments

In this section, we present various numerical examples to illustrate the performance of the optimally-blended rules. We focus on isogeometric analysis and study the convergence of the eigenvalue (EV) and eigenfunction (EF) errors in one and two dimensions. Three or higher dimensions are simple extensions as discussed in Section 3.4. For comparisons with finite element approximations, we refer to [4, 5, 33] where the authors showed that the isogeometric analysis outperforms finite element approximations.

5.1 Convergence study in one dimension

We consider the classic one-dimensional second-order elliptic eigenvalue problem (2.2) with eigenvalues and eigenfunctions $\lambda_{j}=\pi^{2}j^{2},u_{j}=\sqrt{2}\sin(j\pi x),j=1,2,\cdots.$ This is the same example studied in [5, 33] among many others. We assume that once the eigenvalue problem is solved, the numerical eigenvalues $\lambda_{j}^{h}$ and $\tilde{\lambda}_{j}^{h}$ are sorted in ascending order and paired with the true eigenvalues $\lambda_{j}$.

First of all, for $p$-th order isogeometric analysis with maximum continuity, that is $C^{p-1}$, the quadrature rule $G_{p+1}$ approximates both stiffness and mass matrices exactly while $GL_{p+1}$, $G_{p}$, and $O_{p}$ integrate the stiffness matrices exactly but under-integrate the mass matrices. Despite these differences, all of them lead to accurate approximations to the eigenfunctions. Figure 2 shows the plots of the numerical approximations of eigenfunctions $u_{1}$ and $u_{10}$ using $C^{1}$ quadratic isogeometric analysis with two elements for $u_{1}$ and twenty elements for $u_{10}$. As a comparison, Figure 2 shows the numerical approximation of $u_{10}$ using $C^{2}$ cubic isogeometric analysis with twenty elements.

5.1.1 Eigenvalue errors

The analysis of eigenvalue errors done in Section 4.1 is verified numerically in this subsection. The optimally-blended quadrature rules proposed for isogeometric analysis of eigenvalue problem (2.2) in Section 3 yield two additional orders of eigenvalue error convergence.

Figure 4 shows the eigenvalue errors when using $C^{1}$ quadratic isogeometric elements. The domain $\overline{\Omega}=[0,1]$ is discretized uniformly with 20, 40, 80, 160, and 320 elements. The figure shows two extra orders of convergence in the error for the eigenvalues $\lambda_{3}$ and $\lambda_{11}$. With both stiffness and mass matrices integrated exactly by $G_{3}$, the eigenvalue errors converge at the rate of $h^{4}$, while a slight modification of the quadrature rule, which is easily realized by optimally-blended quadrature rule $O_{2}$, leads to a convergence rate of $h^{6}$.

Further inspection of Figure 4 reveals that the eigenvalue errors are generally smaller in the case of $GL_{3}$ than those of $G_{3}$ and smaller for $G_{3}$ than for $G_{2}$. This observation confirms the theoretical analysis we present in Section 3 . The leading order terms in the Taylor expansions of the dispersion error (3.9) for $GL_{3}$, (3.6) for $G_{3}$, and (3.12) for $G_{2}$ have absolute values of 1/2880, 1/1440, and 1/720, respectively. Figure 4 shows the results when using $C^{2}$ cubic isogeometric analysis while Figure 5 shows these when using $C^{3}$ quartic and $C^{4}$ quintic isogeometric analysis. A slight extra-superconvergence is observed in Figure 5 for quartic and quintic isogeometric elements in coarse meshes. These numerical results confirm the eigenvalue error analysis we discuss in Section 4.1.

5.1.2 Eigenfunction errors

In this subsection, we verify numerically the analysis of Section 4.2. While Figures 2 and 2 show several sample plots of approximate and exact eigenfunctions, we show in this subsection convergence rates of the eigenfunction errors.

As the analysis predicts, the error in the eigenfunctions for the optimally-blended schemes does not exhibit extra orders of superconvergence. For the $H^{1}$-seminorm, all schemes yield a convergence of order $p$ for $p$-th order isogeometric elements with maximum continuity $C^{p-1}$ at element interfaces.

Figure 7 shows the $H^{1}$-seminorm errors of the eigenfunctions $u_{3}$ and $u_{11}$ for $C^{1}$ quadratic isogeometric analysis with quadrature rules $G_{3},GL_{3},G_{2}$, and $O_{2}$. The error convergence rate is two and the differences in the errors are fairly small. More precisely, Table 2 shows the errors for $u_{3}$. These numbers are the data for the left plot in Figure 7. Their differences are in the order of $10^{-6}$ for the case with mesh size $1/20$ and of a scale of $10^{-7}$ for the case with mesh size $1/40$. As a consequence, they have the same convergence rates. Similar results are observed for other eigenfunctions and details are omitted here.

For completeness, Figure 7 shows the eigenfunction errors in $L^{2}$-norm for $C^{1}$ quadratic case and Figure 8 for higher order cases. Using Lemma 3 and Theorem 3, we obtain an estimation for eigenfunction errors in $L^{2}$-norm $\|u_{j}-\tilde{u}_{j}^{h}\|_{0,\Omega}\leq Ch^{p+1}$ (see also (A.4) in the A). All the schemes yield errors which are of order $p+1$ for the $p$-th order isogeometric analysis of continuity $C^{p-1}$. These errors have the optimal convergence orders. We observe a superconvergent rate of eigenfunction error in $L^{2}$-norm for quintic elements on coarse meshes. Again, the differences in errors when different quadrature rules are utilized are small. These numerical results confirm the eigenfunction error analysis in Section 4.2.

5.2 Convergence study in two dimensions

The properties in two and higher dimensions are simple extensions of one dimensional cases when using tensor product meshes [39, 40]. Now, we consider the two-dimensional problem (2.2) with eigenvalues and eigenfunctions $\lambda_{jk}=\pi^{2}(j^{2}+k^{2}),u_{jk}=2\sin(j\pi x)\sin(k\pi y),j,k=1,2,\cdots.$ Again, we assume that once we solve the eigenvalue problem, the numerical eigenvalues $\lambda_{j}^{h}$ and $\tilde{\lambda}_{j}^{h}$ are sorted in ascending order and paired with the true eigenvalues $\lambda_{j}$.

We focus on studying the convergence behavior when using the optimal blending rules. Figure 9 shows the eigenvalue errors for $\lambda_{2,2}$ when using quadrature rules for $C^{1}$ quadratic and $C^{2}$ cubic isogeometric elements. The mesh configurations for quadratic and cubic cases are $4\times 4,8\times 8,16\times 16,32\times 32,$ and $64\times 64$. The eigenvalue errors of the optimal blending schemes are of order $2p+2$, while the standard quadratures yield order $2p$. These observations verify the analysis we discussed in Section 4.

5.3 An error estimator

We derive an error estimator in A and study the numerical performance. The effectivity index of the error estimator $R(\tilde{\lambda}_{j}^{h},\tilde{u}^{h}_{j})$ is defined in (A.9). Table 3 shows the effectivity indices of the estimations of the first and the fourth eigenfunction errors in energy-norm when using isogeomtric elements with $p=1,2,3$ in one dimension. The domain $\overline{\Omega}=[0,1]$ is discretized uniformly with $N=5,10,20,40$ elements. Table 3 shows that the effectivity indices are close to one and they converge to 1 at an almost quadratic rate, which confirms the discussion in A.

5.4 Computational efficiency

We now present the computational cost (time) while using the optimally-blended quadrature rules with a comparison to the case while using standard quadrature rules.

Table 4 shows the computational time of the matrix assembly using Gauss and optimal rules and the total simulation time for the eigenvalue problem (2.2) in two dimensions. Herein, the matrices are assembled using the tensor product structure with the corresponding one-dimensional matrices. The total time is dominated by the solution of the generalized matrix eigenvalue problem (2.16) which is the most time-consuming part of the simulation. As Table 4 shows, the total time required to solve the eigenvalue problem is similar for both Gauss and optimal rules. The cost of the assembly of the tensor product matrix does not increase much even when we evaluate two quadrature rules instead of one when using the optimally-blended rules. Moreover, nonstandard quadrature rules (with a minimal number of quadrature points), which are equivalent to the optimally-blended quadrature rules in the sense of producing the same mass and stiffness matrices, are developed in [41] and they reduce computational cost; see also the reduced rules in [42].

6 Concluding remarks

In this paper, we derive blending schemes to optimize the dispersion errors of isogeometric analysis. These optimally-blended quadrature rules lead to two extra order of convergence in the eigenvalue errors. Utilizing this superconvergence result, we derive an asymptotically-exact a posteriori error estimator. The optimal blending schemes are established for $p$ up to order seven for maximum continuity spaces. We will seek to generalize these results in the future work to arbitrary polynomial and variable continuity orders. Another future work would be the generalization of the optimal blending technique to numerical methods for the differential eigenvalue problem with interfaces.

7 Acknowledgments

This publication was made possible in part by the CSIRO Professorial Chair in Computational Geoscience at Curtin University and the Deep Earth Imaging Enterprise Future Science Platforms of the Commonwealth Scientific Industrial Research Organisation, CSIRO, of Australia. Additional support was provided by the European Union’s Horizon 2020 Research and Innovation Program of the Marie Skłodowska-Curie grant agreement No. 777778, the Mega-grant of the Russian Federation Government (N 14.Y26.31.0013), the Curtin Institute for Computation, and Institute for Geoscience Research (TIGeR). The authors thank Eric Chung’s constructive comments on the analysis.