Optimal spectral approximation of $2n$-order differential operators by mixed isogeometric analysis

We approximate the spectra of a class of $2n$-order differential operators using isogeometric analysis in mixed formulations. This class includes a wide range of differential operators such as those arising in elliptic, biharmonic, Cahn-Hilliard, Swift-Hohenberg, and phase-field crystal equations. The spectra of the differential operators are approximated by solving differential eigenvalue problems in mixed formulations, which require auxiliary parameters. The mixed isogeometric formulation when applying classical quadrature rules leads to an eigenvalue error convergence of order $2p$ where $p$ is the order of the underlying B-spline space. We improve this order to be $2p+2$ by applying optimally-blended quadrature rules developed in \cite{20,52} and this order is an optimum in the view of dispersion error. We also compare these results with the mixed finite elements and show numerically that mixed isogeometric analysis leads to significantly better spectral approximations.


Introduction
The finite element method (FEM) is a widely used and highly effective numerical technique for approximate solutions of boundary value problems. The theory of FEM has been extensively developed during the last 60 years. Nowadays, many different variants of FEM are used for the solutions of various complex linear and nonlinear problems. A special group of FEM techniques is the mixed finite element methods [6,15,16,40,48]. The word "mixed" in this case indicates that there are extra functional spaces used to approximate different solution variables. Traditionally these could be the function and its gradient where each is approximated with a different discrete representation. The literature on mixed finite element methods is quite vast and ranges from the first studies in the 1970s by Brezzi [15], Babuška [7], Crouzeix and Raviart [25] to recent contributions [21,22,40,45,62]. A large amount of research has been devoted to various stabilization techniques for the mixed methods [21,30,33,39,46,49] as well as their error estimates [5,38].
The mixed formulation is naturally used for problems with two independent variables such as velocity and pressure in the Stokes equations. In other popular cases, the second variable is the first or the second derivative of the original variable and approximating it directly in the problem formulation has physical interest. For example, in elasticity, both the stress and the displacement can enter the formulation at the same time. Another example is the Darcy flow equation where the mixed variational formulation is posed in terms of the function spaces L 2 (Ω)/R and H(div, Ω) for the pressure and velocity, respectively (though other options are possible as well, see, e.g., [49]). The mixed finite element framework allows to preserve mass conservation, a property that is important in fluid flow problems [31,43,49] and makes the mixed methods competitive numerical techniques in many engineering applications. The mixed problem discretization leads to a linear algebraic system of the saddle point form The mixed isogeometric analysis is a relatively unexplored topic. Again, most of the work in this area focuses on the fluid flow problems where various isogeometric formulations were applied to the Stokes problem [17,37,42,54]. Recent advances include mixed isogeometric formulations for elasticity [34] and poromechanics [13,32]. In this field, coupling the fluid pressure and the solid deformation, mixed isogeometric formulations violate the inf-sup condition and suffer from numerical instabilities in the incompressible and the nearly incompressible limit. To overcome the numerical instabilities, the projection methods [36] or subdivision-stabilized NURBS discretization can be incorporated [32].
The solution of high-order partial differential equations (PDEs) attracted a lot of attention in recent years. An important subclass is the 2n-order PDE, which for n = 1, 2, 3 reduces to different classical PDEs including Laplacian, Allen-Cahn, biharmonic, Cahn-Hilliard, Swift-Hohenberg, phase-field crystal and other problems. Previous work of numerical analysis typically focused on each of these problems separately (e.g., [14,47,59,61]). Only a limited number of investigations exists for a general 2n-order problem (e.g., for phase-field models [60]). For example, non-degeneracy and uniqueness of its periodic solutions were studied in [58].
Dispersion-minimizing methods based on modified integration rules for reducing dispersion error have been developed previously for classical FEM [2,41] and isogeometric analysis [20,52]. The dispersion error is reduced by blending two standard quadrature rules or using special quadrature rules [28]. For the standard finite and isogeometric elements, these optimal dispersion methods lead to two additional orders of error convergence (superconvergence) in the eigenvalues, while the eigenfunction errors do not degenerate.
In this paper, we utilize the mixed FEM framework for a general 2n-order linear differential eigenvalue problem. We develop the mixed isogeometric framework for these eigenvalue problems and present error analysis for both eigenvalue and eigenfunctions. Optimal blending rules for the mixed isogeometric discretizations of the 2n-order problem are presented up to n = 3.
The rest of this paper is organized as follows. Section 2 describes the differential eigenvalue problems under consideration. Section 3 presents the mixed isogeometric formulation. In Section 4, we present the optimally-blended rules and their eigenvalue error analysis. Section 5 shows numerical examples to demonstrate the performance of the method. Concluding remarks are given in Section 6.

