# Optimal reduction of numerical dispersion for wave propagation problems. Part 1: Application to 1-D isogeometric elements

A Idesman

► To cite this version:

A Idesman. Optimal reduction of numerical dispersion for wave propagation problems. Part 1: Application to 1-D isogeometric elements. Computer Methods in Applied Mechanics and Engineering, 2017, 317, pp.970-992. 10.1016/j.cma.2017.01.014. hal-01495260

HAL Id: hal-01495260

https://hal.science/hal-01495260

Submitted on 24 Mar 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimal reduction of numerical dispersion for wave propagation problems. Part 1: Application to 1-D isogeometric elements

A. Idesman

Department of Mechanical Engineering, Texas Tech University, Lubbock, TX 79409-1021, USA

A numerical technique with the optimal coefficients of the stencil equation has been suggested. Based on this approach, new high-order isogeometric elements with the reduced dispersion error have been developed for wave propagation problems in the 1-D case. By the modification of the mass and stiffness matrices, the order of the dispersion error is improved from order $2p$ (the conventional elements) to order $4p$ (the new elements) where $p$ is the order of the polynomial approximations. It was shown that the new approach yields the maximum order of the dispersion error for the stencil equations related to the high-order isogeometric elements. The analysis of the dispersion error of the high-order isogeometric elements with the lumped mass matrix has also shown that independent of the procedures for the calculation of the lumped mass matrix, the second order of the dispersion error cannot be improved with the conventional stiffness matrix. However, the dispersion error for the lumped mass matrix can be improved from the second order to order $2p$ by the modification of the stiffness matrix. The numerical examples confirm the computational efficiency of the new high-order isogeometric elements with reduced dispersion. We have also shown that numerical results obtained by the new and conventional high-order isogeometric elements may include spurious oscillations due to the dispersion error. These oscillations can be quantified and filtered by the two-stage time-integration technique recently developed in our papers. The approach developed in the paper can be directly applied to other space-discretization techniques with similar stencil equations.

Keywords: Isogeometric elements; High-order elements with reduced dispersion; Wave propagation; Numerical dispersion

In the paper the reduction of the numerical dispersion error for wave propagation problems with the application to isogeometric elements in the 1-D case (Part 1) and the 2-D case (Part 2; see [1]) is considered. Wave propagation in an isotropic homogeneous medium is described by the following scalar wave equation in domain $\Omega$:

$$ \frac{\partial^2 u}{\partial t^2} - c^2 \nabla^2 u = 0, \qquad (1) $$

with the boundary conditions $n \cdot \nabla u = g_1$ on $\Gamma^t$ and $u = g_2$ on $\Gamma^u$, and the initial conditions $u(x, t = 0) = g_3$, $v(x, t = 0) = g_4$ in $\Omega$. Here, $u$ is the field variable, $v = \dot{u}$ is the velocity, $c$ is the wave velocity, $t$ is the time, $\Gamma^t$ and $\Gamma^u$ denote the natural and essential boundaries, $g_i$ ($i = 1, 2, 3, 4$) are the given functions, $n$ is the outward unit normal
on $\Gamma^t$. The application of the continuous Galerkin approach and the space discretization (e.g., the finite elements, spectral elements, isogeometric elements; see [2–4] and others) to Eq. (1) leads to a system of ordinary differential equations in time

$$ M \ddot{U} + c^2 K U = 0, \qquad (2) $$

with

$$ M = \sum_{e} M^{e}, \quad K = \sum_{e} K^{e}, \qquad (3) $$

where $U(t)$ is the vector of the field variable, the global mass $M$ and stiffness $K$ matrices have a banded structure and are obtained by the summation of the corresponding local (element $\Omega^e$) matrices $M^e$ and $K^e$:

$$ M^e = \int_{\Omega^e} N^T N d\Omega^e, \qquad (4) $$

$$ K^e = \int_{\Omega^e} \left[ \frac{\partial N}{\partial \mathbf{x}} \right]^T \left[ \frac{\partial N}{\partial \mathbf{x}} \right] d\Omega^e. \qquad (5) $$

Here, $N$ and $\frac{\partial N}{\partial x}$ are the shape matrix and its derivative with respect to the physical coordinate $\mathbf{x}$; see [2,3,5]. Due to the space discretization, the exact solution to Eq. (2) contains the numerical dispersion error. Usually the analysis of the numerical dispersion error and its improvement for many space-discretization techniques such as the finite elements, spectral elements, isogeometric elements and others start with the analysis and modifications of the elemental mass and stiffness matrices; see [3,4,6–25]. For example, one simple and effective finite-element technique for acoustic and elastic wave propagation problems is based on the calculation of the mass matrix $M$ in Eq. (2) as a weighted average of the consistent and lumped mass matrices; see [13–16] and others. For the 1-D case and the linear finite elements, this approach reduces the error in the wave velocity for harmonic waves from the second order to the fourth order of accuracy. However, for harmonic wave propagation in the 2-D and 3-D cases, these results are not valid (nevertheless, in the multi-dimensional case, the averaged mass matrix yields more accurate results compared with the standard mass matrix).

Some useful techniques for the reduction of the numerical dispersion of the finite element solutions for the Helmholtz equation and for the time-harmonic elastic waves have been developed for a given wave number; e.g., see [26,27]. In these approaches, the parameters that improve the numerical dispersion depend on the wave number. Therefore, for general wave propagation problems in the time domain that include many other wave numbers, these approaches may not be very efficient.

An interesting technique with implicit and explicit time-integration methods is suggested in [17] for acoustic waves. It is based on the modified integration rule for the calculation of the mass and stiffness matrices for the linear finite elements. In contrast to the averaged mass matrix, the use of the modified integration rule increases the accuracy for the phase velocity from the second order to the fourth order in the general multi-dimensional case of acoustic waves. A similar improvement in the order of the numerical dispersion error for the linear elements in the 1-D case yields the selective mass scaling technique developed in [28]. Using the mimetic finite difference approach, similar results for acoustic waves and electric waves in the 2-D case have been obtained in [12,29] for the linear finite elements with reduced dispersion and reduced numerical anisotropy.

The dispersion reduction technique for the high-order finite elements have been suggested in [7] for acoustic waves. This technique is based on the calculation of the mass matrix $M$ in Eq. (2) as a weighted average of the consistent and lumped mass matrices for the high-order finite elements. It was also shown in [7] that the same results can be obtained with the modified integration rule for the mass matrix. With this technique, the dispersion error is improved from the order $2p$ to the order $2p + 2$ ($p$ is the order of polynomial approximations). However, even for the quadratic elements, the modified stencil equation derived in [7] in the 1-D case is based on the replacement of the mid-side degree of freedom within the 1-D element by the degrees of freedom related to the nodes located at the ends of the same element (it seems that such a replacement is impossible in the 2-D or 3-D cases due to the belonging of the same mid-side degree of freedom to different elements).

New isogeometric elements with a higher order of continuity across elements are suggested in [3,5,30–32] for dynamics problems. It has been shown in these papers that the isogeometric elements yield more accurate numerical
results for wave propagation problems compared with the high-order finite elements. The modification of the non-diagonal mass matrix in [24,25] for the isogeometric elements allows the increase in the order of the dispersion error from order $2p$ to order $2p + 2$ in the 1-D case and for one specific direction of harmonic waves in the 2-D case (these techniques do not improve the order of the dispersion error in the general 2-D case). It is also necessary to mention that in contrast to the high order finite and spectral elements, the stencil equation for the isogeometric elements on uniform meshes is the same for all internal degrees of freedom (for the high order finite and spectral elements there are several stencil equations with different structures depending on the location of nodes). This simplifies the analysis of the numerical dispersion error for the isogeometric elements. However, along with the advantages, there are some issues with the isogeometric elements. For example, the lumped mass matrix for the isogeometric elements decreases the order of the dispersion error to the second order of accuracy; e.g., see [30]. Therefore, a special iterative procedure for linear dynamics problems was developed in [30] in order to obtain the same order of the dispersion error for the lumped and consistent mass matrices.

Despite the improved accuracy of the above-mentioned techniques, all of them may yield spurious high-frequency oscillations at high-frequency and impact loadings. These spurious oscillations may lead to divergent results at mesh refinement; e.g., see [33,34]. The existing techniques for the suppression of the spurious oscillations are based on the introduction of numerical dissipation (or artificial viscosity) for all time increments at the integration of Eq. (2); e.g., see [2,35] and many others. However, the existing techniques do not quantify the range of the spurious frequencies and also may yield inaccurate results. Recently a special approach for the accurate solution of acoustic and elastic wave propagation problems has been developed in [22,23,33,36–40]. This approach includes the two-stage time-integration techniques with the quantification and the filtering of spurious oscillations for different space-discretization methods, an a priori global error estimator in time, the linear finite elements with reduced dispersion. The numerical results in [22,23] showed that the reduction of the dispersion error for the linear finite elements leads to a huge reduction in the computation time.

The derivation of the analytical formulas for the numerical dispersion error is difficult and very often intractable due to very long and complicated non-linear expressions with a large number of different parameters. This is an issue for high-order space-discretization methods in the multidimensional case. For example, in the recent paper [41] the numerical dispersion analysis of the conventional isogeometric elements for elastic waves in the 2-D case was presented numerically (without analytical expressions). To simplify the analytical study and improvement of the numerical dispersion error for high-order space-discretization methods, this paper starts the analysis with the corresponding stencil equation with unknown coefficients. It is important to note that all these unknown coefficients can be found by the minimization of the numerical dispersion error. Therefore, the accuracy of the approach considered in the paper cannot be improved with other coefficients of the stencil equation and is optimal. With the optimal coefficients the order of the numerical dispersion error is increased from order $2p$ to order $4p$. In Sections 1.1 and 1.3 we will consider the dispersion analysis of the conventional high-order isogeometric elements in the 1-D case. Based on the dispersion analysis, in Sections 1.2 and 1.4 we will develop new quadratic and cubic isogeometric elements with reduced dispersion in the 1-D case. We will also show that by the modification of the stiffness matrix, the dispersion error of the high-order isogeometric elements with the lumped mass matrix can be improved; see Sections 1.2.2 and 1.4.1. The numerical examples in Section 2 will show the computational efficiency of the new approach. For the derivation of many analytical expressions presented below we have used the computational program “Mathematica”.

# 1. Dispersion analysis in the 1-D case

Inserting time-harmonic solutions

$$u(x, t) = \exp(i \omega t) u(x) \qquad (6)$$

into the wave equation (1) in the 1-D case ($\nabla^2 u = \frac{\partial^2 u}{\partial x^2}$) leads to its reduction to the Helmholtz equation