Problem statement
Let Ω ⊂ R d , d = 1, 2, 3, be a bounded domain with Lipschitz boundary ∂Ω. We consider a class of 2n-order linear differential eigenvalue problems: Find λ and non-zero u satisfying Lu = λu in Ω (2.1) and u is subject to appropriate boundary conditions with the differential operator defined as where ∆ is the Laplacian and a m ∈ L ∞ (Ω), m = 0, 1, · · · , n with n being a positive integer. For simplicity, we assume that Ω = [0, 1] d ∈ R d , d = 1, 2, 3 and a m are constants in the following discussions. The operator L in the general equation (2.1) covers many high-order differential operators arising in sciences and engineering. In particular: • For n = 1, a 0 = 0, a 1 = 1, L = −∆ and (2.1) reduces to the Laplacian (or linearized Allen-Cahn [3]) eigenvalue problem.

Mixed formulations
To motivate the presentation of the mixed formulation for any n, we start with n = 2, a 0 = a 1 = 0, a 2 = 1 which (2.1) becomes the biharmonic eigenvalue problem. In this case, the differential equation (2.1) can be recast in a mixed form as a system of equation of second-order −∆u = µ and − ∆µ = λu (3.1) and both u and µ are subject to appropriate boundary conditions. This system of equations is referred to as problem with two unknown fields; see for example [15,23,38]. This new auxiliary parameter has physical meanings. For example, in structural mechanics, the new unknown µ represents the bending moment [57], while in fluid dynamics, when the Stokes equations for viscous steady flow is transformed using stream function this represents the vorticity [4]. If the domain Ω has smooth boundary or it is convex polygonal domain, then the eigenvalue problem (3.1) has infinitely many solutions (λ j , (µ j , u j )) such that [8,26,44] with the Kronecker delta defined as δ jk = 1 when j = k and zero otherwise. Here, the eigenfunctions u j are orthonormal in L 2 (Ω). If (λ j , (µ j , u j )) is an eigenpair of (3.1), then (λ j , u j ) is an eigenpair of (2.1) and µ = −∆u. Hence the regularity of (µ j , u j ) can be inferred from the regularity properties of problem (2.1); see for example [4]. 4

Mixed formulation at continuous level
For arbitrary positive integer n, we set for m = 1, 2, · · · , n − 1 with ψ 0 = u. These auxiliary unknowns allow us to recast the differential equation (2.1) into the mixed form as a system of equations of second-order . The regularity of (ψ 1 j , ψ 2 j , · · · , ψ n−1 j , u j ) can be inferred from the regularity properties of problem (2.1) and we assume sufficient regularity of the problem (2.1). Now, we present the mixed variational formulation for (2.1) at the continuous level. We denote the standard L 2 (Ω)-norm as · 0,Ω ≡ · ≡ · L 2 (Ω) . We adopt the standard Sobolev spaces of integer index s, H s (Ω) and H s 0 (Ω), equipped with the norm · s,Ω ≡ · H s (Ω) ; see [1,23].
We define the bilinear forms where ∇ is the gradient operator. The bilinear form a(·, ·) and b(·, ·) is usually referred as stiffness and mass, respectively. These bilinear forms can be written alternatively as At the continuous level, for simplicity, we assume that the differential equation (2.1) is subject to simply supported boundary conditions The mixed formulations for (2.1) or (3.4) is: Find the eigenpairs (λ, (ψ 1 , ψ 2 , · · · , ψ n−1 , u)) with λ ∈ R and ψ m , u ∈ H 1 0 (Ω), m = 1, 2, · · · , n − 1, satisfying Remark 1. One can also consider other boundary conditions. For example, for n = 2, a 0 = a 1 = 0, a 2 = 1, with homogeneous Dirichlet boundary conditions the mixed variational formulation of the corresponding (2.1) is: Find the eigenpairs (3.10) Herein, n denotes the outward unit normal to the boundary ∂Ω. We refer the readers to [4] for details.

Mixed formulation at discrete level
At discrete level, we specify a finite dimensional approximation space is the span of the B-spline or Lagrange (for FEM) basis functions φ p a of order p. Consequently, the mixed isogeometric analysis (or FEM) of (2.1) with simply supported boundary conditions (3.7) is to seek The definition of the B-spline basis functions in one dimension is as follows. Let X = {x 0 , x 1 , · · · , x m } be a knot vector with knots x j , that is, a nondecreasing sequence of real numbers which are called knots. The j-th B-spline basis function of degree p, denoted as θ j p (x), is defined as [27,51] (3.12) In this paper, for isogeometric analysis, we utilize the B-splines on uniform tensorproduct meshes with non-repeating knots, that is, the B-splines with maximum continuity on uniform meshes, while for finite element method, we utilize the standard Lagrange basis functions. We approximate the eigenfunctions as a linear combination of the Bspline (or Lagrange) basis functions and substitute all the basis functions for V p h in (3.11). This leads to a matrix eigenvalue problem, which is then solved numerically. We give more details of the structures of matrix eigenvalue problem in the following.

Quadrature rules
In practice, we evaluate the integrals involved in the bilinear forms a(·, ·) and b(·, ·) numerically, that is using quadrature rules. On a reference elementK, a quadrature rule is of the form (3.13) whereˆ l are the weights,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 map σ such that K = σ(K), which leads to the correspondence between the functions on K andK. Assuming J K is the corresponding Jacobian of the mapping, (3.13) induces a quadrature rule over the element K given by where l,K = det(J K )ˆ l and n l,K = σ(n l ). For simplicity, we denote by G l the l−point Gauss-Legendre quadrature rule, by L l the l−point Gauss-Lobatto quadrature rule, by R l the l−point Gauss-Radau quadrature rule, and by O p the optimal blending scheme for the p-th order isogeometric analysis with maximum continuity. In one dimension, G l , L l , and R l fully integrate polynomials of order 2l − 1, 2l − 3, and 2l − 2, respectively [9][10][11].
Applying quadrature rules to (3.11), we obtain the approximate form l,K } specifying two (possibly different) quadrature rules. Using quadrature rules, we can write the matrix eigenvalue problem as · · · a n−2 M a n−1 M + a n K , and U , Ψ j , j = 1, · · · , n − 1 are the corresponding representation of the eigenvector as the coefficients of the basis functions. Similar to the standard second-order eigenvalue problem, K and M are referred to as the stiffness and mass matrices resulting from (3.5), respectively. This matrix eigenvalue problem (3.18) has a similar structure with the one obtained by hybrid high-order discretization; see [19,Eqn. 3.13].

Optimally blended quadrature rules
The optimally-blended rules are developed and analyzed for isogeometric analysis in [12,20,28,52] for p ≤ 7 and generalized to arbitrary order p in [29]. We denote the following blendings of Gauss-Legendre rule G p , G p+1 and Gauss-Lobatto rule L p+1 as where τ gg and τ gl are blending parameters.   Table 1 shows the optimal blending parameters for p ≤ 4; see also [20,29]. For optimal blending parameters of higher order p and blending among other quadrature rules, we refer to [29]. These optimally-blended quadrature rules improve spectrum errors significantly, which we show in the next section.