$$\frac{\partial^2 u}{\partial x^2} + k^2 u = 0, \qquad (7)$$

where $\omega$ is the angular velocity, $k = \omega/c$ is the wave number, $i = \sqrt{-1}$. Eq. (7) allows the following exact solutions:

$$u(x) = \exp(\pm ikx). \qquad (8)$$
<figure><img src="image_2.png" /><figcaption>Fig. 1. The spatial locations of the degrees of freedom $u_i$ ($i = A-2, A-1, A, A+1, A+2$) contributing to the stencil equation for the degree of freedom $u_A$ of the new and conventional quadratic isogeometric elements. The three elements shown are used for the calculation of the stencil equation for the degree of freedom $u_A$.</figcaption></figure>

After the space discretization, Eq. (7) reduces to

$$ (K - k^2 M) U = 0. \qquad (9) $$

For the dispersion analysis presented below we use uniform meshes. For the calculation of the mass and stiffness matrices of the conventional high-order isogeometric elements we use univariate B-spline basis functions defined recursively starting with $p = 0$

$$ N_{i,p}(\xi) = \begin{cases} 1, & \text{if } \xi_i \le \xi < \xi_{i+1} \\ 0, & \text{otherwise.} \end{cases} \qquad (10) $$

For $p \ge 1$:

$$ N_{i,p}(\xi) = \frac{\xi - \xi_i}{\xi_{i+p} - \xi_i} N_{i,p-1}(\xi) + \frac{\xi_{i+p+1} - \xi}{\xi_{i+p+1} - \xi_{i+1}} N_{i+1,p-1}(\xi), \qquad (11) $$

where a knot vector {$\xi_1 = 0, \dots, \xi_{n+p+1} = 1$} is a set of non-decreasing real numbers representing coordinates in the parametric space of the curve, $p$ is the order of the B-spline, $n$ is the number of basis functions.

### 1.1. Dispersion analysis for the conventional quadratic isogeometric elements

Let us first start with the dispersion analysis of the conventional quadratic isogeometric elements (the detailed dispersion analysis for the conventional quadratic isogeometric elements can be also found in [3,5]). The elemental mass $M^e$ and stiffness $K^e$ matrices of a typical quadratic interior element are

$$ M^e = \frac{h}{120} \begin{pmatrix} 6 & 13 & 1 \\ 13 & 54 & 13 \\ 1 & 13 & 6 \end{pmatrix}, \qquad K^e = \frac{1}{6h} \begin{pmatrix} 2 & -1 & -1 \\ -1 & 2 & -1 \\ -1 & -1 & 2 \end{pmatrix}, \qquad (12) $$

where $h = 1/n_{el} = 1/(N - p)$, $n_{el}$ is the number of elements, $N$ is the number of control points, and $p = 2$ is the order of the basis functions for the quadratic elements. Considering the discretized Helmholtz equation (9) on an infinite line with the sequence of $h$-spaced control points $x_A = hA$, the stencil equation can be calculated with the help of Eqs. (3) and (12) for three neighboring elements and has the following form (see Fig. 1 and [3]):

$$ k^2 h^2 (u_{A-2} + 26u_{A-1} + 66u_A + 26u_{A+1} + u_{A+2}) + 20[u_{A-2} + 2(u_{A-1} - 3u_A + u_{A+1}) + u_{A+2}] = 0. \qquad (13) $$

This equation allows the following solutions (similar to Eq. (8)):

$$ u_A = \exp(\pm ik_h hA), \qquad (14) $$

where $k_h$ is the numerical wave number.

**Remark.** In the dispersion analysis we do not include the solutions of Eq. (13) related to the spurious evanescent waves (see [3]). The numerical examples in [3] show that these solutions do not affect the numerical results.

Inserting Eq. (14) into Eq. (13) we can find the following relation between the exact and numerical wave numbers $k$ and $k_h$:

$$ \frac{k}{k_h} = \frac{1}{(k_h h)} \frac{2\sqrt{5}\sqrt{-2 \cos(k_h h) - \cos(2(k_h h)) + 3}}{\sqrt{26 \cos(k_h h) + \cos(2(k_h h)) + 33}}. \qquad (15) $$
<figure><img src="image_3.png" /><figcaption>Fig. 2. The numerical dispersion errors $e = \frac{k_h-k}{k}$ (a) and $\ln(|e|)$ (b) of the quadratic isogeometric elements as a function of the mesh size $k_h h$ (a) and $\ln(k_h h)$ (b). Curve 1 corresponds to the conventional quadratic isogeometric elements with the non-diagonal mass matrix. Curves 2 and 3 correspond to the new quadratic isogeometric elements with the 8th and 6th order of the numerical dispersion error (the non-diagonal mass matrix). The slopes of the curves in (b) show the order of convergence at the small mesh size $h$.</figcaption></figure>

The numerical dispersion error $e = \frac{k_h-k}{k} = \frac{k_h}{k} - 1$ calculated with the help of Eq. (15) is plotted in Fig. 2(a); see curve 1. The same results are also plotted in Fig. 2(b) in the logarithmic scale in order to show the order of convergence at mesh refinement. As can be seen, the dispersion error decreases with the decrease in the mesh size $k_h h$. Expanding the right-hand side of Eq. (15) into a Taylor series at small $h \ll 1$, we get

$$ \frac{k}{k_h} = 1 + \frac{(k_h h)^4}{1440} + \frac{(k_h h)^6}{6720} + \frac{23(k_h h)^8}{4147200} + O((k_h h)^{10}), \quad (16) $$

i.e., the conventional quadratic isogeometric elements yield the fourth order of the dispersion error for $(k/k_h - 1)$; see also Fig. 2(b).

## 1.2. A new approach for the quadratic isogeometric elements with reduced dispersion

For the derivation of new quadratic isogeometric elements with reduced dispersion let us start with the stencil equation with unknown coefficients for the degree of freedom $u_A$ that can be written as follows:

$$ k^2 h^2 [m_1(u_{A-2} + u_{A+2}) + m_2(u_{A-1} + u_{A+1}) + m_3 u_A] - [k_1(u_{A-2} + u_{A+2}) + k_2(u_{A-1} + u_{A+1}) + k_3 u_A] = 0, \quad (17) $$

where unknown coefficients $m_j$ and $k_j$ ($j = 1, 2, 3$) are calculated in terms of the elemental mass and stiffness matrices, respectively (at this point, these matrices are not defined). For example, comparing Eqs. (17) and (13), we can find that $m_1 = 1, m_2 = 26, m_3 = 66, k_1 = -20, k_2 = -40, k_3 = 120$ for the conventional quadratic isogeometric elements. Similar to the conventional isogeometric elements we assumed the symmetry for coefficients $m_j$ and $k_j$ in the stencil equation (17) for the degrees of freedom $u_j$ ($j = A - 1, A + 1, A - 2, A + 2$) symmetrically located with respect to the degree of freedom $u_A$. The idea of the new approach is to find the optimal coefficients of the stencil equation (17) that reduce the order of the dispersion error. Inserting Eq. (14) into Eq. (17) we can find the following relation between the exact and numerical wave numbers $k$ and $k_h$:

$$ \frac{k}{k_h} = \frac{1}{(k_h h)} \frac{\sqrt{2k_1 \cos(2(k_h h)) + 2k_2 \cos(k_h h) + k_3}}{\sqrt{2m_1 \cos(2(k_h h)) + 2m_2 \cos(k_h h) + m_3}}. \quad (18) $$

Expanding the right-hand side of Eq. (18) into a Taylor series at small $h \ll 1$ and analyzing the first terms of a Taylor series we can find that the following constraint $k_3 = -2k_1 - 2k_2$ should be met in order to get at least the first order of the dispersion error. In this case Eq. (18) can be rewritten as:

$$ \frac{k}{k_h} = \frac{1}{(k_h h)} \frac{\sqrt{2\sqrt{k_1 \cos(2(k_h h)) + k_2 \cos(k_h h) - k_1 - k_2}}}{\sqrt{2m_1 \cos(2(k_h h)) + 2m_2 \cos(k_h h) + m_3}}. \quad (19) $$
Expanding again the right-hand side of Eq. (19) into a Taylor series at small $h \ll 1$ and equating the first term to one and the following three terms to zero, we can find the coefficients $m_j$ and $k_j$ ($j = 1, 2, 3$) from a system of four algebraic equations. They are

$$m_1 = \frac{23a_1}{2358}, m_2 = \frac{344a_1}{1179}, m_3 = a_1, k_1 = -\frac{155a_1}{786}, k_2 = -\frac{320a_1}{393}, k_3 = \frac{265a_1}{131}, \quad (20)$$

where $a_1$ is the arbitrary coefficient. The numerical dispersion error $e = \frac{k_h-k}{k}$ calculated with the help of Eqs. (19) and (20) is plotted in Fig. 2; see curve 2. Expanding the right-hand side of Eq. (19) into a Taylor series at small $h \ll 1$ with the coefficients given by Eq. (20) we get

$$\frac{k}{k_h} = 1 - \frac{79(k_h h)^8}{9525600} - \frac{2633(k_h h)^{10}}{2640496320} + O((k_h h)^{12}), \quad (21)$$

i.e., the stencil equation (17) with the coefficients given by Eq. (20) yields the 8th order of the dispersion error; see also curve 2 in Fig. 2(b).

**Remark.** The new approach also improves the order of the local truncation error in space by four orders (similar to the improvement of the order of the numerical dispersion error for the new quadratic isogeometric elements); see Appendix.

Let us find the elemental mass and stiffness matrices that yield the stencil equation (17) with the coefficients $m_j$ and $k_j$ ($j = 1, 2, 3$) given by Eq. (20). Similar to the conventional mass and stiffness matrices, let us assume the following form of the new mass and stiffness matrices:

$$M^e = h \begin{pmatrix} m_{11} & m_{12} & m_{13} \\ m_{12} & m_{22} & m_{12} \\ m_{13} & m_{12} & m_{11} \end{pmatrix}, \quad K^e = \frac{1}{h} \begin{pmatrix} k_{11} & k_{12} & k_{13} \\ k_{12} & k_{22} & k_{12} \\ k_{13} & k_{12} & k_{11} \end{pmatrix}, \quad (22)$$

where we take into account the symmetry of the coefficients of the mass and stiffness matrices for the degrees of freedom contributing to these matrices. In this case the matrices $M^e$ and $K^e$ depend on 8 unknown coefficients $m_{11}, m_{12}, m_{13}, m_{22}$ and $k_{11}, k_{12}, k_{13}, k_{22}$. With the help of (3), the coefficients of the stencil equation (17) can be expressed in terms of the coefficients of the matrices $M^e$ and $K^e$ as follows:

$$m_1 = m_{13}, m_2 = 2m_{12}, m_3 = 2m_{11} + m_{22}, k_1 = k_{13}, k_2 = 2k_{12}, k_3 = 2k_{11} + k_{22}. \quad (23)$$

Solving simultaneously Eqs. (20) and (23) we can find the mass and stiffness matrices for the quadratic isogeometric elements with reduced dispersion. They have the following form:

$$M^e = \frac{h}{120} \begin{pmatrix} a_2 & \frac{6880a_1}{393} & \frac{460a_1}{393} \\ \frac{6880a_1}{393} & 120a_1 - 2a_2 & \frac{6880a_1}{393} \\ \frac{460a_1}{393} & \frac{6880a_1}{393} & a_2 \end{pmatrix}, \quad K^e = \frac{1}{6h} \begin{pmatrix} a_3 & -\frac{320a_1}{131} & -\frac{155a_1}{131} \\ -\frac{320a_1}{131} & \frac{1590a_1}{131} - 2a_3 & -\frac{320a_1}{131} \\ -\frac{155a_1}{131} & -\frac{320a_1}{131} & a_3 \end{pmatrix}, \quad (24)$$

where $a_j$ ($j = 1, 2, 3$) are three arbitrary coefficients. The dispersion error is independent of the values of these coefficients $a_j$. Some coefficients $a_j$ can be found from the following conditions. Let us assume that the velocity is the same for the entire domain. In this case we can calculate the kinetic energy for one element with the help of the mass matrix and analytically for the considered element. Equating these two expressions we get

$$V_0^T M^e V_0 = v_0^2 h, \quad (25)$$

where $V_0 = (v_0, v_1, v_2)^T$ is the vector with three equal components. Calculating the left-hand side of Eq. (25) with the help of Eq. (24) we get that $a_1 = 131/210$. Next, let us assume that the displacement is the same for the entire domain. In this case strains, stresses and forces are zero. The forces for one element can be calculated with the help of the stiffness matrix:

$$K^e U_0 = 0, \quad (26)$$
where $U_0 = (u_0, u_0, u_0)^T$ is the vector with three equal components. Calculating the left-hand side of Eq. (26) with the help of Eq. (24) we get that $a_3 = 95/42$. For $a_1 = 131/210$ and $a_3 = 95/42$, Eq. (24) can be rewritten as follows:

$$ M^e = \frac{h}{120} \begin{pmatrix} a_2 & \frac{688}{63} & \frac{46}{63} \\ \frac{688}{63} & \frac{524}{7} - 2a_2 & \frac{688}{63} \\ \frac{46}{63} & \frac{688}{63} & a_2 \end{pmatrix}, \qquad K^e = \frac{1}{6h} \begin{pmatrix} \frac{95}{42} & \frac{32}{21} & \frac{31}{42} \\ \frac{32}{21} & \frac{64}{21} & \frac{32}{21} \\ \frac{31}{42} & \frac{32}{21} & \frac{95}{42} \end{pmatrix}. \qquad (27) $$

**Remark.** All unknown coefficients of the stencil equation (17) have been found from the minimization of the order of the dispersion error; see Eq. (20). The presence of the arbitrary coefficient $a_1$ in Eq. (20) is explained by the form of the homogeneous discretized Helmholtz equation (9); i.e., the multiplication of the mass and stiffness matrices by any coefficient does not change this equation and the stencil equation (17). The unknown coefficient $a_1$ can be found from Eq. (25). This leads to the following constraint for the coefficients $m_i$ of the stencil equation (17):

$$ 2m_1 + 2m_2 + m_3 = 1. \qquad (28) $$

In this case it follows from Eqs. (20) and (28) that $a_1 = 131/210$. Because all coefficients of the stencil equation can be found from the analysis of the dispersion error, the 8th order of the dispersion error in Eq. (21) is maximum possible for the considered form of the stencil equation (that is related to the support of basis functions by a specific number of elements) for all quadratic isogeometric elements.

### 1.2.1. The 6th order of accuracy with the modified mass matrix and the conventional stiffness matrix

It is interesting to note that in our papers [22,23] related to the linear finite elements with reduced dispersion, in the 1-D case we have modified just the mass matrix for the reduction of the dispersion error. Below we will show that it is possible to reduce the order of the dispersion error for the isogeometric elements with the conventional stiffness matrix and the modified mass matrix (partly similar to [24]). In order to find the corresponding modifications of the mass matrix, let us consider the stencil equation (17) with $k_1 = -20$, $k_2 = -40$, $k_3 = 120$; see Eq. (13). Similar to the previous section, with the help of Eqs. (14) and (17) we find the dispersion equation, expand it into a Taylor series and equate the first term of a Taylor series to unity, and the following two terms to zero. From these three algebraic equations we can find that $m_1 = 7/720$, $m_2 = 19/90$, $m_3 = 67/120$. The numerical dispersion error $e = \frac{k_h-k}{k}$ in this case is plotted in Fig. 2; see curve 3. A Taylor series of the dispersion equation with these coefficients is

$$ \frac{k}{k_h} = 1 + \frac{11(k_h h)^6}{120960} + \frac{(k_h h)^8}{345600} + O\left((k_h h)^{10}\right), \qquad (29) $$

i.e., we can get the 6th order of the dispersion error in this case (see also Fig. 2(b)). Finally, using the first three equalities in Eq. (23), we can find the following form of the modified mass matrix:

$$ M = \frac{h}{120} \begin{pmatrix} \frac{67 - a_1}{2} & \frac{38}{3} & \frac{7}{6} \\ \frac{38}{3} & a_1 & \frac{38}{3} \\ \frac{7}{6} & \frac{38}{3} & \frac{67 - a_1}{2} \end{pmatrix} \qquad (30) $$

with one arbitrary coefficient $a_1$. The dispersion error is independent of the value of this coefficient.

**Remark.** We should note that at $a_1 = 54\frac{2}{3}$ the mass matrix in Eq. (30) coincides with that in [24]. The coefficients $m_1 = 7/720$, $m_2 = 19/90$, $m_3 = 67/120$ of the stencil equation are independent of $a_1$ and also coincide with those in [24]. This means that curves 3 in Fig. 2 correspond to the dispersion error of the isogeometric elements developed in [24]. In contrast to [24], we use a completely different derivation and also show that the 6th order of the dispersion error in [24] for the quadratic isogeometric elements with the conventional stiffness matrix is maximum possible (because the above-mentioned coefficients $m_i$ are optimal).
<figure><img src="image_4.png" /><figcaption>Fig. 3. The numerical dispersion errors $e = \frac{k_h-k}{k}$ (a) and $\ln(|e|)$ (b) of the quadratic and cubic isogeometric elements with the lumped mass matrix as a function of the mesh size $k_h h$ (a) and $\ln(k_h h)$ (b). Curves 1 and 2 correspond to the conventional quadratic and cubic isogeometric elements with the lumped mass matrix. Curves 3 and 4 correspond to the new quadratic and cubic isogeometric elements with the lumped mass matrix. The slopes of the curves in (b) show the order of convergence at the small mesh size $h$.</figcaption></figure>

### 1.2.2. Numerical dispersion of the quadratic isogeometric elements with the lumped mass matrix

The approach presented in the previous sections can be also used for the analysis and the reduction of the dispersion error for the lumped mass matrix. In this case we should consider the stencil equation (17) with $m_1 = m_2 = 0$. For the conventional stiffness matrix with $k_1 = 20, k_2 = 40, k_3 = -60$, the coefficient $m_3$ can be found from the analysis of the dispersion equation (as described in the previous sections) and is equal to $m_3 = 1$. The numerical dispersion error $e = \frac{k_h-k}{k}$ for the lumped mass matrix and the conventional stiffness matrix is plotted in Fig. 3(a); see curve 1. The same results are also plotted in Fig. 3(b) in the logarithmic scale in order to show the order of convergence at mesh refinement. A Taylor series of the dispersion equation has the following form

$$ \frac{k}{k_h} = 1 - \frac{(k_h h)^2}{8} + \frac{43(k_h h)^4}{5760} + O\left((k_h h)^6\right), \qquad (31) $$

i.e., the dispersion error has the second order of accuracy for the lumped mass matrix (see also Fig. 3(b)). Using the first three equalities in Eq. (23), we can find the following form of the lumped mass matrix:

$$ M^e = \frac{h}{120} \begin{pmatrix} \frac{120 - a_1}{2} & 0 & 0 \\ 0 & a_1 & 0 \\ 0 & 0 & \frac{120 - a_1}{2} \end{pmatrix} \qquad (32) $$

with one arbitrary coefficient $a_1$. The dispersion error is independent of the value of this coefficient. It is important to note that according to the analysis of the stencil equation with the lumped mass matrix, there are no additional degrees of freedom that allow the reduction of the dispersion error with the conventional stiffness matrix. Independent of the techniques for the calculations of the lumped mass matrix, the best accuracy corresponds to the second order of the dispersion error.

However, the accuracy of the results with the lumped mass matrix can be improved by the modification of the stiffness matrix. Let us consider again the stencil equation (17) with $m_1 = m_2 = 0$ and undefined coefficients $k_j$ ($j = 1, 2, 3$). Using the approach similar to that in Section 1.2 and equating the first term of a Taylor series of the dispersion equation to unity, and the following one term to zero, we can find that $m_3 = 3/4a_1, k_1 = a_1/16, k_2 = -a_1, k_3 = -2k_1 - 2k_2 = 15/8a_1$ with one arbitrary coefficient $a_1$. The numerical dispersion error $e = \frac{k_h-k}{k}$ for the lumped mass matrix and the modified stiffness matrix is plotted in Fig. 3; see curve 3. In this case, a Taylor series of the dispersion equation has the following form

$$ \frac{k}{k_h} = 1 - \frac{(k_h h)^4}{180} + \frac{(k_h h)^6}{2016} + O\left((k_h h)^8\right), \qquad (33) $$
<figure><img src="image_5.png" /><figcaption>Fig. 4. The spatial locations of the degrees of freedom $u_i$ ($i = A-3, A-2, A-1, A, A+1, A+2, A+3$) contributing to the stencil equation for the degree of freedom $u_A$ of the new and conventional cubic isogeometric elements. The four elements shown are used for the calculation of the stencil equation for the degree of freedom $u_A$.</figcaption></figure>

i.e., we can get the fourth order for the dispersion error (see also Fig. 3(b)). Using Eq. (23), we can find the following form of the lumped mass matrix and the modified stiffness matrix:

$$ M^e = \frac{h}{120} \begin{pmatrix} \frac{1}{2}(15a_1 - a_2) & 0 & 0 \\ 0 & a_2 & 0 \\ 0 & 0 & \frac{1}{2}(15a_1 - a_2) \end{pmatrix}, \quad K^e = \frac{1}{6h} \begin{pmatrix} \frac{1}{16}(15a_1 - 8a_3) & -\frac{a_1}{2} & \frac{a_1}{16} \\ -\frac{a_1}{2} & a_3 & -\frac{a_1}{2} \\ \frac{a_1}{16} & -\frac{a_1}{2} & \frac{1}{16}(15a_1 - 8a_3) \end{pmatrix} \quad (34) $$

with three arbitrary coefficients $a_j$ ($j = 1, 2, 3$) that do not affect the dispersion error. Some coefficients can be found with the help of Eqs. (25) and (26) and they are $a_1 = a_3 = 8$. With these coefficients, Eq. (34) can be rewritten as:

$$ M^e = \frac{h}{120} \begin{pmatrix} 60 - \frac{a_2}{2} & 0 & 0 \\ 0 & a_2 & 0 \\ 0 & 0 & 60 - \frac{a_2}{2} \end{pmatrix}, \quad K^e = \frac{1}{6h} \begin{pmatrix} \frac{7}{2} & -4 & \frac{1}{2} \\ -4 & 8 & -4 \\ \frac{1}{2} & -4 & \frac{7}{2} \end{pmatrix}. \quad (35) $$

### 1.3. Dispersion analysis for the conventional cubic isogeometric elements

The analysis of the dispersion error for the conventional cubic isogeometric elements with $p = 3$ can be done similar to Section 1.1. The mass and stiffness matrices of a typical interior cubic element are

$$ M^e = \frac{h}{5040} \begin{pmatrix} 20 & 129 & 60 & 1 \\ 129 & 1188 & 933 & 60 \\ 60 & 933 & 1188 & 129 \\ 1 & 60 & 129 & 20 \end{pmatrix}, \quad K^e = \frac{1}{120h} \begin{pmatrix} 6 & 7 & -12 & -1 \\ 7 & 34 & -29 & -12 \\ -12 & -29 & 34 & 7 \\ -1 & -12 & 7 & 6 \end{pmatrix}. \quad (36) $$

These matrices lead to the following stencil equation calculated with the help of four elements; see Fig. 4:

$$ k^2 h^2 (u_{A-3} + 120u_{A-2} + 1191u_{A-1} + 2416u_A + 1191u_{A+1} + 120u_{A+2} + u_{A+3}) \\ + 42(u_{A-3} + 24u_{A-2} + 15u_{A-1} - 80u_A + 15u_{A+1} + 24u_{A+2} + u_{A+3}). \quad (37) $$

Inserting Eq. (14) into Eq. (37) we can find the following relation between the exact and numerical wave numbers $k$ and $k_h$:

$$ \frac{k}{k_h} = \frac{1}{(k_h h)} \frac{\sqrt{42\sqrt{-15 \cos(k_h h) - 24 \cos(2(k_h h)) - \cos(3(k_h h)) + 40}}}{\sqrt{1191 \cos(k_h h) + 120 \cos(2(k_h h)) + \cos(3(k_h h)) + 1208}}. \quad (38) $$

The numerical dispersion error $e = \frac{k_h - k}{k} = \frac{k_h}{k} - 1$ calculated with the help of Eq. (38) is plotted in Fig. 5(a); see curve 1. The same results are also plotted in Fig. 5(b) in the logarithmic scale in order to show the order of convergence at mesh refinement. As can be seen, the dispersion error decreases with the decrease in the mesh size $k_h h$. Expanding the right-hand side of Eq. (38) into a Taylor series at small $h \ll 1$, we get

$$ \frac{k}{k_h} = 1 + \frac{(k_h h)^6}{60480} + \frac{(k_h h)^8}{120960} + \frac{(k_h h)^{10}}{1064448} + \frac{16519(k_h h)^{12}}{475517952000} + O\left((k_h h)^{14}\right), \quad (39) $$
<figure><img src="image_6.png" /><figcaption>Fig. 5. The numerical dispersion errors $e = \frac{k_h-k}{k}$ (a) and $\ln(|e|)$ (b) of the cubic isogeometric elements as a function of the mesh size $k_h h$ (a) and $\ln(k_h h)$ (b). Curves 1 and 2 correspond to the conventional and new cubic isogeometric elements with the non-diagonal mass matrix. The slopes of the curves in (b) show the order of convergence at the small mesh size $h$.</figcaption></figure>

i.e., the conventional cubic isogeometric elements yield the 6th order of the dispersion error for $(k/k_h-1)$; see also Fig. 5(b).

## 1.4. A new approach for the cubic isogeometric elements with reduced dispersion

The derivation of new cubic isogeometric elements with reduced dispersion is similar to that for the quadratic elements in Section 1.2. Let us start with the stencil equation with unknown coefficients for the degree of freedom $u_A$ that can be written as follows:

$$k^2 h^2 [m_2(u_{A-2} + u_{A+2}) + m_3(u_{A-1} + u_{A+1}) + m_4 u_A + m_1(u_{A-3} + u_{A+3})] - [k_2(u_{A-2} + u_{A+2}) + k_3(u_{A-1} + u_{A+1}) + k_4 u_A + k_1(u_{A-3} + u_{A+3})] = 0, \quad (40)$$

where the unknown coefficients $m_j$ and $k_j$ ($j = 1, 2, 3, 4$) are calculated in terms of the elemental mass and stiffness matrices, respectively (at this point, these matrices are not defined). For example, comparing Eqs. (40) and (37), we can find that $m_1 = 1$, $m_2 = 120$, $m_3 = 1191$, $m_4 = 2416$, $k_1 = -42$, $k_2 = -1008$, $k_3 = -630$, $k_4 = 3360$ for the conventional cubic isogeometric elements. Similar to the conventional isogeometric elements we assumed the symmetry for the coefficients $m_j$ and $k_j$ in the stencil equation (40) for the degrees of freedom $u_j$ ($j = A-1, A+1, A-2, A+2, A-3, A+3$) symmetrically located with respect to the degree of freedom $u_A$. Inserting Eq. (14) into Eq. (40) we can find the following relation between the exact and numerical wave numbers $k$ and $k_h$:

$$\frac{k}{k_h} = \frac{1}{(k_h h)} \frac{\sqrt{2k_2 \cos(2(k_h h)) + 2k_3 \cos(k_h h) + k_4 + 2k_1 \cos(3(k_h h))}}{\sqrt{2m_2 \cos(2(k_h h)) + 2m_3 \cos(k_h h) + m_4 + 2m_1 \cos(3(k_h h))}} \quad (41)$$

Expanding the right-hand side of Eq. (41) into a Taylor series at small $h \ll 1$ and analyzing the first terms of a Taylor series we can find that the following constraint $k_4 = -2k_1 - 2k_2 - 2k_3$ should be met in order to get at least the first order of the dispersion error. Expanding the right-hand side of Eq. (41) with $k_4 = -2k_1 - 2k_2 - 2k_3$ into a Taylor series at small $h \ll 1$ and equating the first term to unity and the following five terms to zero, we can find the coefficients $m_j$ and $k_j$ ($j = 1, 2, 3, 4$) from a system of six algebraic equations. They are

$$m_1 = a_1, \quad m_2 = 36774a_1/619, \quad m_3 = 329913a_1/619, \quad m_4 = 725308a_1/619, \quad (42)$$

$$k_1 = -49483a_1/1857, \quad k_2 = -261954a_1/619, \quad k_3 = -263655a_1/619, \quad k_4 = 3252620a_1/1857, \quad (43)$$

where $a_1$ is the arbitrary coefficient. The numerical dispersion error $e = \frac{k_h-k}{k}$ calculated with the help of Eqs. (41)-(43) is plotted in Fig. 5; see curve 2. Expanding the right-hand side of Eq. (41) into a Taylor series at small $h \ll 1$ with the coefficients given by Eqs. (42) and (43) we get

$$\frac{k}{k_h} = 1 - \frac{12741(k_h h)^{12}}{389701312000} + O((k_h h)^{14}). \quad (44)$$
Eq. (44) corresponds to the 12th order of the dispersion error; see also Fig. 5(b).

**Remark.** The new approach also improves the order of the local truncation error in space by six orders (similar to the improvement of the order of the numerical dispersion error for the new cubic isogeometric elements); see Appendix.

Let us find the elemental mass and stiffness matrices that yield the stencil equation (40) with the coefficients $m_j$ and $k_j$ ($j = 1, 2, 3, 4$) given by Eqs. (42) and (43). Similar to the conventional mass and stiffness matrices, let us assume the following form of the new mass and stiffness matrices:

$$ M^e = h \begin{pmatrix} m_{11} & m_{12} & m_{13} & m_{14} \\ m_{12} & m_{22} & m_{23} & m_{13} \\ m_{13} & m_{23} & m_{22} & m_{12} \\ m_{14} & m_{13} & m_{12} & m_{11} \end{pmatrix}, \qquad K^e = \frac{1}{h} \begin{pmatrix} k_{11} & k_{12} & k_{13} & k_{14} \\ k_{12} & k_{22} & k_{23} & k_{13} \\ k_{13} & k_{23} & k_{22} & k_{12} \\ k_{14} & k_{13} & k_{12} & k_{11} \end{pmatrix}, \qquad (45) $$

where we take into account the symmetry of the coefficients of the mass and stiffness matrices for the degrees of freedom contributing to these matrices. In this case the matrices $M^e$ and $K^e$ depend on 12 unknown coefficients $m_{11}, m_{12}, m_{13}, m_{14}, m_{22}, m_{23}$ and $k_{11}, k_{12}, k_{13}, k_{14}, k_{22}, k_{23}$. With the help of (3), the coefficients of the stencil equation (40) can be expressed in terms of the coefficients of the matrices $M^e$ and $K^e$ as follows:

$$ m_1 = m_{14}, \quad m_2 = 2m_{13}, \quad m_3 = 2m_{12} + m_{23}, \quad m_4 = 2(m_{11} + m_{22}), \qquad (46) $$

$$ k_1 = k_{14}, \quad k_2 = 2k_{13}, \quad k_3 = 2k_{12} + k_{23}, \quad k_4 = 2(k_{11} + k_{22}). \qquad (47) $$

Solving simultaneously Eqs. (42), (43), (46) and (47) we can find the mass and stiffness matrices for the cubic isogeometric elements with reduced dispersion. They have the following form:

$$ M^e = h \begin{pmatrix} a_1 & \frac{329913a_2}{1238} - \frac{a_3}{2} & \frac{18387a_2}{619} & a_2 \\ \frac{329913a_2}{1238} - \frac{a_3}{2} & \frac{362654a_2}{619} - a_1 & a_3 & \frac{18387a_2}{619} \\ \frac{18387a_2}{619} & a_3 & \frac{362654a_2}{619} - a_1 & \frac{329913a_2}{1238} - \frac{a_3}{2} \\ a_2 & \frac{18387a_2}{619} & \frac{329913a_2}{1238} - \frac{a_3}{2} & a_1 \end{pmatrix}, \qquad (48) $$

$$ K^e = \frac{1}{h} \begin{pmatrix} \frac{1626310a_2}{1857} - a_4 & \frac{-619a_5 - 263655a_2}{1238} & \frac{-130977a_2}{619} & \frac{-49483a_2}{1857} \\ \frac{-619a_5 - 263655a_2}{1238} & a_4 & a_5 & \frac{130977a_2}{619} \\ \frac{-130977a_2}{619} & a_5 & a_4 & \frac{-619a_5 - 263655a_2}{1238} \\ \frac{-49483a_2}{1857} & \frac{-130977a_2}{619} & \frac{-619a_5 - 263655a_2}{1238} & \frac{1626310a_2}{1857} - a_4 \end{pmatrix}, \qquad (49) $$

where $a_j$ ($j = 1, 2, 3, 4, 5$) are five arbitrary coefficients. The dispersion error is independent of the values of these coefficients. Some coefficients $a_j$ can be found from Eqs. (25) and (26) with $U_0 = (u_0, u_0, u_0, u_0)^T$ and $V_0 = (v_0, v_0, v_0, v_0)^T$ where $U_0$ and $V_0$ are the vectors with four equal components. Using Eq. (25) and the mass matrix from Eq. (48) we can find that $a_2 = 619/1459920$. Using Eq. (26) and the stiffness matrix from Eq. (49) we can find that $a_4 = (25029 - 69520a_5)/139040$. For $a_2 = 619/1459920$ and $a_4 = (25029 - 69520a_5)/139040$, Eqs. (48) and (49) can be rewritten as follows:

$$ M^e = \frac{h}{5040} \begin{pmatrix} 5040a_1 & \frac{989739}{1738} - 2520a_3 & \frac{55161}{869} & \frac{1857}{869} \\ \frac{989739}{1738} - 2520a_3 & \frac{1087962}{869} - 5040a_1 & 5040a_3 & \frac{55161}{869} \\ \frac{55161}{869} & 5040a_3 & \frac{1087962}{869} - 5040a_1 & \frac{989739}{1738} - 2520a_3 \\ \frac{1857}{869} & \frac{55161}{869} & \frac{989739}{1738} - 2520a_3 & 5040a_1 \end{pmatrix}, \qquad (50) $$
$$ K^e = \frac{1}{120h} \begin{pmatrix} 60a_5 + \frac{239399}{10428} & -60a_5 - \frac{37665}{3476} & -\frac{1701}{158} & -\frac{7069}{5214} \\ -60a_5 - \frac{37665}{3476} & \frac{75087}{3476} - 60a_5 & 120a_5 & -\frac{1701}{158} \\ -\frac{1701}{158} & 120a_5 & \frac{75087}{3476} - 60a_5 & -60a_5 - \frac{37665}{3476} \\ -\frac{7069}{5214} & -\frac{1701}{158} & -60a_5 - \frac{37665}{3476} & 60a_5 + \frac{239399}{10428} \end{pmatrix}. \quad (51) $$

**Remark.** Similar to the explanations in Section 1.2, the 12th order of the dispersion error in Eq. (44) is the maximum order for all cubic isogeometric elements.

### 1.4.1. Numerical dispersion of the cubic isogeometric elements with the lumped mass matrix

The analysis of the numerical dispersion of the cubic isogeometric elements with the lumped mass matrix is similar to that for the quadratic elements in Section 1.2.2. Let us consider the stencil equation (40) with $m_1 = m_2 = m_3 = 0$. For the conventional stiffness matrix with $k_1 = -42, k_2 = -1008, k_3 = -630, k_4 = 3360$, the coefficient $m_4$ can be found by the analysis of the dispersion equation (as described in the previous Sections) and is equal to $m_4 = 1$. In this case, the dispersion equation has the following form

$$ \frac{k}{k_h} = \frac{1}{(k_h h)} \frac{\sqrt{-15 \cos(k_h h) - 24 \cos(2(k_h h)) - \cos(3(k_h h)) + 40}}{2\sqrt{15}}. \quad (52) $$

The numerical dispersion error $e = \frac{k_h - k}{k}$ given by Eq. (52) is plotted in Fig. 3; see curve 2. It is interesting to note that this numerical dispersion error for the cubic elements is greater than that for the quadratic elements with the lumped mass matrix; see curve 1 in Fig. 3. A Taylor series of the right-hand side of Eq. (52) yields

$$ \frac{k}{k_h} = 1 - \frac{(k_h h)^2}{6} + \frac{(k_h h)^4}{80} + O\left((k_h h)^6\right), \quad (53) $$

i.e., the dispersion error has the second order of accuracy for the lumped mass matrix (see also Fig. 3(b)). Using Eq. (46), we can find the following form of the lumped mass matrix:

$$ M^e = h \begin{pmatrix} a_1 & -\frac{a_2}{2} & 0 & 0 \\ -\frac{a_2}{2} & \frac{1}{2} - a_1 & a_2 & 0 \\ 0 & a_2 & \frac{1}{2} - a_1 & -\frac{a_2}{2} \\ 0 & 0 & -\frac{a_2}{2} & a_1 \end{pmatrix} \quad (54) $$

with two arbitrary coefficients $a_1$ and $a_2$. The dispersion error is independent of the values of these coefficients. Due to Eq. (3) for the calculation of the global mass matrix, the coefficient $a_2$ in Eq. (54) should be zero in order to get the diagonal global mass matrix (or the elemental mass matrix in Eq. (54) should be modified for the elements close to the boundary). It is important to note that according to the analysis of the stencil equation (40) with the lumped mass matrix, there are no additional degrees of freedom that allow the reduction of the dispersion error with the conventional stiffness matrix. Independent of the techniques for the calculations of the lumped mass matrix, the best accuracy corresponds to the second order of the dispersion error.

However, the accuracy of the numerical results with the lumped mass matrix can be improved by the modification of the stiffness matrix. Let us consider again the stencil equation (40) with $m_1 = m_2 = m_3 = 0$ and undefined coefficients $k_j$ ($j = 1, 2, 3, 4$). Using the approach similar to that in Section 1.4 and equating the first term of a Taylor series of the dispersion equation to unity, and the following three terms to zero, we can find that $m_4 = 90a_1$, $k_1 = -a_1$, $k_2 = 27a_1/2$, $k_3 = -135a_1$, $k_4 = -2k_1 - 2k_2 - 2k_3 = 245a_1$ with one arbitrary coefficient $a_1$. The numerical dispersion error $e = \frac{k_h - k}{k}$ for the lumped mass matrix and the modified stiffness matrix is plotted in Fig. 3;
see curve 4. In this case, a Taylor series of the dispersion equation has the following form

$$ \frac{k}{k_h} = 1 - \frac{(k_h h)^6}{1120} + \frac{(k_h h)^8}{7200} + O\left((k_h h)^{10}\right), \qquad (55) $$

i.e., we can get the 6th order for the dispersion error (see also Fig. 3(b)). Using Eqs. (46) and (47), we can find the following form of the lumped mass matrix and the modified stiffness matrix:

$$ M^e = h \begin{pmatrix} a_1 & -\frac{a_2}{2} & 0 & 0 \\ -\frac{a_2}{2} & -a_1 - 45a_3 & a_2 & 0 \\ 0 & a_2 & -a_1 - 45a_3 & -\frac{a_2}{2} \\ 0 & 0 & -\frac{a_2}{2} & a_1 \end{pmatrix}, \qquad (56) $$

$$ K^e = \frac{1}{h} \begin{pmatrix} \frac{1}{2}(-2a_4 - 245a_3) & \frac{1}{2}(135a_3 - a_5) & -\frac{27a_3}{4} & a_3 \\ \frac{1}{2}(135a_3 - a_5) & a_4 & a_5 & -\frac{27a_3}{4} \\ -\frac{27a_3}{4} & a_5 & a_4 & \frac{1}{2}(135a_3 - a_5) \\ a_3 & -\frac{27a_3}{4} & \frac{1}{2}(135a_3 - a_5) & \frac{1}{2}(-2a_4 - 245a_3) \end{pmatrix} \qquad (57) $$

with five arbitrary coefficients $a_j$ ($j = 1, 2, 3, 4, 5$) that do not affect the dispersion error. The coefficients $a_3 = -1/90$ and $a_4 = 1/40(27 - 20a_5)$ can be found with the help of Eqs. (25) and (26). With these coefficients, Eqs. (56) and (57) can be rewritten as:

$$ M^e = \frac{h}{5040} \begin{pmatrix} 5040a_1 & -2520a_2 & 0 & 0 \\ -2520a_2 & 2520(1 - 2a_1) & 5040a_2 & 0 \\ 0 & 5040a_2 & 2520(1 - 2a_1) & -2520a_2 \\ 0 & 0 & -2520a_2 & 5040a_1 \end{pmatrix}, \qquad (58) $$

$$ K^e = \frac{1}{120h} \begin{pmatrix} 60a_5 + \frac{247}{3} & -30(2a_5 + 3) & 9 & -\frac{4}{3} \\ -30(2a_5 + 3) & 81 - 60a_5 & 120a_5 & 9 \\ 9 & 120a_5 & 81 - 60a_5 & -30(2a_5 + 3) \\ -\frac{4}{3} & 9 & -30(2a_5 + 3) & 60a_5 + \frac{247}{3} \end{pmatrix}. \qquad (59) $$

## 2. Numerical examples

In order to show the advantages of the new high-order isogeometric elements with reduced dispersion, below we solve three 1-D wave propagation problems by the conventional and new high-order isogeometric elements. The first two problems have continuous solutions and will be solved without the filtering stage. The third problem is related to impact loading and will be solved by the two-stage time-integration technique with the filtering stage; see [33,36].

### 2.1. A standing wave in 1-D elastic bar

Similar to [42], the problem of a standing wave is used for the demonstration of the order of convergence of the new numerical technique. Let us consider an elastic bar of length $L = 1$. The wave velocity is chosen to be $c = 1$. A standing wave in the 1-D case can be described by the following exact solution to the wave equation:

$$ u(x, t) = \sin\left(\frac{a\pi x}{L}\right) \cos\left(\frac{a\pi t}{L}\right), \qquad (60) $$
<figure><img src="image_7.png" /><figcaption>Fig. 6. The error in velocity $e_v = |v_{exact}(x = L/2, t = T) - v_{num}(x = L/2, t = T)|$ as a function of the mesh size $h$ in the logarithmic scale. Curves 1 and 2 correspond to the conventional quadratic and cubic isogeometric elements. Curves 3 and 4 correspond to the new quadratic and cubic isogeometric elements. The slopes of the curves at small $h$ show the order of convergence of the corresponding techniques.</figcaption></figure>

where $a = 5$ is used. The initial conditions at time $t = 0$ and the boundary conditions in terms of displacements at $x = 0$ and $x = L$ are selected according to the exact solution, Eq. (60); i.e., the initial displacements are $u(x, t = 0) = \sin(\frac{ax}{L})$, the initial velocities are zero and the two ends of the bar are fixed. The observation time is selected to be $T = L/a$. According to Eq. (60), the displacements at this observation time change the sign and the velocities at all points become zero again. The problem is solved by the conventional and new isogeometric elements on meshes with uniformly spaced control points. For the time integration the trapezoidal rule is used with very small time increments at which the error in time is very small and can be neglected. This means that the difference between numerical and analytical solutions is only related to the space-discretization error. Fig. 6 shows the convergence of the error in the velocity $e_v = |v_{exact}(x = L/2, t = T) - v_{num}(x = L/2, t = T)|$ in the center of the bar at mesh refinement where $v_{exact}(x = L/2, t = T)$ and $v_{num}(x = L/2, t = T)$ are the exact and numerical velocities in the center of the bar at the observation time $T$. We should mention that the maximum error in the velocity occurs in the center of the bar. $h$ in Fig. 6 is the distance between uniformly spaced control points. At the same $h$, the meshes with quadratic and cubic isogeometric elements include the same number of degrees of freedom. The results in Fig. 6 are plotted in the logarithmic scale. Therefore, the slopes of the curves at small $h$ shown in Fig. 6 correspond to the order of convergence (the order of accuracy) of the conventional and new isogeometric elements. The new isogeometric elements yield much more accurate results than those obtained by the conventional elements. Moreover, at the same number of degrees of freedom (at the same $h$), the new quadratic elements are more accurate compared to the conventional quadratic and cubic elements (see curves 1, 2 and 3 in Fig. 6). The results in Fig. 6 are in good agreement with the theoretical order of accuracy of the conventional and new isogeometric elements reported in the previous sections of the paper.

## 2.2. Propagation of sinusoidal pulse in 1-D elastic bar

Let us consider an elastic bar of length $L = 4$ with free ends. The wave velocity is chosen to be $c = 1$. Zero initial displacements $u(0, t) = 0$ and the following distribution of the initial velocity are applied: $v(x, 0) = 0$ at $x < 1.6$ and $x > 2$ as well as $v(x, 0) = \sin(\Omega(x - 1.6) - \pi/2) + 1$ at $1.6 \le x \le 2$ with frequency $\Omega = 5\pi$; see Fig. 7. According to the exact D'Alembert's solution to this problems, the initial velocity sinusoidal pulse splits into two velocity sinusoidal pulses of the same frequency $\Omega$ and a half amplitude that propagate in the opposite directions with wave velocity $c = 1$. For example, the exact solutions at the observation times $T = 1$ and $T = 17$ are shown in Fig. 8; see curves 1 (it can be checked that the exact solutions at the observation times $T = 1$ and $T = 17$ coincide). Due to the small shift of the center of the initial velocity pulse (located at $x = 1.8$; see Fig. 7) from the center of the bar the exact solution is not symmetric with respect the center of the bar.

The problem is solved by the conventional and new isogeometric elements on a mesh with uniformly spaced control points with 51, 101 and 201 degrees of freedom. For the time integration we use the trapezoidal rule with very small time increments at which the error in time is very small and can be neglected. This means that the difference between numerical and analytical solutions is only related to the space-discretization error. The numerical results for the velocity distribution at the observation times $T = 1$ and $T = 17$ obtained by the conventional quadratic and cubic isogeometric elements and by the new quadratic isogeometric elements with 51 degrees of freedom are
<figure><img src="image_8.png" /><figcaption>Fig. 7. The distribution of the initial velocity in a 1-D bar of length L = 4.</figcaption></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 8. The velocity distribution in the 1-D bar at the observation times $T = 1$ (a, c, e) and $T = 17$ (b, d, f). (a, b), (c, d) and (e, f) correspond to the numerical results obtained by the conventional quadratic and cubic isogeometric elements and by the new quadratic isogeometric elements with reduced dispersion, respectively (51 degrees of freedom are used for all results). Curves 1 and 2 correspond to the exact and numerical results, respectively.</figcaption></figure>

shown in Fig. 8. As can be seen, despite the same exact solutions at the observation times $T = 1$ and $T = 17$, the numerical solutions at these times are different. At the smaller observation time $T = 1$, the numerical results are close to the exact solutions; see Fig. 8(a, c, e). At the larger observation time $T = 17$ the numerical results are much less accurate for all elements due to the numerical dispersion error (the space-discretization error is an increasing function of time even at the exact time-integration). The numerical dispersion leads to the shift in space of
<figure><img src="image_15.png" /><figcaption>Fig. 9. The comparison of the velocity distribution in the 1-D bar at the observation time $T = 17$. Curves 1, 2, 3 and 4 correspond to the exact solution and the numerical results obtained by the conventional quadratic and cubic isogeometric elements and by the new quadratic isogeometric elements with reduced dispersion, respectively (51 degrees of freedom are used for all results); see also Fig. 8.</figcaption></figure>

the center of the sinusoidal pulse as well as to the appearance of small spurious oscillations before and after the main oscillation; compare curves 1 and 2 in Fig. 8(b, d, f). The results also clearly show that at the same numbers of degrees of freedom, the new quadratic isogeometric elements with reduced dispersion significantly improve the accuracy of the numerical results compared with the conventional quadratic and cubic isogeometric elements; see also Fig. 9. Despite a higher accuracy, at the same numbers of degrees of freedom the new quadratic isogeometric elements with reduced dispersion require no additional computational costs compared with the conventional quadratic isogeometric elements and require less computation costs compared with the conventional cubic isogeometric elements (due to a smaller bandwidth of the stencil equation and the corresponding matrix of the system of algebraic equations for the quadratic elements compared to those for the cubic elements). The convergence of the numerical results obtained by the new and conventional quadratic isogeometric elements at mesh refinement is shown in Fig. 10. As can be seen from Fig. 10, the new quadratic isogeometric elements converge to the exact solution much faster compared to the conventional quadratic isogeometric elements. Fig. 11 shows that at the larger observation time $T = 17$, the new cubic isogeometric elements with reduced dispersion are much more accurate compared to the conventional cubic isogeometric elements (see curves 2) and more accurate compared to the new quadratic elements at the same number of degrees of freedom (see curves 2 in Figs. 10(b) and 11(b)).

**Remark.** It is interesting to note that at the derivation of the mass and stiffness matrices of the new isogeometric elements, the basis functions have not been used. Therefore, after the solution of the global semidiscrete system related to the new quadratic isogeometric elements with the non-diagonal mass matrix, we have also plotted the numerical results assuming that the global control variables yield the values of the field function $u(x, t)$ at the nodes of a uniform mesh (i.e., similar to the conventional finite elements, at these nodes one basis function equals unity, all other basis functions equal zero). The results with this post-processing procedure are shown in Fig. 12; see curves 1. As can be seen, in this case the field function $u(x, t)$ is very close to that obtained by the new isogeometric elements with the conventional isogeometric basis functions; see curves 2 in Fig. 12.

## 2.3. 1-D elastic bar under suddenly applied load

The 1-D impact problem considered below is a good benchmark wave propagation problem; e.g., see [36]. All low and high frequencies are excited for this problem. Despite a simple analytical solution, this problem is challenging for any numerical approach due to the appearance of spurious oscillations in numerical solutions; e.g., see [33,43] and Figs. 14(a, b) and 15(a, b). To treat this issue, we use the two-stage time integration technique developed in our papers [33,36,39]. The filtering stage of this technique includes the time continuous Galerkin (TCG) method with 10 time increments (5 positive plus 5 negative time increments) the size of which is calculated according to the special formulas; see [33,36]. We should also mention that due to a simple analytical solution, it is easy to compare the accuracy of the numerical solutions of this problem when the spurious oscillations are removed from the results. The numerical solutions of this problem obtained by the conventional high-order finite elements, spectral elements and isogeometric elements have been compared in our paper [39]. The numerical results in [39] also show the computational efficiency of high-order elements compared with the linear elements at the increase in the observation time (the numerical dispersion significantly affects the error increase at large observation times even for the exact time integration).
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 10. The velocity distribution in the 1-D bar at the observation time $T = 17$ obtained by the conventional quadratic isogeometric elements (a, c, e) and by the new quadratic isogeometric elements with reduced dispersion (b, d, f). (a, b), (c, d) and (e, f) correspond to the numerical results obtained on the meshes with 51, 101 and 201 degrees of freedom, respectively. Curves 1 and 2 correspond to the exact and numerical results, respectively.</figcaption></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 11. The velocity distribution in the 1-D bar at the observation time $T = 17$ obtained by the conventional cubic isogeometric elements (a) and by the new cubic isogeometric elements with reduced dispersion (b). 51 degrees of freedom are used for all results. Curves 1 and 2 correspond to the exact and numerical results, respectively.</figcaption></figure>

Let us consider an elastic bar of length $L = 4$ under a suddenly applied load (see Fig. 13). The wave velocity is chosen to be $c = 1$. The following boundary conditions are applied: $u(0, t) = t$ (which corresponds to the velocity $v(0, t) = v_0 = 1$) and $u(4, t) = 0$ (which corresponds to the velocity $v(4, t) = 0$). Zero initial conditions are used; i.e., $u(x, 0) = v(x, 0) = 0$. The observation time is chosen to be $T = 18$. During this time the waves travel within the
<figure><img src="image_24.png" /><figcaption>Fig. 12. The velocity distribution in the 1-D bar at the observation time $T = 17$ obtained by the solution of the global system of semidiscrete equations with 51 (a) and 101 (b) degrees of freedom (the same as for the new quadratic isogeometric elements with reduced dispersion). Curves 1 correspond to the post-processing similar to that for the finite elements (see the text for more explanations). Curves 2 correspond to the new quadratic isogeometric elements with reduced dispersion.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 13. Impact of an elastic bar of length $L = 4$ against a rigid wall.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 14. The velocity distribution along the bar for the 1-D impact problem at time $T = 18$ after several reflections from the ends of the bar (a–d). Curve 1 corresponds to the analytical solution. Curves 2 and 3 correspond to the numerical results without (a, b) and with (c, d) the filtering and are obtained with the conventional quadratic isogeometric elements and the new quadratic isogeometric elements with reduced dispersion, respectively (101 degrees of freedom are used for all results). (b, d) zoom (a, c) in the vicinity of the wave front ($1.8 < x < 2.2$).</figcaption></figure>