Eigenvalue errors
In this section, following the eigenvalue error estimates established in [20,29,52] for the second-order Laplacian eigenvalue problem, we derive a priori eigenvalue error estimates for (3.18), which can be viewed as a generalization to the 2n-order differential eigenvalue problems. The mixed formulation helps us deliver the optimal error convergence rates for the biharmonic eigenvalue problem.
Following the structure in [29], we seek an approximation of eigenfunction u h in the form where U j are the unknown coefficients which corresponds to the the p-th order polynomial approximation which are to be determined, that is the component of the unknown vector U in (3.18). Using the Bloch wave assumption [50], we write where i 2 = −1 and ω is an approximated frequency. In the view of auxiliary fields (3.3), we denote ψ m h = u h and Ψ m = U when m = 0. In the Bloch wave assumption (4.2), jh resembles the spatial variable x. Using the auxiliary fields defined (3.3), this allows us further assume Bloch wave solutions for their derivatives, that is Ψ m,j = ω 2m e ijωh , m = 0, 1, 2, · · · , n − 1, where Ψ m,j is the component of the unknown vector Ψ m in (3.18). We observe that (4.3) recovers (4.2) when m = 0.
The C p−1 B-spline basis function θ j p has a support over p + 1 elements. Thus, for m = 0, 1, 2, · · · , n − 1, we have where a h (·, ·), b h (·, ·) approximate (or exactly-integrate when we use appropriate quadrature rules) bilinear forms and The symmetry of the B-spline basis functions (on uniform meshes and away from the boundaries) further implies that (4.8) Thus, using this symmetry, Euler's formula, Bloch wave assumptions (4.2) and (4.3), we can deduce the following there holds when a h (·, ·) and b h (·, ·) in (4.6) are approximated using the optimally blended quadrature rules.
Proof. We omit the proof as the identity is proved using the same arguments in the proof of Lemma 6 supplied with Lemma 7 in the paper [29].
Theorem 1. Assuming a h (·, ·) and b h (·, ·) in (4.6) are approximated using the optimally blended quadrature rules, there holds Proof. Let θ j p be a test function for each equation in (3.15). We represent the approximated eigenfunctions ψ m h , m = 0, 1, · · · , n − 1, in the the same way as (4.1). Substituting all these terms into (3.15), we obtain which, by using (4.9), is further simplified as Applying Lemma 1, we have , the mixed formulation with optimally-blended rules leads to an improved eigenvalue error of rate |λ h − λ| ≈ O(h 2p+2 ). When standard quadrature rules are applied, the method retains its optimal rates |λ h − λ| ≈ O(h 2p ). This theorem establishes that the mixed formulation for 2n-order differential eigenvalue maintains the theoretical findings established for the approximation method of the secondorder Laplacian eigenvalue problem. For the generalization to multiple dimensions, we refer to [29] for the second-order Laplacian eigenvalue problem.

Numerical experiments
In this section, we present 1D, 2D, and 3D numerical results which cover the spectral approximations of the biharmonic, Cahn-Hilliard, Swift-Hohenberg, and phase-field crystal operators. We limit our studies to simple geometrical domains and uniform meshes to focus our attention on the numerical aspects of the problem. We utilize both mixed isogeometric and finite elements with p = 1, 2, 3, 4. For our simulations, we consider the unitary domain Ω = [0, 1] d . The exact solutions of (2.1) are given in (2.3)-(2.5) with the parameters specified in Section 2.
First of all, we present the matrix eigenvalue problems associated with the operators described above, which are easily obtained from the general representation (3.18). In particular, with slight abuse of notation, we have for biharmonic eigenvalue problem, for Cahn-Hilliard eigenvalue problem, for Swift-Hohenberg eigenvalue problem, and for phase-field crystal eigenvalue problem, respectively. For the fourth-order differential eigenvalue problem, one can symmetrize the matrix on the left-hand side, that is, the matrix eigenvalue problem (5.1) to (5.3) are equivalent to and respectively. Numerical experiments show that it takes less time to solve (5.5) to (5.7) than solving the original problems (5.1) to (5.3) that have non-symmetric matrices on the left-hand side. We use the symmetric systems to solve the differential eigenvalue problems numerically in the following experiments. 11 As the eigenvalues of (2.1) can be large, we present the relative eigenvalue errors, which for a j-th eigenvalue, is defined as We start our accuracy and convergence studies in 1D and then extend them to multiple dimensions. Figure 1 shows the eigenvalue errors for isogeometric elements of maximum continuity for the biharmonic, Cahn-Hilliard and Swift-Hohenberg equations with homogeneous Dirichlet boundary conditions. The exact solutions are in form of (2.3) in 1D. Similar plots have been previously shown in [24] for approximate eigenvalues of second-order PDEs. As we can see in Figure 1, all three fourth-order operators have quite similar errors when the same boundary conditions are used; the same happens also in multiple dimensions, as we show below.
The case of homogeneous Dirichlet boundary conditions has a particular importance for the biharmonic equation where it corresponds to the case of a simply supported plate. For the differential eigenvalue problems that describe the phase separation process, periodic boundary conditions are more relevant. In the following figures, we consider Dirichlet boundary conditions for the biharmonic equation and periodic boundary conditions for Cahn-Hilliard, Swift-Hohenberg, and phase-field crystal equations.  Figure 2 shows the approximation errors for 1D Cahn-Hilliard and phase-field crystal operators with periodic boundary conditions for p = 2, 3, 4. Swift-Hohenberg spectra are found to be very similar to the Cahn-Hilliard results and are omitted here for brevity. Once again, we can see a clear improvement in the spectral accuracy of IGA discretizations with an increase in p. Due to the use of periodic boundary conditions, outlier modes (large spikes in the errors for j/N close to one (high-frequencies) in the spectra 12 of high-order isogeometric discretizations) are absent in this case. This is related to the fact that outliers are caused by the basis functions with support on the boundaries of the domain [53], which are absent in the periodic case.  Figure 3 compares the eigenvalue and eigenfunction errors between the mixed C 0 finite elements and C 1 isogeometric elements. We observe branching of the finite element spectrum which is typical for high-order C 0 discretizations. Notably, B-spline basis functions of the highest p − 1 continuity do not exhibit such branching patterns on uniform meshes. Similar to the standard (non-mixed) formulation, the large spikes in the approximation errors in the middle of the spectra at the transition point between the acoustic and optical branches are absent in mixed isogeometric discretizations.
In Figure 4, we compare the approximation errors for quadratic elements when using mixed isogeometric analysis with standard Gauss quadratures and optimally-blended rules (3.19). In this case, we also use Dirichlet boundary conditions for the biharmonic equation and periodic boundary conditions for Cahn-Hilliard and phase-field crystal equations. Despite this, all three cases exhibit similar behaviour. The optimally-blended rules lead to the convergence rates of order 2p + 2 as predicted by the theory of subsection 3.4. Figure 5 compares the relative approximation errors for the 3D biharmonic eigenvalue problem. In the multidimensional plots shown below, the axes correspond to eigenvalue indices j, l, q in (2.4) and (2.5). Again, using mixed isogeometric analysis with optimally-blended rules leads to smaller errors when compared to the fully-integrated case that employs standard Gauss quadratures. Tables 2 to 5 show the first, second, fourth, and eighth eigenvalue errors when using mixed isogeometric analysis with standard Gauss quadratures and optimally-blended rules for linear, quadratic, and cubic elements. In these tables, we denote by G when using the standard Gauss quadrature rule while by O when the optimally-blended rule is applied. There are different optimally-blended rules and they lead to the same numerical results; see [20,52]. Herein, we use the G p+1 -L p+1 optimally-blended rules. We also denote the convergence rates as ρ p for p-th order elements. For p-th order elements, we  obtain convergence rates of order 2p + 2 when using the optimally-blended rules. These numerical experiments verify our theoretical findings. Figure 6 shows the convergence of the eigenvalues λ j , j = 1, 2, 4, 8 of the 2D bi-14 harmonic problem. The squared eigenvalue errors of the quadratic mixed isogeometric analysis with standard Gauss quadratures have convergence order close to 4 (that is 2p). Using optimally-blended rules leads to smaller eigenvalue errors and convergence order of 6 (that is 2p + 2). Figure 7 shows the convergence of 2D Cahn-Hilliard operator eigenvalues when periodic boundary conditions are employed. Similar results are obtained for the phase-field equations.

Conclusions and future outlook
We present and study a mixed formulation of isogeometric analysis for a set of 2norder differential eigenvalue problems, which includes operators arising from the biharmonic problem, Cahn-Hilliard, Swift-Hohenberg, and phase-field crystal equations. We show that the mixed formulation with standard quadrature rules applied to integrate the inner-products leads to optimal rates (2p) of convergence on eigenvalue errors while the mixed formulation with optimally-blended rules leads to super-convergent 2p + 2 approximated eigenvalues. This work generalizes the theoretical results obtained for Laplacian eigenvalue problem to higher-order differential eigenvalue problems.
ture Science Platforms of the Commonwealth Scientific Industrial Research Organisation, CSIRO, of Australia. Additional support was provided by the European Union's Horizon