bar with two reflections from each end of the bar. The analytical solution to this problem for the velocity at $T = 18$ is: $v_a(x, 18) = 1$ for $0 \le x < 2$ and $v_a(x, 18) = 0$ for $2 < x \le 4$ (at $x = 2$, a jump in the velocity occurs).

The problem is solved by the conventional and new quadratic and cubic isogeometric elements on a mesh with uniformly spaced control points with 101 degrees of freedom. For the time integration at the stage of basic
<figure><img src="image_27.png" /><figcaption>Fig. 15. The velocity distribution along the bar for the 1-D impact problem at time $T = 18$ after several reflections from the ends of the bar (a–d). Curve 1 corresponds to the analytical solution. Curves 2, 3 and 4 correspond to the numerical results without (a, b) and with (c, d) the filtering and are obtained with the conventional cubic isogeometric elements and the new quadratic and cubic isogeometric elements with reduced dispersion, respectively (101 degrees of freedom are used for all results). (b, d) zoom (a, c) in the vicinity of the wave front ($1.8 < x < 2.2$).</figcaption></figure>

computations we use the trapezoidal rule with very small time increments at which the error in time is very small and can be neglected. This means that the difference between numerical and analytical solutions is only related to the space-discretization error. The details related to the filtering stage can be found in our papers [33,36,39]. The numerical results for the velocity distribution at time $T = 18$ obtained by the conventional and new quadratic and cubic isogeometric elements are shown in Figs. 14 and 15. As can be seen, after basic computations all results include spurious high-frequency oscillations; see Figs. 14(a, b) and 15(a, b). These oscillations are smaller for the isogeometric elements with the smaller dispersion error. However, it is difficult to compare the accuracy of the numerical results after the stage of basic computations due to the spurious oscillations. After the filtering stage, the comparison of accuracy is simple due to the simple analytical solution of the 1-D impact problems (curves 1 in Figs. 14 and 15). We can clearly see that the new high-order isogeometric elements with reduced dispersion significantly improve the accuracy of the numerical results compared with the conventional high-order isogeometric elements; compare curves 2 and 3 in Fig. 14(c, d) and curves 2 and 3, 4 in Fig. 15(c, d). We should also mention that the new quadratic isogeometric elements with reduced dispersion yield more accurate results than those obtained by the conventional cubic isogeometric elements; see curves 2 and 3 in Fig. 15(c, d). Despite a higher accuracy, the new high-order isogeometric elements with reduced dispersion do not require additional computational costs compared with the conventional high-order isogeometric elements.

**Remark.** It was shown in [3,5] that uniformly spaced control points yield more accurate results for structural dynamics problems than a linear parametrization with the constant Jacobian determinant. However, even with uniformly spaced control points, the numerical results include spurious oscillations in basic computations (see Figs. 14(a, b) and 15(a, b)) and require the filtering stage.

### 3. Concluding remarks

The Part 1 of paper considers the development of the new high-order isogeometric elements with reduced dispersion in the 1-D case. The new findings of the paper can be summarized as follows:
* By the analysis of the dispersion error of the stencil equation with arbitrary coefficients we have shown that these coefficients can be found by the minimization of the order of the dispersion error. In this case, the order of the dispersion error can be increased from the order $2p$ for the conventional high-order isogeometric elements to the order $4p$ for the new elements ($p$ is the order of the polynomial approximations). Because all coefficients of the stencil equation can be found from the analysis of the dispersion error, the order $4p$ of the dispersion error is maximum possible for the considered form of the stencil equation (that is related to the support of basis functions by a specific number of elements) for the corresponding high-order isogeometric elements. We should mention that we have never seen in the literature such a significant improvement of the order of the dispersion error for the conventional finite elements, the spectral elements or the isogeometric elements. For example, for the linear and high-order finite and isogeometric elements (e.g., see [7,17,22–24]), the order of the dispersion error has been improved from the order $2p$ to the order $2p + 2$ and this leads to a significant decrease in the computation time at a given accuracy. The increase in accuracy for the new isogeometric elements with the order $4p$ of the dispersion error is much higher compared with the known approaches. This will lead to a much more significant reduction in the computation time at a given accuracy.

* Many wave propagation problems require small time increments for the time integration and the use of explicit time-integration methods with the lumped mass matrix is very computationally efficient. However, for the high-order isogeometric elements, the lumped mass matrix yields the second order of the dispersion error (e.g., see [3,5]) and makes the application of the high-order isogeometric elements computationally inefficient. Therefore, a special iterative procedure for linear dynamics problems that requires additional computational costs was suggested in [30] in order to obtain the same order of the dispersion error for the lumped and consistent mass matrices. In the paper we have shown that independent of the procedures for the calculation of the lumped mass matrix, the order of the dispersion error cannot not be improved with the conventional stiffness matrix and no additional iterations. We also showed that the dispersion error for the lumped mass matrix can be improved from the second order to the order $2p$ by the modification of the stiffness matrix (without additional computational costs).

* For the linear basis functions, the conventional isogeometric elements coincide with the linear finite elements. The analysis of the stencil equation for the linear elements shows that the reduction of the numerical dispersion error for the linear elements from order two to order four in the known numerical techniques (e.g., see[13,14,17,22]) is maximum possible for all linear elements and cannot be improved (the dispersion analysis of the linear elements is similar to that for the quadratic and cubic elements presented in Section 1).

* By the solution of wave propagation problems in Section 2 we have shown the computational efficiency of the new high-order isogeometric elements with reduced dispersion. For example, the new quadratic isogeometric elements yield much more accurate results than those obtained by the conventional cubic isogeometric elements; e.g., see Figs. 6, 9 and 15.

* Despite a good accuracy, the conventional and new high-order isogeometric elements yield spurious oscillations under high-frequency and impact loadings. In the paper we have shown that by the application of the two-stage time-integration technique developed in our recent papers [33,36,39] we can get accurate numerical results with the isogeometric elements as well; e.g., see Figs. 14 and 15. This technique quantifies and filters spurious oscillations.

* We should also mention that the approach developed in the paper can be directly applied to other space-discretization techniques that have the same form of the stencil equations (e.g., see Eqs. (17) and (40)); e.g., to the high-order finite-difference method in order to reduce the dispersion error.

The development of the new high-order isogeometric elements with reduced dispersion for wave propagation problems in the 2-D case with the corresponding numerical examples is considered in Part 2 of the paper; see [1].

## Acknowledgments

The research has been supported in part by the Air Force Office of Scientific Research (contract FA9550-16-1-0177) and by Texas Tech University.

## Appendix. Local truncation error in space for the new and conventional isogeometric elements

Based on the semi-discrete system of Eqs. (2) in the time domain and Eq. (3), the stencil equation for the degree of freedom $u_A$ in the time domain for the new and conventional quadratic elements has the following form (the derivation
is similar to that for Eq. (17)):

$$h^2[m_1(\ddot{u}_{A-2} + \ddot{u}_{A+2}) + m_2(\ddot{u}_{A-1} + \ddot{u}_{A+1}) + m_3\ddot{u}_A] + c^2[k_1(u_{A-2} + u_{A+2}) \\ + k_2(u_{A-1} + u_{A+1}) + k_3u_A] = 0, \qquad (A.1)$$

where the coefficients $m_j$ and $k_j$ ($j = 1, 2, 3$) are given by Eq. (20) for the new elements with reduced dispersion and are $m_1 = 1, m_2 = 26, m_3 = 66, k_1 = -20, k_2 = -40, k_3 = 120$ for the conventional quadratic isogeometric elements; see Section 1.1. Similarly, the stencil equation for the degree of freedom $u_A$ in the time domain for the new and conventional cubic elements has the following form (the derivation is similar to that for Eq. (40)):

$$h^2[m_2(\ddot{u}_{A-2} + \ddot{u}_{A+2}) + m_3(\ddot{u}_{A-1} + \ddot{u}_{A+1}) + m_4\ddot{u}_A + m_1(\ddot{u}_{A-3} + \ddot{u}_{A+3})] \\ + c^2[k_2(u_{A-2} + u_{A+2}) + k_3(u_{A-1} + u_{A+1}) + k_4u_A + k_1(u_{A-3} + u_{A+3})] = 0, \qquad (A.2)$$

where the coefficients $m_j$ and $k_j$ ($j = 1, 2, 3, 4$) are given by Eqs. (42) and (43) for the new elements with reduced dispersion and are $m_1 = 1, m_2 = 120, m_3 = 1191, m_4 = 2416, k_1 = -42, k_2 = -1008, k_3 = -630, k_4 = 3360$ for the conventional cubic isogeometric elements; see Section 1.3. Let us analyze the local truncation error in space for the stencil equations (A.1) and (A.2) with the corresponding coefficients $m_i$ and $k_i$. The local truncation error is the residual of the stencil equation when the unknowns $u_i$ ($i = A-3, A-2, A-1, A, A+1, A+2, A+3$) correspond to the exact solution. Because at the derivation of the coefficients $m_i$ and $k_i$ for the new approach (see Eq. (20) and Eqs. (42) and (43)) the basis functions have not been used we assume that the control variables $u_i$ yield the values of the field function $u(x, t)$ at the nodes of a uniform mesh (similar to the conventional finite elements); see also the Remark and Fig. 12 in Section 2.2. Under this assumption, let us expand $u_i$ ($i = A-3, A-2, A-1, A+1, A+2, A+3$) into a Taylor series at small $h \ll 1$ as follows:

$$u_{A\pm l} = u_A + \frac{\partial u_A}{\partial x}(\pm l h) + \frac{\partial^2 u_A}{\partial x^2}(\pm l h)^2 + \frac{\partial^3 u_A}{\partial x^3}(\pm l h)^3 + \frac{\partial^4 u_A}{\partial x^4}(\pm l h)^4 + \dots \qquad (A.3)$$

with $l = 1, 2$ for the quadratic isogeometric elements and $l = 1, 2, 3$ for the cubic isogeometric elements. The exact solution $u_A$ to Eq. (1) meets the following equations:

$$\frac{\partial^2 u_A}{\partial t^2} - c^2 \frac{\partial^2 u_A}{\partial x^2} = 0, \qquad (A.4)$$

$$\frac{\partial^{2j+2} u_A}{\partial x^{2j} \partial t^2} - c^2 \frac{\partial^{2j+2} u_A}{\partial x^{2j+2}} = 0, \qquad (A.5)$$

with $j = 1, 2, 3, 4, \dots$. Inserting Eqs. (A.3)–(A.5) into the stencil equation (A.1) we will get the following local truncation error in space $e_{new}$ for the new quadratic isogeometric elements:

$$e_{new} = \frac{79c^2a_1}{2971080} \frac{\partial^{10}u_A}{\partial x^{10}} h^{10} + O(h^{12}) \qquad (A.6)$$

and $e_{conv}$ for the conventional quadratic isogeometric elements:

$$e_{conv} = -\frac{c^2 \partial^6 u_A}{6 \partial x^6} h^6 + O(h^8), \qquad (A.7)$$

i.e., the new approach improves the order of the local truncation error in space by four orders (similar to the improvement of the order of the numerical dispersion error for the new quadratic isogeometric elements).

Similarly, inserting Eqs. (A.3)–(A.5) into the stencil equation (A.2) we will get the following local truncation error in space $e_{new}$ for the new cubic isogeometric elements:

$$e_{new} = \frac{38223c^2a_1}{247847600} \frac{\partial^{14}u_A}{\partial x^{14}} h^{14} + O(h^{16}) \qquad (A.8)$$

and $e_{conv}$ for the conventional cubic isogeometric elements:

$$e_{conv} = \frac{c^2 \partial^8 u_A}{6 \partial x^8} h^8 + O(h^{10}), \qquad (A.9)$$
i.e., the new approach improves the order of the local truncation error in space by six orders (similar to the improvement of the order of the numerical dispersion error for the new cubic isogeometric elements).

## References

[1] A.V. Idesman, B. Dey, Optimal reduction of numerical dispersion for wave propagation problems. Part 2: Application to 2-D isogeometric elements, Comput. Methods Appl. Mech. Engrg. (2016) 1-38 submitted for publication.

[2] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1987.

[3] T. Hughes, A. Reali, G. Sangalli, Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS, Comput. Methods Appl. Mech. Engrg. 197 (49-50) (2008) 4104-4124.

[4] G. Seriani, S.P. Oliveira, Dispersion analysis of spectral element methods for elastic wave propagation, Wave Motion 45 (6) (2008) 729-744.

[5] J. Cottrell, A. Reali, Y. Bazilevs, T. Hughes, Isogeometric analysis of structural vibrations, Comput. Methods Appl. Mech. Engrg. 195 (41-43) (2006) 5257-5296.

[6] H. Ahmadian, M. Friswell, J. Mottershead, Minimization of the discretization error in mass and stiffness formulations by an inverse method, Internat. J. Numer. Methods Engrg. 41 (2) (1998) 371-387.

[7] M. Ainsworth, H.A. Wajid, Optimally blended spectral-finite element scheme for wave propagation and nonstandard reduced integration, SIAM J. Numer. Anal. 48 (1) (2010) 346-371.

[8] H.P. Cherukuri, Dispersion analysis of numerical approximations to plane wave motions in an isotropic elastic solid, Comput. Mech. 25 (4) (2000) 317-328.

[9] W. Dauksher, A.F. Emery, Solution of elastostatic and elastodynamic problems with Chebyshev spectral finite elements, Comput. Methods Appl. Mech. Engrg. 188 (1) (2000) 217-233.

[10] D. Gabriel, J. Plesek, R. Kolman, F. Vales, Dispersion of elastic waves in the contact-impact problem of a long cylinder, J. Comput. Appl. Math. 234 (2010) 1930-1936.

[11] M.N. Guddati, B. Yue, Modified integration rules for reducing dispersion error in finite element method, Comput. Methods Appl. Mech. Engrg. 193 (2004) 275-287.

[12] V. Gyrya, K. Lipnikov, M-adaptation method for acoustic wave equation on square meshes, J. Comput. Acoust. 20 (2012) 1250022-1:23.

[13] S. Krenk, Dispersion-corrected explicit integration of the wave equation, Comput. Methods Appl. Mech. Engrg. 191 (2001) 975-987.

[14] K.J. Marfurt, Accuracy of finite difference and finite element modeling of the scalar and elastic wave equation, Geophysics 49 (1984) 533-549.

[15] R. Mullen, T. Belytschko, Dispersion analysis of finite element semidiscretizations of the two-dimensional wave equation, Internat. J. Numer. Methods Engrg. 18 (1982) 11-29.

[16] G. Seriani, S.P. Oliveira, Optimal blended spectral-element operators for acoustic wave modeling, Geophysics 72 (5) (2007) 95-106.

[17] B. Yue, M.N. Guddati, Dispersion-reducing finite elements for transient acoustics, J. Acoust. Soc. Am. 118 (4) (2005) 2132-2141.

[18] F.I. Zyserman, P.M. Gauzellino, Dispersion analysis of a nonconforming finite element method for the three-dimensional scalar and elastic wave equations, Finite Elem. Anal. Des. 41 (13) (2005) 1309-1326.

[19] X. Li, T. Zhu, M. Zhang, G. Long, Seismic scalar wave equation with variable coefficients modeling by a new convolutional differentiator, Comput. Phys. Comm. 181 (11) (2010) 1850-1858.

[20] Z.C. He, A.G. Cheng, G.Y. Zhang, Z.H. Zhong, G.R. Liu, Dispersion error reduction for acoustic problems using the edge-based smoothed finite element method (ES-FEM), Internat. J. Numer. Methods Engrg. 86 (11) (2011) 1322-1338.

[21] D. Wang, R. Tezaur, J. Toivanen, C. Farhat, Overview of the discontinuous enrichment method, the ultra-weak variational formulation, and the partition of unity method for acoustic scattering in the medium frequency regime and performance comparisons, Internat. J. Numer. Methods Engrg. 89 (4) (2012) 403-417.

[22] A. Idesman, M. Schmidt, J.R. Foley, Accurate finite element modeling of linear elastodynamics problems with the reduced dispersion error, Comput. Mech. 47 (2011) 555-572.

[23] A.V. Idesman, D. Pham, Finite element modeling of linear elastodynamics problems with explicit time-integration methods and linear elements with the reduced dispersion error, Comput. Methods Appl. Mech. Engrg. 271 (2014) 86-108.

[24] D. Wang, W. Liu, H. Zhang, Novel higher order mass matrices for isogeometric structural vibration analysis, Comput. Methods Appl. Mech. Engrg. 260 (2013) 92-108.

[25] D. Wang, W. Liu, H. Zhang, Superconvergent isogeometric free vibration analysis of Euler-Bernoulli beams and Kirchhoff plates with new higher order mass matrices, Comput. Methods Appl. Mech. Engrg. 286 (2015) 230-267.

[26] I. Harari, C.L. Nogueira, Reducing dispersion of linear triangular elements for the Helmholtz equation, J. Eng. Mech. 128 (3) (2002) 351-358.

[27] I. Harari, R. Ganel, E. Grosu, Stabilized finite elements for time-harmonic elastic waves, Comput. Methods Appl. Mech. Engrg. 200 (21-22) (2011) 1774-1786.

[28] A. Tkachuk, M. Bischoff, Direct and sparse construction of consistent inverse mass matrices: general variational formulation and application to selective mass scaling, Internat. J. Numer. Methods Engrg. 101 (6) (2015) 435-469.

[29] V. Bokil, N. Gibson, V. Gyrya, D. McGregor, Dispersion reducing methods for edge discretizations of the electric vector wave equation, J. Comput. Phys. 287 (2015) 88-109.

[30] F. Auricchio, L. Beirao da Veiga, T. Hughes, A. Reali, G. Sangalli, Isogeometric collocation for elastostatics and explicit dynamics, Comput. Methods Appl. Mech. Engrg. 249-252 (2012) 2-14.

[31] D.J. Benson, Y. Bazilevs, E. De Luycker, M.-C. Hsu, M. Scott, T.J.R. Hughes, T. Belytschko, A generalized finite element formulation for arbitrary basis functions: From isogeometric analysis to XFEM, Internat. J. Numer. Methods Engrg. 83 (6) (2010) 765-785.

[32] D. Benson, Y. Bazilevs, M.-C. Hsu, T. Hughes, A large deformation, rotation-free, isogeometric shell, Comput. Methods Appl. Mech. Engrg. 200 (1316) (2011) 1367-1378.
[33] A.V. Idesman, Accurate time integration of linear elastodynamics problems, Comput. Model. Eng. Sci. 71 (2) (2011) 111–148.

[34] A. Idesman, S.P. Mates, Accurate finite element simulation and experimental study of elastic wave propagation in a long cylinder under impact loading, Impact Eng. 71 (2014) 1–16.

[35] K.J. Bathe, Finite Element Procedures, Prentice-Hall Inc., Upper Saddle River, New Jersey, 1996.

[36] A.V. Idesman, H. Samajder, E. Aulisa, P. Seshaiyer, Benchmark problems for wave propagation in elastic materials, Comput. Mech. 43 (6) (2009) 797–814.

[37] A.V. Idesman, A new exact, closed-form a-priori global error estimator for second- and higher-order time integration methods for linear elastodynamics, Internat. J. Numer. Methods Engrg. 88 (2011) 1066–1084.

[38] A.V. Idesman, D. Pham, Accurate finite element modeling of acoustic waves, Comput. Phys. Comm. 185 (2014) 2034–2045.

[39] A.V. Idesman, D. Pham, J. Foley, M. Schmidt, Accurate solutions of wave propagation problems under impact loading by the standard, spectral and isogeometric high-order finite elements. comparative study of accuracy of different space-discretization techniques, Finite Elem. Anal. Des. 88 (2014) 67–89.

[40] A.V. Idesman, Accurate finite element modeling of wave propagation in composite and functionally graded materials, Compos. Struct. 117 (2014) 298–308.

[41] L. Dede, C. Jaggli, A. Quarteroni, Isogeometric numerical dispersion analysis for two-dimensional elastic wave propagation, Comput. Methods Appl. Mech. Engrg. 284 (2015) 320–348.

[42] C.-S. Chou, C. Shu, Y. Xing, Optimal energy conserving local discontinuous Galerkin methods for second-order wave equation in heterogeneous media, J. Comput. Phys. 272 (2014) 88–107.

[43] G.M. Hulbert, T.J.R. Hughes, Space-time finite element methods for second-order hyperbolic equations, Comput. Methods Appl. Mech. Engrg. 84 (3) (1990) 327–348.