# A boundary penalization technique to remove outliers from isogeometric
analysis on tensor-product meshes

Quanling Deng  Victor M. Calo  Department of Mathematics, University of
Wisconsin–Madison, Madison, WI 53706, USA  School of Electrical Engineering,
Computing and Mathematical Sciences, Curtin University, Perth, WA 6845,
Australia  Mineral Resources, Commonwealth Scientific and Industrial Research
Organisation (CSIRO), Kensington, Perth, WA 6152, Australia

###### Abstract

We introduce a boundary penalization technique to improve the spectral
approximation of isogeometric analysis (IGA). The method removes the outliers
appearing in the high-frequency region of the approximate spectrum when using
the $C^{p-1},p$-th ($p\geq 3$) order isogeometric elements. We focus on the
classical Laplacian (Dirichlet) eigenvalue problem in 1D to illustrate the
idea and then use the tensor-product structure to generate the stiffness and
mass matrices for multidimensional problems. To remove the outliers, we
penalize the higher-order derivatives from both the solution and test spaces
at the domain boundary. Intuitively, we construct a better approximation by
weakly imposing features of the exact solution. Effectively, we add terms to
the variational formulation at the boundaries with minimal extra-computational
cost. We then generalize the idea to remove the outliers for the isogeometric
analysis of the Neumann eigenvalue problem (for $p\geq 2$). The boundary
penalization does not change the test and solution spaces. In the limiting
case, when the penalty goes to infinity, we perform the dispersion analysis of
$C^{2}$ cubic elements for the Dirichlet eigenvalue problem and $C^{1}$
quadratic elements for the Neumann eigenvalue problem. We obtain the
analytical eigenpairs for the resulting matrix eigenvalue problems. Numerical
experiments show optimal convergence rates for the eigenvalues and
eigenfunctions of the discrete operator.

###### keywords:

isogeometric analysis , spectral error , eigenvalue , outlier , boundary
penalization , dispersion analysis

††journal: CMAME

##  1 Introduction

Isogeometric analysis is a Galerkin finite element method, introduced in 2005
[1, 2]. The method is popular and efficient. The primary motivation for its
introduction was to unify the classical finite element analysis with computer-
aided design and analysis tools. Isogeometric analysis uses basis functions
with higher-order continuity (smoother), for example, the B-splines or non-
uniform rational basis splines (NURBS). This alternative set of basis
functions improves the numerical approximations of various partial
differential equations.

Cottrell et al. [3] studied with isogeometric analysis structural vibrations
(an eigenvalue problem) and compared the spectral approximation properties of
classical finite elements against isogeometric analysis for the first time.
The isogeometric elements improved the accuracy of the spectral approximation
significantly in the entire discrete spectrum. Maximum continuity basis
functions remove the optical branch from the discrete spectrum of the finite
element approximations. Nevertheless, a thin layer of outliers appears in the
high-frequency range of the spectrum. Since then, our community has been
seeking a method to remove the outliers.

There is rich literature in this line of work; see [4, 5, 6, 7, 8, 9, 10, 11,
12] and references therein. In [5], the authors further explored the
advantages of isogeometric analysis over finite elements on the spectral
approximation. They used the Pythagorean eigenvalue error theorem (see a proof
in [13]) to relate the eigenvalue and eigenfunction errors. A generalization
of this theorem that accounts for inexact quadrature in the inner products
associated with the stiffness and mass matrices was developed in [10] and this
generalized theorem was applied to develop the optimally-blended quadrature
rules. The work [14, 15] developed quadrature blending rules to reduce the
dispersion errors for finite element analysis of the Helmholtz equation. In
[7, 10], the authors generalized the idea to isogeometric elements using the
dispersion analysis and reduced the spectral errors for the differential
eigenvalue problems. The introduction of a single quadrature rule reduced the
computational cost for the blended rules [9] for $C^{1}$ quadratic
isogeometric elements. The paper [11] studied the error reductions, stopping
bands in the finite element numerical spectrum, and outliers in the
isogeometric spectrum with variable continuities. Also, error reductions for
other differential operators are possible [16, 17]. This sequence of results
delivers superconvergent eigenvalue errors. However, the outliers still
pollute the high-frequencies of the spectrum, particularly in multiple
dimensions.

In this paper, we propose a boundary penalization technique to remove the
outliers from the isogeometric spectrum. We consider the spectral
approximation of the Laplace operator. We first present the idea in 1D and
then generalize the construction to 2D and 3D using tensor-product
discretizations. We penalize high-order derivatives near the boundary to
remove the outliers. Intuitively, we add to the variational form near the
boundary terms that build the numerical method features of the exact solution.
These terms weakly impose the differential equation $-\Delta u=\lambda u$ and
its higher-order derivatives, $(-\Delta)^{\alpha}u=\lambda^{\alpha}u$, at the
domain boundaries. These terms assume high regularity for the eigenmodes. We
only use $(-\Delta)^{\alpha}u=\lambda^{\alpha}u$ for the degrees of freedom
near the boundary nodes. We impose these conditions weakly by penalization. We
show that this technique removes the outliers in the numerical spectrum
approximated by smooth isogeometric elements. As the penalty goes to infinity,
we impose the condition $(-\Delta)^{\alpha}u=\lambda^{\alpha}u$ more strongly.
In the limiting case, the penalty on $(-\Delta)^{\alpha}u=\lambda^{\alpha}u$
at the boundaries becomes an extra condition for the matrix system. For
$C^{2}$ cubic elements for the Dirichlet eigenvalue problem and $C^{1}$
quadratic elements for the Neumann eigenvalue problem, based on the work [18],
we perform the dispersion analysis (which unifies with the spectrum analysis
in [4]) near the boundary and give exact eigenpairs for the resulting matrix
problems. For higher-order elements, a reconstruction of the basis functions
near the boundaries could lead to exact eigenpairs of the matrix problems. In
[19] the authors proposed the strong imposition of higher-order derivatives
for 1D problems to remove outliers, although they did not discuss the
extension for multidimensional problems. However, this is not the focus of
this work; we will report these results in future work.

The structure of the paper is the following. Section 2 presents the problem
and its Galerkin discretization. In Section 3, we discuss how to remove the
outliers for isogeometric analysis. Section 4 seeks to remove the outliers
from the isogeometric analysis of the Neumann eigenvalue problem. In Section
5, we perform the dispersion analysis for the nodes near the boundaries and
give exact eigenpairs of the resulting matrix eigenvalue problems. Section 6
collects numerical examples that demonstrate the performance of the proposed
ideas. We present concluding remarks in Section 7.

##  2 Problem setting: Dirichlet eigenvalue problem

The classical second-order Dirichlet eigenvalue problem reads: find the
eigenpairs $(\lambda,u)\in\mathbb{R}^{+}\times H^{1}_{0}(\Omega)$ with
$\|u\|_{\Omega}=1$ such that

$\displaystyle-\Delta u$ $\displaystyle=\lambda u\quad\text{in}\quad\Omega,$
(2.1) $\displaystyle u$ $\displaystyle=0\quad\text{on}\quad\partial\Omega,$

where $\Delta=\nabla^{2}$ is the Laplacian,
$\Omega=[0,1]^{d}\subset\mathbb{R}^{d},d=1,2,3$, is a bounded open domain with
Lipschitz boundary. The eigenvalue problem (2.1) has a countable set of
eigenvalues $\lambda_{j}\in\mathbb{R}^{+}$ (see, for example, [20, Sec. 9.8])

$0<\lambda_{1}<\lambda_{2}\leq\lambda_{3}\leq\cdots$

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

$(u_{j},u_{k})=\delta_{jk},$ (2.2)

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

###  2.1 Galerkin discretization

To discretize (2.1) with finite and isogeometric elements, for simplicity, we
first assume that $\Omega$ is a cube and we use a uniform tensor product mesh
of size $h_{x}>0,h_{y}>0,h_{z}>0$ on $\Omega$. We denote each element as $E$
and its collection as $\mathcal{T}_{h}$ such that
$\bar{\Omega}=\cup_{E\in\mathcal{T}_{h}}E$. Let
$h=\max_{E\in\mathcal{T}_{h}}\text{diameter}(E)$. The variational formulation
of (2.1) at the continuous level is to find $\lambda\in\mathbb{R}^{+}$ and
$u\in H^{1}_{0}(\Omega)$ such that

$a(w,u)=\lambda b(w,u),\quad\forall\ w\in H^{1}_{0}(\Omega),$ (2.3)

where $a(w,v)=(\nabla w,\nabla v)$ and $b(w,v)=(w,v)$. Herein, we use the
standard notation for the Hilbert and Sobolev spaces. We denote by
$(\cdot,\cdot)$ and $\|\cdot\|$ the $L^{2}$-inner product and its norm,
respectively. From (2.2), the normalized eigenfunctions are also orthogonal in
the energy inner product

$a(u_{j},u_{k})=\lambda_{j}b(u_{j},u_{k})=\lambda_{j}\delta_{jk}.$ (2.4)

We specify a finite dimensional approximation space $V^{h}_{p}\subset
H^{1}_{0}(\Omega)$ where $V^{h}_{p}=\text{span}\\{\theta^{j}_{p}\\}$ is the
span of the B-spline basis functions $\theta^{j}_{p}$. The definition of the
B-spline basis functions in 1D is as follows. Let
$X=\\{x_{0},x_{1},\cdots,x_{m}\\}$ be a knot vector with knots $x_{j}$, which
is a non-decreasing sequence of real numbers. The $j$-th B-spline basis
function of degree $p$, denoted as $\theta^{j}_{p}(x)$, is defined as [21, 22]

$\displaystyle\theta^{j}_{0}(x)$ $\displaystyle=\begin{cases}1,\quad\text{if}\
x_{j}\leq x<x_{j+1}\\\ 0,\quad\text{otherwise}\\\ \end{cases}$ (2.5)
$\displaystyle\theta^{j}_{p}(x)$
$\displaystyle=\frac{x-x_{j}}{x_{j+p}-x_{j}}\theta^{j}_{p-1}(x)+\frac{x_{j+p+1}-x}{x_{j+p+1}-x_{j+1}}\theta^{j+1}_{p-1}(x).$

The span of these basis functions generates a finite-dimensional subspace of
the $H^{1}(\Omega)$ (see [23, 24] for details):

$V^{h}_{p}=\text{span}\\{\Theta_{j}^{p}\\}_{j=1}^{N_{h}}=\begin{cases}\text{span}\\{\theta^{j_{x}}_{p_{x}}(x)\\}_{j_{x}=1}^{N_{x}},&\text{in
1D},\\\
\text{span}\\{\theta^{j_{x}}_{p_{x}}(x)\theta^{j_{y}}_{p_{y}}(y)\\}_{j_{x},j_{y}=1}^{N_{x},N_{y}},&\text{in
2D},\\\
\text{span}\\{\theta^{j_{x}}_{p_{x}}(x)\theta^{j_{y}}_{p_{y}}(y)\theta^{j_{z}}_{p_{z}}(z)\\}_{j_{x},j_{y},j_{z}=1}^{N_{x},N_{y},N_{z}},&\text{in
3D},\\\ \end{cases}$ (2.6)

where $p_{x},p_{y},p_{z}$ specify the approximation order in each dimension,
respectively. $N_{x},N_{y},N_{z}$ are the total numbers of basis functions in
each dimension, and $N_{h}$ is the total number of degrees of freedom.

The isogeometric analysis of (2.1) seeks $\lambda^{h}\in\mathbb{R}^{+}$ and
$u^{h}\in V^{h}_{p}$ such that

$a(w^{h},u^{h})=\lambda^{h}b(w^{h},u^{h}),\quad\forall\ w^{h}\in V^{h}_{p}.$
(2.7)

We approximate the eigenfunctions as a linear combination of the B-spline
basis functions and substitute all the B-spline basis functions for $w_{h}$ in
(2.7). This approximation leads to the following matrix eigenvalue problem

$KU=\lambda^{h}MU,$ (2.8)

where
$K_{kl}=a(\theta_{p}^{k},\theta_{p}^{l}),M_{kl}=b(\theta_{p}^{k},\theta_{p}^{l}),$
and $U$ is the corresponding representation of the eigenvector as the
coefficients of the B-spline basis functions. When we restrict the continuity
of the B-splines to $C^{0}$, the method becomes a finite element with
Bernstein basis functions.

##  3 Boundary penalization technique to remove outliers

We consider $C^{p-1},p$-th order isogeometric elements. Spectral error
outliers start to appear when using $C^{2}$ cubic isogeometric elements for
the Dirichlet eigenvalue problem. For the Neumann eigenvalue problem, they
start to appear when using $C^{1}$ quadratic isogeometric elements. We first
describe the key idea, and then we present the idea using modified bilinear
forms in 1D.

The motivating idea stems from the use of high-order numerical time-marching
schemes for time-dependent problems. To illustrate the idea, we consider the
wave equation $u_{tt}(x,t)-\Delta u(x,t)=f(x,t)$ with an initial solution and
velocity. For a sufficiently regular solution $u$ at $x=0$, we can approximate
the initial acceleration as $\Delta u(x,0)+f(x,0)$. This consistency condition
allows us to compute the initial acceleration for the time integration
schemes, for example, the HHT-$\alpha$ [25] and generalized-$\alpha$ [26]
methods. Similarly, one can determine the initial higher-order time-derivative
values for the heat equation by taking time-derivatives when the solution is
regular enough. This class of consistency conditions is common for
initializing many higher-order time-marching methods.

Now, we borrow the idea for the eigenvalue problem (2.1). Consequently, we
expect that the function $-\Delta u$ also vanishes at the boundary. Figure 1
shows the $C^{2}$ cubic and $C^{3}$ quartic splines with their second
derivatives. Imposing the homogeneous Dirichlet boundary conditions removes
the first and last basis functions. However, there are still two basis
functions that have non-zero second derivatives at each boundary. Thus, we
modify the discrete system to impose the second derivatives’ homogeneity (zero
value) weakly.

Figure 1: $C^{2}$ cubic and $C^{3}$ quartic splines with their second
derivatives.

For $C^{2}$ cubic elements, we remove the outliers from the isogeometric
approximations by penalizing the Laplacian of the approximate solution at the
boundaries. For quintic $C^{4}$ and sextic $C^{5}$ isogeometric elements,
there are two more outliers in the spectrum. Applying a similar argument, we
obtain that the bi-Laplacian of the function also has to be zero on the
boundary. Thus, we remove these additional outliers by penalizing the higher-
order derivatives of the approximate solutions at the boundaries. This
procedure extends to higher-order elements by following a similar
construction. Herein, we remark that the extra regularity we require from the
solution in these weak constraints is of the same order as the one required by
operator splitting schemes that are typically used on tensor-product
discretizations to accelerate the computation (see, for example, [27, 28]).

Firstly, we denote

$\alpha=\lfloor\frac{p-1}{2}\rfloor=\begin{cases}\frac{p-1}{2},&p\quad\text{is
odd},\\\ \frac{p-2}{2},&p\quad\text{is even}.\\\ \end{cases}$ (3.1)

We present the idea in 1D. We obtain the 1D matrix problem and then apply the
tensor-product structure to construct an approximate matrix problem for
multidimensional problems.

In 1D, we consider the problem (2.1) on the unit interval $\Omega=[0,1]$. We
propose the boundary penalization to the isogeometric analysis of (2.1) in 1D
as: Find $\tilde{\lambda}^{h}\in\mathbb{R}$ and $\tilde{u}^{h}\in V^{h}_{p}$
such that

$\tilde{a}(w^{h},\tilde{u}^{h})=\tilde{\lambda}^{h}b(w^{h},\tilde{u}^{h}),\quad\forall\
w^{h}\in V^{h}_{p},$ (3.2)

where for $w,v\in V^{h}_{p}$

$\displaystyle\tilde{a}(w,v)$ $\displaystyle=\int_{0}^{1}w^{\prime}v^{\prime}\
dx+\sum_{\ell=1}^{\alpha}\eta_{a,\ell}\pi^{2}h^{6\ell-3}s^{(2\ell)}(w,v),$
(3.3a) $\displaystyle\tilde{b}(w,v)$ $\displaystyle=\int_{0}^{1}wv\
dx+\sum_{\ell=1}^{\alpha}\eta_{b,\ell}h^{6\ell-1}s^{(2\ell)}(w,v),$ (3.3b)
$\displaystyle s^{(2\ell)}(w,v)$
$\displaystyle=w^{(2\ell)}(0)v^{(2\ell)}(0)+w^{(2\ell)}(1)v^{(2\ell)}(1),$
(3.3c)

and $\eta_{a,\ell},\eta_{b,\ell}$ are parameters that penalize the (jumps of)
product of higher-order derivatives at the domain boundaries. As default we
set $\eta_{a,\ell}=\eta_{b,\ell}=1$. We refer to this penalty technique near
the boundaries as a discrete correction (DC) technique while the overall
method as a discretely corrected isogeometric analysis (DC-IGA).

######  Remark 1.

For the bilinear form $\tilde{a}(\cdot,\cdot)$, we set the scaling
$h^{6\ell-3}$, where $6\ell-3=(4\ell-1)+2(\ell-1)$. The part $(4\ell-1)$
contributes to the scaling that is consistent the scaling of the original
inner-product of the derivatives $(w^{\prime},v^{\prime})$. The part
$2(\ell-1)$ contributes to the consistency errors from the
$(-\Delta)^{\ell}u=\lambda^{\ell}u$. As a consequence, we rewrite
$\lambda^{\ell}=\lambda\lambda^{\ell-1}$, where $\lambda^{\ell-1}$ scales like
$\mathcal{O}(h^{2(1-\ell)})$ in the high-frequency region. Thus, the bilinear
form $\tilde{b}(\cdot,\cdot)$ has a scaling with two extra orders with respect
to $h$. We validate this observation in the numerical section.

As in (2.8), this leads to the matrix eigenvalue problem

$\tilde{K}^{1D}U=\tilde{\lambda}^{h}\tilde{M}^{1D}U,$ (3.4)

where
$\tilde{K}_{kl}=\tilde{a}(\theta_{p}^{k},\theta_{p}^{l}),\tilde{M}_{kl}=\tilde{b}(\theta_{p}^{k},\theta_{p}^{l}),$
and $U$ is the corresponding representation of the eigenvector as the
coefficients of the B-spline basis functions. Using the tensor-product
structure and introducing the outer product, the corresponding 2D matrix
problem is

$(\tilde{K}^{1D}_{x}\otimes\tilde{M}^{1D}_{y}+\tilde{M}^{1D}_{x}\otimes\tilde{K}^{1D}_{y})U=\tilde{\lambda}^{h}\tilde{M}^{1D}_{x}\otimes\tilde{M}^{1D}_{y}U,$
(3.5)

where $\otimes$ is the Kronecker product, and 3D matrix problem is

$(\tilde{K}^{1D}_{x}\otimes\tilde{M}^{1D}_{y}\otimes\tilde{M}^{1D}_{z}+\tilde{M}^{1D}_{x}\otimes\tilde{K}^{1D}_{y}\otimes\tilde{M}^{1D}_{z}+\tilde{M}^{1D}_{x}\otimes\tilde{M}^{1D}_{y}\otimes\tilde{K}^{1D}_{z})U=\tilde{\lambda}^{h}\tilde{M}^{1D}_{x}\otimes\tilde{M}^{1D}_{y}\otimes\tilde{M}^{1D}_{z}U,$
(3.6)

where $\tilde{K}^{1D}_{q},\tilde{M}^{1D}_{q},q=x,y,z$ are 1D matrices
generated as in (3.4) from the modified bilinear forms in (3.3). We refer to
[15, 29] for similar extensions to multiple dimensions.

Alternatively, for multidimensional problem (2.1), to remove the outliers, one
can penalize the boundary nodes similarly with the following bilinear forms

$\displaystyle\tilde{a}(w,v)$ $\displaystyle=(\nabla w,\nabla
v)+\sum_{\ell=1}^{\alpha}\eta_{a,\ell}\pi^{2}h^{6\ell-3}s^{(2\ell)}(w,v),$
(3.7a) $\displaystyle\tilde{b}(w,v)$
$\displaystyle=(w,v)+\sum_{\ell=1}^{\alpha}\eta_{b,\ell}h^{6\ell-1}s^{(2\ell)}(w,v),$
(3.7b) $\displaystyle s^{(2\ell)}(w,v)$
$\displaystyle=(\Delta^{\ell}w,\Delta^{\ell}v)_{\partial\Omega},$ (3.7c)

where $(\cdot,\cdot)_{\partial\Omega}$ denotes the $L^{2}$ inner product on
the boundary ${\partial\Omega}$. These bilinear forms lead to matrix
eigenvalue problems which are slightly different (i.e., modified entries near
the boundaries) from the problems (3.5) in 2D and (3.6) in 3D. The advantage
is that on the penalty terms, one does not require a tensor-product mesh (but
may require on the original isogeometric basis functions). In the numerical
simulation, we observe that these bilinear forms remove the outliers. However,
the spectral errors at high-frequency region are larger. Thus, for
multidimensional problem (2.1), we focus on the (3.5) in 2D and (3.6) in 3D.

######  Remark 2.

The new scheme corresponds to a consistent boundary penalization. The
correction is consistent in the sense that the extra terms vanish for smooth
solutions. For $p=1,2$, $\alpha=0$, thus, the extra terms become zeros. There
are no outliers and hence no modification to the bilinear forms. The
additional terms penalize the jumps of the spatial operators at the boundary
nodes. The jump terms consider the differences of the $2\ell$-th derivatives
of the basis functions and the Dirichlet boundary values (zeros, in this case,
hence omitted in the modified bilinear forms).

Our numerical results show that for $\ell\geq 1$, we can set $\eta_{a,\ell}=0$
or $\eta_{b,\ell}=0$ and assign the other parameter appropriately to remove
the outliers from the spectrum. However, the eigenfunction errors remain
large. When $\eta_{a,\ell}>0$ or $\eta_{b,\ell}>0$, the eigenfunction errors
decrease in the high-frequency region. By default, we set
$\eta_{a,\ell}=\eta_{b,\ell}=1$ for all $\ell$. The modifications to the
variational forms are local to the basis functions with support on the
elements contiguous to the boundary. With these modifications, we observe
optimally convergent errors in both eigenvalues and eigenfunctions.

##  4 Neumann eigenvalue problem

Now, we consider the Neumann eigenvalue problem

$\displaystyle-\Delta u$ $\displaystyle=\lambda u\quad\text{in}\quad\Omega,$
(4.1) $\displaystyle\nabla u\cdot\boldsymbol{n}$
$\displaystyle=0\quad\text{on}\quad\partial\Omega.$

For this problem, the outliers appear in the spectrum of isogeometric elements
when using $C^{1}$ quadratic elements. We describe the bilinear forms for the
1D case to obtain the matrix problem (3.4) that remove the outliers from the
discrete solution.

We denote

$\beta=\lfloor\frac{p}{2}\rfloor=\begin{cases}\frac{p}{2},&p\quad\text{is
even},\\\ \frac{p-1}{2},&p\quad\text{is odd}.\\\ \end{cases}$ (4.2)

Extending the arguments we propose for the Dirichlet eigenvalue problem (2.1),
we introduce a boundary penalization to the isogeometric discretization of
(4.1) in 1D as: Find $\hat{\lambda}^{h}\in\mathbb{R}$ and $\hat{u}^{h}\in
V^{h}_{p}$ such that

$\hat{a}(w^{h},\hat{u}^{h})=\hat{\lambda}^{h}b(w^{h},\hat{u}^{h}),\quad\forall\
w^{h}\in V^{h}_{p},$ (4.3)

where for $w,v\in V^{h}_{p}$

$\displaystyle\hat{a}(w,v)$ $\displaystyle=\int_{0}^{1}w^{\prime}v^{\prime}\
dx+\sum_{\ell=1}^{\beta}\hat{\eta}_{a,\ell}\pi^{2}h^{6\ell-5}s^{(2\ell-1)}(w,v),$
(4.4a) $\displaystyle\hat{b}(w,v)$ $\displaystyle=\int_{0}^{1}wv\
dx+\sum_{\ell=1}^{\beta}\hat{\eta}_{b,\ell}h^{6\ell-3}s^{(2\ell-1)}(w,v),$
(4.4b) $\displaystyle s^{(2\ell-1)}(w,v)$
$\displaystyle=w^{(2\ell-1)}(0)v^{(2\ell-1)}(0)+w^{(2\ell-1)}(1)v^{(2\ell-1)}(1),$
(4.4c)

and $\hat{\eta}_{a,\ell},\hat{\eta}_{b,\ell}$ penalize the jumps. By default,
we set $\hat{\eta}_{a,\ell}=\hat{\eta}_{b,\ell}=1$. We use the scaling orders
by following the arguments in Remark 1 for the Dirichlet eigenvalue problem.
Similarly, for multidimensional problem (4.1), we use (3.5) and (3.6) for 2D
and 3D matrix problems, respectively.

##  5 Infinite penalization and boundary dispersion analysis

In this section, we consider the limiting case when the penalty goes to
$\infty$. In such a case, the penalized terms on
$(-\Delta)^{\alpha}u=\lambda^{\alpha}u$ at the boundaries become strong
boundary conditions for the resulting matrix problems. In particular, for the
Dirichlet eigenvalue problem (2.1), we present the dispersion analysis near
the boundaries for $C^{2}$ cubic isogeometric elements. For the Neumann
eigenvalue problem (4.1), we present the dispersion analysis near the
boundaries for $C^{1}$ quadratic isogeometric elements. In these two cases on
uniform meshes, we obtain the exact eigenvalues and eigenvectors from the
resulting matrix eigenvalue problems. To give the exact eigenpairs, we first
present the following result.

######  Theorem 1 (Analytical eigenvalues and eigenvectors, case 1).

Let $G^{(\xi)}=(G^{(\xi)}_{j,k})$ be a square matrix with dimension $n$ such
that

$G^{(\xi)}_{j,j+k}=\xi_{|k|},\quad|k|\leq m,\quad j=1,\cdots,n.$ (5.1)

Let $H^{(\xi)}=(H^{(\xi)}_{j,k})$ be a zero square matrix with dimension $n$.
For $m\geq 2$, we modify $H^{(\xi)}$ with

$H^{(\xi)}_{j,k}=H^{(\xi)}_{n-j+1,n-k+1}=\xi_{j+k},\quad k=1,\cdots,m-j,\quad
j=1,\cdots,m-1.$ (5.2)

Let $A=G^{(\mu)}-H^{(\mu)},B=G^{(\nu)}-H^{(\nu)}.$ Assume that $B$ is
invertible. Then, the generalized matrix eigenvalue problem $AX=\lambda BX$
has eigenpairs $(\lambda_{j},X_{j})$ with $X_{j}=(X_{j,1},\cdots,X_{j,n})^{T}$
where

$\lambda_{j}=\frac{\mu_{0}+2\sum_{l=1}^{m}\mu_{l}\cos(lj\pi
h)}{\nu_{0}+2\sum_{l=1}^{m}\nu_{l}\cos(lj\pi h)},\quad X_{j,k}=c\sin(j\pi
kh),\quad h=\frac{1}{n+1},\ j,k=1,2,\cdots,n,$ (5.3)

where $0\neq c\in\mathbb{R}.$

###### Proof.

The matrices defined here are Toeplitz matrices with modified entries near the
boundaries. From Bloch wave assumption and with the modifications near the
boundaries in mind, we seeks eigenvectors of the form $\sin(j\pi kh)$. Using
the trigonometric identity
$\sin(\alpha\pm\beta)=\sin(\alpha)\cos(\beta)\pm\cos(\alpha)\sin(\beta)$, one
can verify that each row of the problem $AX=\lambda BX$, i.e.,
$\sum_{k=1}^{n}A_{ik}X_{j,k}=\lambda\sum_{k=1}^{n}B_{ik}X_{j,k},i=1,\cdots,n$,
reduces to $\alpha_{0}+2\sum_{l=1}^{m}\alpha_{l}\cos(lj\pi
h)=\lambda\big{(}\beta_{0}+2\sum_{l=1}^{m}\beta_{l}\cos(lj\pi h)\big{)}$. This
is independent of the row number $i$, that is, all the boundary and internal
rows lead to the same expression. The problem has at most $n$ eigenpairs and
the $n$ eigenvectors are linearly independent. Thus, the $n$ eigenpairs
$(\lambda_{j},X_{j})$ given in (5.3) are the analytical solutions to the
problem $AX=\lambda BX$. This completes the proof. ∎

The eigenpairs given in this theorem correspond to the Dirichlet eigenvalue
problem. Similarly, we have the following result for the Neumann eigenvalue
problem. We omit the proof for brevity.

######  Theorem 2 (Analytical eigenvalues and eigenvectors, case 2).

Let $G^{(\xi)}$ be defined in (5.1). Let $H^{(\xi)}=(H^{(\xi)}_{j,k})$ be a
zero square matrix with dimension $n$. We modify $H^{(\xi)}$ with

$H^{(\xi)}_{j,k}=H^{(\xi)}_{n-j+1,n-k+1}=\xi_{j+k-1},\quad
k=1,\cdots,m-j+1,\quad j=1,\cdots,m.$ (5.4)

Let $A=G^{(\mu)}+H^{(\mu)},B=G^{(\nu)}+H^{(\nu)}.$ Assume that $B$ is
invertible. Let $h=\frac{1}{n}$. Then, the generalized matrix eigenvalue
problem $AX=\lambda BX$ has eigenpairs $(\lambda_{j},X_{j})$ with
$X_{j}=(X_{j,1},\cdots,X_{j,n})^{T}$ where

$\lambda_{j+1}=\frac{\mu_{0}+2\sum_{l=1}^{m}\mu_{l}\cos(lj\pi
h)}{\nu_{0}+2\sum_{l=1}^{m}\nu_{l}\cos(lj\pi h)},\qquad
X_{j+1,k}=c\cos\big{(}j\pi(k-\frac{1}{2})h\big{)}$ (5.5)

with $k=1,2,\cdots,n,\ j=0,1,\cdots,n-1,0\neq c\in\mathbb{R}.$

The matrix $G$ in (5.1) is a Toeplitz matrix while the matrices $H$ defined in
(5.2) and (5.4) are zero matrices with modified entries at the first and last
few rows (near boundaries). All these matrices are symmetric and persymmetric.
We refer to [18, Section 2] for more details.

###  5.1 $C^{2}$ cubic elements for the Dirichlet eigenvalue problem

First, we solve the problem (2.1) in 1D by the isogeometric elements (2.7)
with $C^{2}$ cubic B-splines. We discretize the unit interval $\Omega=[0,1]$
with a uniform grid $\mathcal{T}_{h}=\\{x_{j}=jh\\}_{j=0}^{N}$, where $h=1/N$
is the grid size. The $C^{2}$ cubic isogeometric elements lead to the
stiffness and mass matrices (after applying the homogeneous Dirichlet boundary
conditions)

$\displaystyle K$
$\displaystyle=\frac{1}{h}\begin{bmatrix}\frac{3}{2}&\frac{3}{80}&-\frac{1}{4}&-\frac{1}{80}\\\\[5.69046pt]
\frac{3}{80}&\frac{27}{40}&-\frac{1}{30}&-\frac{47}{240}&-\frac{1}{120}\\\\[5.69046pt]
-\frac{1}{4}&-\frac{1}{30}&\frac{2}{3}&-\frac{1}{8}&-\frac{1}{5}&-\frac{1}{120}\\\\[5.69046pt]
-\frac{1}{80}&-\frac{47}{240}&-\frac{1}{8}&\frac{2}{3}&-\frac{1}{8}&-\frac{1}{5}&-\frac{1}{120}\\\\[5.69046pt]
&-\frac{1}{120}&-\frac{1}{5}&-\frac{1}{8}&\frac{2}{3}&-\frac{1}{8}&-\frac{1}{5}&-\frac{1}{120}\\\\[5.69046pt]
&&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\\\[5.69046pt]
&&&-\frac{1}{120}&-\frac{1}{5}&-\frac{1}{8}&\frac{2}{3}&-\frac{1}{8}&-\frac{47}{240}&-\frac{1}{80}\\\\[5.69046pt]
&&&&-\frac{1}{120}&-\frac{1}{5}&-\frac{1}{8}&\frac{2}{3}&-\frac{1}{30}&-\frac{1}{4}\\\\[5.69046pt]
&&&&&-\frac{1}{120}&-\frac{47}{240}&-\frac{1}{30}&\frac{27}{40}&\frac{3}{80}\\\\[5.69046pt]
&&&&&&-\frac{1}{80}&-\frac{1}{4}&\frac{3}{80}&\frac{3}{2}\\\\[5.69046pt]
\end{bmatrix},$ (5.6) $\displaystyle M$
$\displaystyle=h\begin{bmatrix}\frac{31}{140}&\frac{5}{32}&\frac{29}{840}&\frac{1}{3360}\\\\[5.69046pt]
\frac{5}{32}&\frac{183}{560}&\frac{283}{1260}&\frac{239}{10080}&\frac{1}{5040}\\\\[5.69046pt]
\frac{29}{840}&\frac{283}{1260}&\frac{151}{315}&\frac{397}{1680}&\frac{1}{42}&\frac{1}{5040}\\\\[5.69046pt]
\frac{1}{3360}&\frac{239}{10080}&\frac{397}{1680}&\frac{151}{315}&\frac{397}{1680}&\frac{1}{42}&\frac{1}{5040}\\\\[5.69046pt]
&\frac{1}{5040}&\frac{1}{42}&\frac{397}{1680}&\frac{151}{315}&\frac{397}{1680}&\frac{1}{42}&\frac{1}{5040}\\\\[5.69046pt]
&&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\\\[5.69046pt]
&&&\frac{1}{5040}&\frac{1}{42}&\frac{397}{1680}&\frac{151}{315}&\frac{397}{1680}&\frac{239}{10080}&\frac{1}{3360}\\\\[5.69046pt]
&&&&\frac{1}{5040}&\frac{1}{42}&\frac{397}{1680}&\frac{151}{315}&\frac{283}{1260}&\frac{29}{840}\\\\[5.69046pt]
&&&&&\frac{1}{5040}&\frac{239}{10080}&\frac{283}{1260}&\frac{183}{560}&\frac{5}{32}\\\\[5.69046pt]
&&&&&&\frac{1}{3360}&\frac{29}{840}&\frac{5}{32}&\frac{31}{140}\\\\[5.69046pt]
\end{bmatrix},$

which are of dimension $(N+1)\times(N+1)$.

For the purpose of dispersion analysis, one assumes that the component of the
eigenvector $U_{j}$ of (2.8) takes the form (see, e.g., [3, Sec. 5] and [4,
Sec. 4])

$U_{j,k}=\sin(\omega_{j}kh),\quad k=1,2,\cdots,N+1,$ (5.7)

where $\omega_{j}=j\pi$ is the eigenfrequency associated with the $j$-th mode.
For an internal node with $5\leq k\leq N-3$, the $k$-th row of the matrix
problem (2.8) with $K$ and $M$ defined in (5.6) leads to the dispersion
relation

$\lambda^{h}_{j}h^{2}=-42+\frac{1008\big{(}52+49\cos(\omega_{j}h)+4\cos(2\omega_{j}h)\big{)}}{1208+1191\cos(\omega_{j}h)+120\cos(2\omega_{j}h)+\cos(3\omega_{j}h)}.$
(5.8)

With $\omega_{j}h<1$, a Taylor expansion of (5.8) leads to the dispersion
error representation

$\Lambda^{h}=\Lambda+\frac{1}{30240}\Lambda^{4}+\mathcal{O}(\Lambda^{5}),$
(5.9)

where $\Lambda=(\omega_{j}h)^{2}$ and $\Lambda^{h}=h^{2}\lambda^{h}_{j}$. The
dispersion error is given by subtracting both sides $\Lambda$ and dividing
both sides by $h^{2}$. Similarly, the first four (similarly, for the last
four) rows of matrix problem (2.8) with $K$ and $M$ defined in (5.6) lead to
the dispersion relations

$\displaystyle\lambda^{h}_{j}h^{2}$
$\displaystyle=-42+\frac{2016\big{(}10+11\cos(\omega_{j}h)+2\cos(2\omega_{j}h)\big{)}}{430+526\cos(\omega_{j}h)+116\cos(2\omega_{j}h)+\cos(3\omega_{j}h)},$
(5.10) $\displaystyle\lambda^{h}_{j}h^{2}$
$\displaystyle=-42+\frac{4032\big{(}40+76\cos(\omega_{j}h)+47\cos(2\omega_{j}h)+4\cos(3\omega_{j}h)\big{)}}{3841+7066\cos(\omega_{j}h)+4532\cos(2\omega_{j}h)+478\cos(3\omega_{j}h)+4\cos(4\omega_{j}h)},$
$\displaystyle\lambda^{h}_{j}h^{2}$
$\displaystyle=-42+\frac{1008\big{(}57+96\cos(\omega_{j}h)+108\cos(2\omega_{j}h)+\cdots\big{)}}{1355+2324\cos(\omega_{j}h)+2536\cos(2\omega_{j}h)+\cdots},$
$\displaystyle\lambda^{h}_{j}h^{2}$
$\displaystyle=-42+\frac{84\big{(}96\sin(2\omega_{j}h)+1176\sin(3\omega_{j}h)+\cdots\big{)}}{3\sin(\omega_{j}h)+239\sin(2\omega_{j}h)+\cdots},$

which have the following error representations

$\displaystyle\Lambda^{h}$
$\displaystyle=\frac{1302}{1073}+\mathcal{O}(\Lambda),$ (5.11)
$\displaystyle\Lambda^{h}$
$\displaystyle=\frac{518}{1769}+\mathcal{O}(\Lambda),$
$\displaystyle\Lambda^{h}$
$\displaystyle=\frac{42}{941}+\mathcal{O}(\Lambda),$
$\displaystyle\Lambda^{h}$
$\displaystyle=\frac{14}{13439}+\mathcal{O}(\Lambda),$

respectively. These dispersion relations and errors near the boundary are
inconsistent with (5.8) and (5.9) for the internal nodes. This inconsistency
leads to the outliers in the approximate spectrum. The boundary penalization
technique proposed in Section 3 recovers the consistency weakly. In the
limiting case, when the penalty goes to infinity for $C^{2}$ cubic elements,
we recover the consistency strongly. We perform the analysis as follows.

In the limiting case when the penalty goes to $\infty$, the penalized terms in
(3.3) at the boundaries $x=0$ and $x=1$ become the conditions on the unknown
vector components of the generalized matrix eigenvalue problem (2.8)

$3U_{1}-U_{2}=0,\qquad 3U_{N+1}-U_{N}=0.$ (5.12)

We apply these conditions to reduce the matrix eigenvalue problem (2.8) from
dimension $(N+1)\times(N+1)$ to dimension $(N-1)\times(N-1)$. In particular,
we multiply the first column by $\frac{1}{3}$ and add the result to the second
column, then multiply the first row by $\frac{1}{3}$ and add the result to the
second row. Similarly, we multiply the $(N+1)$-th column by $\frac{1}{3}$ and
add the result to the $N$-th column, then multiply the $(N+1)$-th row by
$\frac{1}{3}$ and add the result to the $N$-th row. Removing the first and
last rows and columns, we obtain

$\displaystyle\tilde{K}$
$\displaystyle=\frac{1}{h}\begin{bmatrix}\frac{13}{15}&-\frac{7}{60}&-\frac{1}{5}&-\frac{1}{120}\\\\[5.69046pt]
-\frac{7}{60}&\frac{2}{3}&-\frac{1}{8}&-\frac{1}{5}&-\frac{1}{120}\\\\[5.69046pt]
-\frac{1}{5}&-\frac{1}{8}&\frac{2}{3}&-\frac{1}{8}&-\frac{1}{5}&-\frac{1}{120}\\\\[5.69046pt]
-\frac{1}{120}&-\frac{1}{5}&-\frac{1}{8}&\frac{2}{3}&-\frac{1}{8}&-\frac{1}{5}&-\frac{1}{120}\\\\[5.69046pt]
&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\\\[5.69046pt]
&&-\frac{1}{120}&-\frac{1}{5}&-\frac{1}{8}&\frac{2}{3}&-\frac{1}{8}&-\frac{1}{5}\\\\[5.69046pt]
&&&-\frac{1}{120}&-\frac{1}{5}&-\frac{1}{8}&\frac{2}{3}&-\frac{7}{60}\\\\[5.69046pt]
&&&&-\frac{1}{120}&-\frac{1}{5}&-\frac{7}{60}&\frac{13}{15}\\\\[5.69046pt]
\end{bmatrix},$ (5.13) $\displaystyle\tilde{M}$
$\displaystyle=h\begin{bmatrix}\frac{41}{90}&\frac{17}{72}&\frac{1}{42}&\frac{1}{5040}\\\\[5.69046pt]
\frac{17}{72}&\frac{151}{315}&\frac{397}{1680}&\frac{1}{42}&\frac{1}{5040}\\\\[5.69046pt]
\frac{1}{42}&\frac{397}{1680}&\frac{151}{315}&\frac{397}{1680}&\frac{1}{42}&\frac{1}{5040}\\\\[5.69046pt]
\frac{1}{5040}&\frac{1}{42}&\frac{397}{1680}&\frac{151}{315}&\frac{397}{1680}&\frac{1}{42}&\frac{1}{5040}\\\\[5.69046pt]
&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\\\[5.69046pt]
&&\frac{1}{5040}&\frac{1}{42}&\frac{397}{1680}&\frac{151}{315}&\frac{397}{1680}&\frac{1}{42}\\\\[5.69046pt]
&&&\frac{1}{5040}&\frac{1}{42}&\frac{397}{1680}&\frac{151}{315}&\frac{17}{72}\\\\[5.69046pt]
&&&&\frac{1}{5040}&\frac{1}{42}&\frac{17}{72}&\frac{41}{90}\\\\[5.69046pt]
\end{bmatrix},$

which are of dimensions $(N-1)\times(N-1)$. Applying Theorem 1, we have the
following result.

######  Theorem 3 (Analytical solutions).

Let $\tilde{K}$ and $\tilde{M}$ be defined in (5.13). The problem
$\tilde{K}U=\lambda^{h}\tilde{M}U$ has analytical eigenpairs
$(\lambda_{j}^{h},U_{j})$ where

$\lambda^{h}_{j}=-42N^{2}+\frac{1008N^{2}\big{(}52+49\cos(j\pi h)+4\cos(2j\pi
h)\big{)}}{1208+1191\cos(j\pi h)+120\cos(2j\pi h)+\cos(3j\pi h)},\quad
U_{j,k}=c\sin(j\pi kh)$ (5.14)

with $h=\frac{1}{N},j,k=1,2,\cdots,N-1,0\neq c\in\mathbb{R}$.

######  Remark 3.

For the matrix eigenvalue problem with matrices $\tilde{K}$ and $\tilde{M}$ as
defined in (5.13), the corresponding dispersion relations for both the
internal and boundary nodes are the same as in (5.8). That is, the dispersion
relations and errors are consistent for all the eigenvector components. In
such a case, the analytical approximate eigenvalues are given in (5.14), and
it is a smooth function in terms of the index $j$; therefore, no outliers. For
multidimensional problems, the analytical eigenpairs can be obtained by
applying Theorem 2.8 in [18]. For higher-order elements, one requires a
reconstruction of the basis functions near the boundaries to recover this
consistency. We will report these results shortly.

###  5.2 $C^{1}$ quadratic elements for the Neumann eigenvalue problem

We now present the boundary dispersion analysis for the Neumann eigenvalue
problem (4.1) when using $C^{1}$ quadratic elements. On a uniform grid with
$N$ elements in 1D, the isogeometric analysis with $C^{1}$ quadratic elements
for (4.1) leads to the following stiffness and mass matrices

$\displaystyle K$
$\displaystyle=\frac{1}{h}\begin{bmatrix}\frac{4}{3}&-1&-\frac{1}{3}\\\\[5.69046pt]
-1&\frac{4}{3}&-\frac{1}{6}&-\frac{1}{6}\\\\[5.69046pt]
-\frac{1}{3}&-\frac{1}{6}&1&-\frac{1}{3}&-\frac{1}{6}\\\\[5.69046pt]
&-\frac{1}{6}&-\frac{1}{3}&1&-\frac{1}{3}&-\frac{1}{6}\\\\[5.69046pt]
&&\ddots&\ddots&\ddots&\ddots&\ddots&\\\\[5.69046pt]
&&&-\frac{1}{6}&-\frac{1}{3}&1&-\frac{1}{3}&-\frac{1}{6}\\\\[5.69046pt]
&&&&-\frac{1}{6}&-\frac{1}{3}&1&-\frac{1}{6}&-\frac{1}{3}\\\\[5.69046pt]
&&&&&-\frac{1}{6}&-\frac{1}{6}&\frac{4}{3}&-1\\\\[5.69046pt]
&&&&&&-\frac{1}{3}&-1&\frac{4}{3}\\\ \end{bmatrix},$ (5.15) $\displaystyle M$
$\displaystyle=h\begin{bmatrix}\frac{1}{5}&\frac{7}{60}&\frac{1}{60}\\\\[5.69046pt]
\frac{7}{60}&\frac{1}{3}&\frac{5}{24}&\frac{1}{120}\\\\[5.69046pt]
\frac{1}{60}&\frac{5}{24}&\frac{11}{20}&\frac{13}{60}&\frac{1}{120}\\\\[5.69046pt]
&\frac{1}{120}&\frac{13}{60}&\frac{11}{20}&\frac{13}{60}&\frac{1}{120}\\\\[5.69046pt]
&&\ddots&\ddots&\ddots&\ddots&\ddots&\\\\[5.69046pt]
&&&\frac{1}{120}&\frac{13}{60}&\frac{11}{20}&\frac{13}{60}&\frac{1}{120}\\\\[5.69046pt]
&&&&\frac{1}{120}&\frac{13}{60}&\frac{11}{20}&\frac{5}{24}&\frac{1}{60}\\\\[5.69046pt]
&&&&&\frac{1}{120}&\frac{5}{24}&\frac{1}{3}&\frac{7}{60}\\\\[5.69046pt]
&&&&&&\frac{1}{60}&\frac{7}{60}&\frac{1}{5}\end{bmatrix},$

which are of dimensions $(N+2)\times(N+2)$.

Similarly, we assume that the internal component of the eigenvector $U_{j}$ of
(2.8) with these matrices takes the form

$U_{j,k}=\cos(\omega_{j}\big{(}k-\frac{3}{2})h\big{)},\quad k=1,2,\cdots,N+2,$
(5.16)

where $\omega_{j}=j\pi$ is the eigenfrequency associated with the $j$-th mode.
For an internal node with $4\leq k\leq N-1$, the $k$-th row of the matrix
problem (2.8) with $K$ and $M$ defined in (5.15) leads to the dispersion
relation

$\lambda^{h}_{j}h^{2}=-20+\frac{240\big{(}3+2\cos(\omega_{j}h)\big{)}}{33+26\cos(\omega_{j}h)+\cos(2\omega_{j}h)}.$
(5.17)

With $\omega_{j}h<1$, a Taylor expansion of (5.17) leads to the dispersion
error representation

$\Lambda^{h}=\Lambda+\frac{1}{720}\Lambda^{3}+\mathcal{O}(\Lambda^{4}).$
(5.18)

Similarly, the first two (similarly, for the last two) rows of matrix problem
(2.8) with $K$ and $M$ defined in (5.15) lead to the dispersion relations

$\displaystyle\lambda^{h}_{j}h^{2}$
$\displaystyle=\frac{40\sin^{2}(\omega_{j}h)}{15+24\cos(\omega_{j}h)+\cos(2\omega_{j}h)},$
(5.19) $\displaystyle\lambda^{h}_{j}h^{2}$
$\displaystyle=-20+\frac{200}{9+\cos(\omega_{j}h)},$

which have the error representations

$\displaystyle\Lambda^{h}=\Lambda+\frac{1}{60}\Lambda^{2}+\mathcal{O}(\Lambda^{3}),\qquad\Lambda^{h}=\Lambda-\frac{1}{30}\Lambda^{2}+\mathcal{O}(\Lambda^{3}),$
(5.20)

respectively. These dispersion relations and errors near the boundary are
inconsistent with (5.17) and (5.18) for the internal nodes. This inconsistency
leads to the outliers in the approximate spectrum. The boundary penalization
technique proposed in Section 4 recovers the consistency weakly. In the
limiting case, when the penalty goes to infinity for $C^{1}$ quadratic
elements, we recover consistency strongly. In the limiting case, when the
penalty goes to infinity, the penalized terms in (4.4) at the boundaries $x=0$
and $x=1$ become the conditions on the unknown vector components

$U_{1}-U_{2}=0,\qquad U_{N+2}-U_{N+1}=0.$ (5.21)

We apply these conditions to reduce the matrix eigenvalue problem (2.8) from
dimension $(N+2)\times(N+2)$ to dimension $N\times N$. In particular, we add
the first column to the second column, then add the first row to the second
row. Similarly, we add the $(N+2)$-th column to the $(N+1)$-th column, then
add the $(N+2)$-th row to the $(N+1)$-th row. Removing the first and last rows
and columns, we obtain

$\displaystyle\tilde{K}$
$\displaystyle=\frac{1}{h}\begin{bmatrix}\frac{2}{3}&-\frac{1}{2}&-\frac{1}{6}\\\\[5.69046pt]
-\frac{1}{2}&1&-\frac{1}{3}&-\frac{1}{6}\\\\[5.69046pt]
-\frac{1}{6}&-\frac{1}{3}&1&-\frac{1}{3}&-\frac{1}{6}\\\\[5.69046pt]
&\ddots&\ddots&\ddots&\ddots&\ddots&\\\\[5.69046pt]
&&-\frac{1}{6}&-\frac{1}{3}&1&-\frac{1}{3}&-\frac{1}{6}\\\\[5.69046pt]
&&&-\frac{1}{6}&-\frac{1}{3}&1&-\frac{1}{2}\\\\[5.69046pt]
&&&&-\frac{1}{6}&-\frac{1}{2}&\frac{2}{3}\\\ \end{bmatrix},$ (5.22)
$\displaystyle\tilde{M}$
$\displaystyle=h\begin{bmatrix}\frac{23}{30}&\frac{9}{40}&\frac{1}{120}\\\\[5.69046pt]
\frac{9}{40}&\frac{11}{20}&\frac{13}{60}&\frac{1}{120}\\\\[5.69046pt]
\frac{1}{120}&\frac{13}{60}&\frac{11}{20}&\frac{13}{60}&\frac{1}{120}\\\\[5.69046pt]
&\ddots&\ddots&\ddots&\ddots&\ddots&\\\\[5.69046pt]
&&\frac{1}{120}&\frac{13}{60}&\frac{11}{20}&\frac{13}{60}&\frac{1}{120}\\\\[5.69046pt]
&&&\frac{1}{120}&\frac{13}{60}&\frac{11}{20}&\frac{9}{40}\\\\[5.69046pt]
&&&&\frac{1}{120}&\frac{9}{40}&\frac{23}{30}\\\ \end{bmatrix},$

which are of dimensions $N\times N$. Applying Theorem 2, we have the following
result.

######  Theorem 4 (Analytical solutions).

Let $\tilde{K}$ and $\tilde{M}$ be defined in (5.22). The problem
$\tilde{K}U=\lambda^{h}\tilde{M}U$ has analytical eigenpairs
$(\lambda_{j}^{h},U_{j})$ where

$\lambda^{h}_{j+1}=-20N^{2}+\frac{240N^{2}\big{(}3+2\cos(j\pi
h)\big{)}}{33+26\cos(j\pi h)+\cos(2j\pi h)},\quad
U_{j,k}=c\cos\big{(}j\pi(k-\frac{1}{2})h\big{)}$ (5.23)

with $h=\frac{1}{N},k=1,2,\cdots,N,j=0,1,2,\cdots,N-1,0\neq c\in\mathbb{R}$.

Similarly, as for the $C^{2}$ cubic elements for the Dirichlet eigenvalue
problem, the corresponding dispersion relations for both the internal and
boundary nodes are consistent, and outliers are eliminated. These analytical
eigenvector components lead to optimal eigenfunction errors.

##  6 Numerical examples

In this section, we present numerical simulations of problem (2.1) in 1D, 2D,
and 3D using isogeometric elements. Once we solve the eigenvalue problem, we
sort the discrete eigenvalues in ascending order and pair them with the exact
eigenvalues. We focus on the numerical approximation properties of the
eigenvalues. In the case of 1D, however, we also report the eigenfunction (EF)
errors in $H^{1}$ semi-norm (energy norm). To have comparable scales, we
collect the relative eigenvalue errors and scale the energy norm by the
corresponding eigenvalues [10, 7]. The relative eigenvalue errors are

$e_{\lambda_{j}}=\frac{\lambda^{h}_{j}-\lambda_{j}}{\lambda_{j}},\qquad\tilde{e}_{\lambda_{j}}=\frac{\tilde{\lambda}^{h}_{j}-\lambda_{j}}{\lambda_{j}}.$
(6.1)

From the minimax principle [13], $e_{\lambda_{j}}>0$. To see the trend after
removing the outliers, we do not apply the absolute values of the eigenvalue
errors. We define the scaled eigenfunction error in the energy norm as

$e_{u_{j}}=\frac{|u_{j}-u_{j}^{h}|_{1,\Omega}}{\lambda_{j}},\qquad\tilde{e}_{u_{j}}=\frac{|u_{j}-\tilde{u}_{j}^{h}|_{1,\Omega}}{\lambda_{j}}.$
(6.2)

###  6.1 Dirichlet eigenvalue problem

####  6.1.1 Numerical study in 1D

We consider $\Omega=[0,1]$. The one dimensional differential eigenvalue
problem (2.1) has the following explicit eigenvalues and eigenfunctions

$\lambda_{j}=j^{2}\pi^{2},\quad\text{and}\quad u_{j}=\sqrt{2}\sin(j\pi
x),\quad j=1,2,\cdots,$

respectively.

Figure 2: Eigenvalue and eigenfunction errors when using $C^{2}$ cubic (left)
and $C^{3}$ quartic (right) elements with $N=100$ in 1D.

Figure 3: Eigenvalue and eigenfunction errors when using $C^{4}$ quintic
(left) and $C^{5}$ sextic (right) elements with $N=100$ in 1D.

Figures 2 and 3 show the relative eigenvalue and eigenfunction errors in
energy norm errors of $C^{2}$ cubic, $C^{3}$ quartic, $C^{4}$ quintic, and
$C^{5}$ sextic isogeometric approximations using the boundary penalization,
respectively. In all the cases, we use $N=100$ uniform elements. We set the
penalty parameters to the default value of 1 in all cases. The solid lines
show the original error curves of the isogeometric elements, while the dotted
lines show the discretely corrected errors. In all the cases, the consistent
boundary penalization removes outliers from the spectrum. Table 1 shows the
optimal error convergence rates of approximate eigenvalues and eigenfunctions
by the discrete correction.

For cubic and quartic elements, our experiments show that the boundary
penalization with slightly larger penalty parameters (in the scale of
$\mathcal{O}(h^{-1})$) also removes the outliers while maintaining accuracy in
the lower-frequency region. If the penalty parameters are large (in the scale
of $\mathcal{O}(h^{-2})$), then we lose accuracy in the lower-frequency
region. We observe similar behavior for quintic and sextic elements.

$p$ | $N$ | $\frac{|\lambda_{1}^{h}-\lambda_{1}|}{\lambda_{1}}$ | $|u_{1}^{h}-u_{1}|_{H^{1}}$ | $\|u_{1}^{h}-u_{1}\|_{L^{2}}$ | $\frac{|\lambda_{6}^{h}-\lambda_{6}|}{\lambda_{6}}$ | $|u_{6}^{h}-u_{6}|_{H^{1}}$ | $\|u_{6}^{h}-u_{6}\|_{L^{2}}$
---|---|---|---|---|---|---|---
| 8 | 1.31E-07 | 1.14E-03 | 2.31E-05 | 2.99E-02 | 4.06E+00 | 1.29E-01
| 16 | 1.93E-09 | 1.38E-04 | 1.38E-06 | 1.60E-04 | 2.45E-01 | 2.91E-03
3 | 32 | 2.98E-11 | 1.71E-05 | 8.48E-08 | 1.63E-06 | 2.42E-02 | 1.27E-04
| 64 | 5.61E-13 | 2.14E-06 | 5.28E-09 | 2.25E-08 | 2.83E-03 | 7.11E-06
| $\rho_{3}$ | 5.95 | 3.02 | 4.03 | 6.76 | 3.48 | 4.7
| 8 | 1.76E-07 | 1.32E-03 | 6.09E-05 | 1.49E-01 | 1.09E+01 | 4.37E-01
4 | 16 | 3.22E-10 | 5.47E-05 | 1.24E-06 | 4.49E-04 | 4.31E-01 | 1.00E-02
| 32 | 5.26E-13 | 2.06E-06 | 2.31E-08 | 8.70E-07 | 1.61E-02 | 1.81E-04
| $\rho_{4}$ | 9.17 | 4.66 | 5.68 | 8.69 | 4.70 | 5.62
Table 1: Relative eigenvalue and eigenfunction errors for $C^{2}$-cubic and
$C^{3}$-quartic elements in 1D.

####  6.1.2 Numerical study in 2D

Let $\Omega=[0,1]\times[0,1]$. The two dimensional differential eigenvalue
problem (2.1) has exact eigenvalues and eigenfunctions

$\lambda_{jk}=(j^{2}+k^{2})\pi^{2},\quad\text{and}\quad u_{jk}=2\sin(j\pi
x)\sin(k\pi y),\quad j,k=1,2,\cdots,$

respectively. Figures 4 and 5 show the relative eigenvalue errors of $C^{2}$
cubic, $C^{3}$ quartic, $C^{4}$ quintic, and $C^{5}$ sextic isogeometric
elements, respectively. We build the basis functions on a uniform mesh with
$40\times 40$ elements. We set the penalty parameters to 1 as a default. This
penalization removes all outliers. More importantly, the corrected relative
eigenvalue error curves remain almost flat when compared to the original ones.
That is, the overall spectral approximation improves significantly. For
brevity, we do not include a full set of 3D results, but the performance is
similar; see the left plot in Figure 6 for an example of $C^{3}$ quartic
isogeometric analysis in 3D. Lastly, the right plot in Figure 6 shows the
comparison of the eigenvalue errors when using a non-uniform mesh with
$20\times 20\times 20$ elements in 3D. The mesh is created using the randomly
generated 1D grid with nodes at $x=$ 0, 0.0475, 0.0492, 0.1555, 0.1961,
0.2461, 0.2939, 0.3274, 0.4039, 0.4379, 0.4834, 0.5607, 0.6022, 0.6276,
0.7009, 0.7252, 0.8128, 0.8729, 0.8997, 0.9632, 1. We observe that the
eigenvalue errors in the high-frequency regions are reduced.

Figure 4: Eigenvalue errors when using $C^{2}$ cubic (left) and $C^{3}$
quartic (right) isogeometric elements with $N=40\times 40$ in 2D.

Figure 5: Eigenvalue errors when using $C^{4}$ quintic (left) and $C^{5}$
sextic (right) isogeometric elements with $N=40\times 40$ in 2D. Figure 6:
Eigenvalue errors when using $C^{3}$ quartic isogeometric analysis with
$N=20\times 20\times 20$ uniform elements (left) and non-uniform elements in
3D. The non-uniform grid is generated randomly in 1D and then the 3D mesh is
constructed based on the tensor-product.

###  6.2 Neumann eigenvalue problem

For simplicity, we report the results for 2D Neumann eigenvalue problem (4.1).
Figures 7 and 8 show the relative eigenvalue errors of $C^{1}$ quadratic,
$C^{2}$ cubic, $C^{3}$ quartic, and $C^{4}$ quintic isogeometric elements,
respectively. We use a uniform mesh with $40\times 40$ elements. We set the
penalty parameters to 1. The boundary penalization removes all the outliers,
and the eigenvalue error curves are almost flat compared to the original ones.

Figure 7: Eigenvalue errors of the Neumann eigenvalue problem when using
$C^{1}$ quadratic (left) and $C^{2}$ cubic (right) isogeometric elements with
$N=40\times 40$ in 2D.

Figure 8: Eigenvalue errors of the Neumann eigenvalue problem when using
$C^{3}$ quartic (left) and $C^{4}$ quintic (right) isogeometric elements with
$N=40\times 40$ in 2D.

###  6.3 Numerical study on condition numbers

We now study the condition numbers to show further advantages of the proposed
method. Due to the symmetry of the stiffness and mass matrices, the condition
numbers of the generalized matrix eigenvalue problems (2.8) and (3.4) are
given by

$\gamma:=\frac{\lambda^{h}_{\max}}{\lambda^{h}_{\min}},\qquad\tilde{\gamma}:=\frac{\tilde{\lambda}^{h}_{\max}}{\tilde{\lambda}^{h}_{\min}},$
(6.3)

where $\lambda^{h}_{\max},\tilde{\lambda}^{h}_{\max}$ are the largest
eigenvalues and $\lambda^{h}_{\min},\tilde{\lambda}^{h}_{\min}$ are the
smallest eigenvalues. The condition number characterizes the stiffness of the
system. We follow the recent work of soft-finite element method (softFEM) [30]
and define the condition number reduction ratio of the method with respect to
IGA as

$\rho:=\frac{\gamma}{\tilde{\gamma}}=\frac{\lambda^{h}_{\max}}{\tilde{\lambda}^{h}_{\max}}\cdot\frac{\tilde{\lambda}^{h}_{\min}}{\lambda^{h}_{\min}}.$
(6.4)

In general, one has $\lambda^{h}_{\min}\approx\tilde{\lambda}^{h}_{\min}$ for
IGA and the proposed method with a sufficient number of elements (in practice,
these methods with only a few elements already lead to good approximations to
the minimum eigenvalues). Thus, the condition number reduction ratio is mainly
characterized by the ratio of the largest eigenvalues. Finally, we define the
condition number reduction percentage as

$\varrho=100\frac{\gamma-\tilde{\gamma}}{\gamma}\,\%=100(1-\rho^{-1})\,\%.$
(6.5)

Table 2 shows the smallest and largest eigenvalues, condition numbers and
their reduction ratios and percentages for 1D, 2D, and 3D problems. We observe
that the condition numbers of the proposed method are significantly smaller.
The condition number of the proposed method reduces by about 32% for $C^{2}$
cubic, 60% for $C^{3}$ quartic, 75% for $C^{4}$ quintic, and 83% for $C^{5}$
sextic elements in 1D, 2D, and 3D. Similar stiffness reductions are observed
in [30] for softFEM and in [31] for IGA with optimally-blended quadrature
rules and discrete corrections.

$d$ | $p$ | $\lambda_{\min}^{h}$ | $\lambda_{\max}^{h}$ | $\tilde{\lambda}_{\max}^{h}$ | $\gamma$ | $\tilde{\gamma}$ | $\rho$ | $\varrho$
---|---|---|---|---|---|---|---|---
| 3 | 9.87 | 5.82E+05 | 3.95E+05 | 5.90E+04 | 4.00E+04 | 1.47 | 32.13%
1 | 4 | 9.87 | 9.80E+05 | 3.95E+05 | 9.93E+04 | 4.00E+04 | 2.48 | 59.69%
| 5 | 9.87 | 1.57E+06 | 4.16E+05 | 1.59E+05 | 4.22E+04 | 3.78 | 73.52%
| 6 | 9.87 | 2.38E+06 | 3.99E+05 | 2.41E+05 | 4.05E+04 | 5.96 | 83.22%
| 3 | 19.74 | 2.91E+05 | 1.98E+05 | 1.47E+04 | 1.00E+04 | 1.47 | 32.16%
2 | 4 | 19.74 | 4.90E+05 | 1.97E+05 | 2.48E+04 | 1.00E+04 | 2.48 | 59.69%
| 5 | 19.74 | 7.86E+05 | 2.01E+05 | 3.98E+04 | 1.02E+04 | 3.91 | 74.45%
| 6 | 19.74 | 1.19E+06 | 1.98E+05 | 6.03E+04 | 1.00E+04 | 6.01 | 83.36%
| 3 | 29.61 | 1.09E+05 | 7.41E+04 | 3.69E+03 | 2.50E+03 | 1.47 | 32.16%
3 | 4 | 29.61 | 1.84E+05 | 7.40E+04 | 6.20E+03 | 2.50E+03 | 2.48 | 59.69%
| 5 | 29.61 | 2.95E+05 | 7.44E+04 | 9.95E+03 | 2.51E+03 | 3.96 | 74.76%
| 6 | 29.61 | 4.46E+05 | 7.41E+04 | 1.51E+04 | 2.50E+03 | 6.02 | 83.40%
Table 2: Minimal and maximal eigenvalues, condition numbers, reduction ratios
and percentages when using IGA and DC-IGA. The polynomial degrees are
$p\in\\{3,4,5,6\\}$. There are 200, $100\times 100$, and $50\times 50\times
50$ elements in 1D, 2D, and 3D, respectively.

##  7 Concluding remarks

We remove the outliers from the discrete spectrum of isogeometric analysis.
The idea is to penalize the boundary elements with extra information from the
higher-order derivatives. By removing the outliers, we improve the overall
spectrum significantly, especially in multiple dimensions. Consequently, we
improve the approximate eigenfunctions in the high-frequency regions. The
methods add penalizations to the discrete systems only for the degrees of
freedom near the boundaries; thus, the boundary penalization implementation is
simple to perform in existing isogeometric analysis codes. Our numerical
simulation also shows that the technique works equally well for both the
Dirichlet and Neumann eigenvalue problems when using the open knots and non-
open knots. However, for curved boundary and domains with geometrical
mappings, a simple direct implementation does not improve much the high-
frequency eigenvalue errors, and further study on the boundary penalty terms
is necessary.

Future work will analyze the dispersion error for higher-order elements to
characterize the approximate eigenvalue errors analytically and remove the
outliers by imposing strong boundary conditions, which reconstructs the basis
functions near the boundary elements. We will also analyze the impact of the
outlier-removal techniques on time marching for time-dependent partial
differential equations, especially the effects of high-frequency error
pollutions in the time-marching schemes. From the study on stiffness reduction
in Section 5.3, we expect to have a larger critical time step allowing for
faster explicit in time isogeometric analysis of engineering problems.

## Acknowledgments

The authors thank professor Thomas J.R. Hughes (UT Austin) for several
discussions and comments which led to an improvement of the paper. This
publication was also made possible in part by the CSIRO Professorial Chair in
Computational Geoscience at Curtin University and the Deep Earth Imaging
Enterprise Future Science Platforms of the Commonwealth Scientific Industrial
Research Organisation, CSIRO, of Australia. This project has received funding
from the European Union’s Horizon 2020 research and innovation program under
the Marie Sklodowska-Curie grant agreement No 777778 (MATHROCKS). The Curtin
Corrosion Centre and the Curtin Institute for Computation kindly provide
ongoing support.

## References

  * [1] T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer methods in applied mechanics and engineering 194 (39) (2005) 4135–4195.
  * [2] J. A. Cottrell, T. J. R. Hughes, Y. Bazilevs, Isogeometric analysis: toward integration of CAD and FEA, John Wiley & Sons, 2009.
  * [3] J. A. Cottrell, A. Reali, Y. Bazilevs, T. J. R. Hughes, Isogeometric analysis of structural vibrations, Computer methods in applied mechanics and engineering 195 (41) (2006) 5257–5296.
  * [4] T. J. R. 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, Computer methods in applied mechanics and engineering 197 (49) (2008) 4104–4124.
  * [5] T. J. R. Hughes, J. A. Evans, A. Reali, Finite element and NURBS approximations of eigenvalue, boundary-value, and initial-value problems, Computer Methods in Applied Mechanics and Engineering 272 (2014) 290–320.
  * [6] M. Bartoň, V. Calo, Q. Deng, V. Puzyrev, Generalization of the Pythagorean eigenvalue error theorem and its application to isogeometric analysis, in: Numerical Methods for PDEs, Springer, 2018, pp. 147–170.
  * [7] V. Calo, Q. Deng, V. Puzyrev, Dispersion optimized quadratures for isogeometric analysis, Journal of Computational and Applied Mathematics 355 (2019) 283–300.
  * [8] V. Calo, Q. Deng, V. Puzyrev, Quadrature blending for isogeometric analysis, Procedia Computer Science 108 (2017) 798–807.
  * [9] Q. Deng, M. Bartoň, V. Puzyrev, V. Calo, Dispersion-minimizing quadrature rules for C1 quadratic isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 328 (2018) 554–564.
  * [10] V. Puzyrev, Q. Deng, V. M. Calo, Dispersion-optimized quadrature rules for isogeometric analysis: modified inner products, their dispersion properties, and optimally blended schemes, Computer Methods in Applied Mechanics and Engineering 320 (2017) 421–443.
  * [11] V. Puzyrev, Q. Deng, V. Calo, Spectral approximation properties of isogeometric analysis with variable continuity, Computer Methods in Applied Mechanics and Engineering 334 (2018) 22–39.
  * [12] Q. Deng, V. Calo, Dispersion-minimized mass for isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 341 (2018) 71–92.
  * [13] G. Strang, G. J. Fix, An analysis of the finite element method, Vol. 212, Prentice-hall Englewood Cliffs, NJ, 1973.
  * [14] K. J. Marfurt, Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations, Geophysics 49 (5) (1984) 533–549.
  * [15] M. Ainsworth, H. A. Wajid, Optimally blended spectral-finite element scheme for wave propagation and nonstandard reduced integration, SIAM Journal on Numerical Analysis 48 (1) (2010) 346–371.
  * [16] Q. Deng, V. Puzyrev, V. Calo, Optimal spectral approximation of 2n-order differential operators by mixed isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 343 (2019) 297–313.
  * [17] Q. Deng, V. Puzyrev, V. Calo, Isogeometric spectral approximation for elliptic differential operators, Journal of Computational Science (2018).
  * [18] Q. Deng, Analytical solutions to some generalized and polynomial eigenvalue problems, Special Matrices 9 (1) (2021) 240–256.  doi:doi:10.1515/spma-2020-0135.
URL https://doi.org/10.1515/spma-2020-0135

  * [19] E. Sande, C. Manni, H. Speleers, Sharp error estimates for spline approximation: Explicit constants, n-widths, and eigenfunction convergence, Mathematical Models and Methods in Applied Sciences 29 (06) (2019) 1175–1205.
  * [20] H. Brezis, Functional analysis, Sobolev spaces and partial differential equations, Universitext, Springer, New York, 2011.
  * [21] C. De Boor, A practical guide to splines, Vol. 27, Springer-Verlag New York, 1978\.
  * [22] L. Piegl, W. Tiller, The NURBS book, Springer Science & Business Media, 1997\.
  * [23] A. Buffa, C. De Falco, G. Sangalli, Isogeometric analysis: new stable elements for the Stokes equation, International Journal for Numerical Methods in Fluids (2010).
  * [24] J. A. Evans, T. J. Hughes, Isogeometric divergence-conforming b-splines for the Darcy–Stokes–Brinkman equations, Mathematical Models and Methods in Applied Sciences 23 (04) (2013) 671–741.
  * [25] H. M. Hilber, T. J. Hughes, R. L. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics, Earthquake Engineering & Structural Dynamics 5 (3) (1977) 283–292.
  * [26] J. Chung, G. Hulbert, A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-$\alpha$ method, Journal of Applied Mechanics 60 (2) (1993) 371–375.
  * [27] P. N. Vabishchevich, Additive Operator-Difference Schemes: Splitting Schemes, Walter de Gruyter, 2013.
  * [28] P. Behnoudfar, V. M. Calo, Q. Deng, P. D. Minev, A variationally separable splitting for the generalized-$\alpha$ method for parabolic equations, International Journal for Numerical Methods in Engineering 121 (5) (2020) 828–841.
  * [29] L. Gao, Kronecker products on preconditioning, Ph.D. thesis, King Abdullah University of Science and Technology (2013).
  * [30] Q. Deng, A. Ern, SoftFEM: revisiting the spectral finite element approximation of elliptic operators, arXiv preprint arXiv:2011.06953 (2020).
  * [31] Q. Deng, V. Calo, Outlier removal for isogeometric spectral approximation with the optimally-blended quadratures, Procedia Computer Science, accepted (2021).


Learning differential equations from data

K. D. Olumoyin

ABSTRACT. Differential equations are used to model problems that originate in
disciplines such as physics, biology, chemistry, and engineering. In recent
times, due to the abundance of data, there is an active search for data-driven
methods to learn Differential equation models from data. However, many
numerical methods often fall short. Advancements in neural networks and deep
learning, have motivated a shift towards data-driven deep learning methods of
learning differential equations from data. In this work, we propose a forward-
Euler based neural network model and test its performance by learning ODEs
such as the FitzHugh-Nagumo equations from data using different number of
hidden layers and different neural network width.

Keywords and phrases: Neural Network, Differential equations, Modeling

1\. INTRODUCTION

In the late 1500s, Danish astronomer, Tycho Brahe used geometry and the
astronomical tools of his time to collect accurate data about planetary
positions. Later, his student Johannes Kepler used Brahe’s records to
formulate the laws of planetary motion [3]. This was one of the earliest
recorded use of a data-driven approach to formulate and discover the governing
laws of our physical and celestial world. On the other hand, we have the
fascinating discovery of the heat equation by Mathematician Joseph Fourier
[2], his was a triumph of intuition and reasoning. And for many decades, such
reasoning steered our knowledge of modeling with differential equations.

In recent years, we are seeing a resurgence of data-driven mathematical
modeling. We are no longer turning to first principles but increasingly to
data-driven approaches to understanding the patterns in complex systems. These
systems have the potential to help us understand finance, epidemiology,
biology, physics, chemistry, and many other fields. These systems are usually
nonlinear and high-dimensional and with applications of mathematical methods
in combination with the enormous data available and the improved computational
resources, we can make more accurate predictions, forecasts, and
extrapolations.

We are interested in the dynamical systems that describes observed dynamics of
biochemical nature such as the Fitzhugh-Nagumo models [6, 5]. In recent years,
there is a shift in the modeling approach from first principle to data-driven
approaches. This is a merge of data usually big data and deep learning in the
analysis and understanding of dynamical systems. This gives us an opportunity
to validate assumptions about underlying physical and biological laws that
govern these observed dynamics.

2\. PRELIMINARY

Advancement in machine learning and data science have made it possible to
learn complex patterns from data. Data-driven approaches including neural
networks and nonlinear regression have been applied to dynamical systems [7,
9, 8]. Data-driven approaches are not replacing good intuition and reasoning
but instead are unlocking new possibilities [1]. In [12], the authors used
Gaussian processes to discover differential equations from data. The
limitation with their approach is that they can only learn linear systems. An
advancement was made in [7], where the authors formulated a sparse
identification of a nonlinear dynamical system for extracting governing
equations from data. Their approach is able to learn nonlinear systems however
it can easy become computationally expensive when scaled to higher dimensions.
It requires some knowledge of the dynamics while constructing the library of
functions that aid the neural network to learn the dynamics. One of the most
successful approaches in the last five years is the physics-informed neural
network (PINN) introduced in [10], where the form of the differential equation
is assumed to be known and the task is to learn the parameters of the
differential equation from data.PINN excels when trained with sparse data
because it uses the governing equations as part of the training process. PINN
is generally limited to low-dimensional parameter identification. It has not
performed well in modeling differential equations with sharp gradients, and it
has become necessary to use collocation points in the differential equation
residuals. This results in a huge computational cost in training. However, in
a true sense of learning differential equations from data, we desire
approaches that do not assume that the dynamics of the differential equation
is known. The following works have succeeded in modeling differential
equations from data, where the dynamics of the differential equation is not
assumed to be known [9] and [8]. In this work, we propose a forward-Euler
feedforward neural network.

3\. METHODOLOGY

In this work, we will do the following:

  1. (1)

Formulate a neural network, This network is used to advance the neural network
prediction of the differential equation from a temporal state to the next
higher temporal state.

  2. (2)

Imposition of initial and boundary conditions so that we can have a well-posed
optimization problem during the training.

  3. (3)

A residual differential equation connection inspired by the forward-Euler
scheme will be used to construct the loss function. If we seek to improve the
network, we will have to use higher-order temporal schemes to design the loss
function. The use of higher-order temporal schemes will require changes in the
network architecture. In this work, we want to demonstrate that a loss
function built using a forward Euler temporal scheme in the network
formulation is capable to learn the form of a differential equation from data
to high accuracy.

Our goal is to develop a neural network that can be used in an inverse
problem, where data is available. Neural network approaches are suitable for
differential equations because they can scale to high dimensions. They can
also solve forward and inverse problems with sparse data. As such, the
formulation of such a network can be used as a new numerical solver for
various time-dependent differential equations given specific initial and
boundary conditions.

In totality, our approach provides an approximation of a discrete solution
$\mathbf{u}(x,t;\theta)$ denoted as $\mathbf{u}^{pred}$ that satisfies (3) for
specifically given an initial condition. Here $\theta$ denotes the neural
network trainable parameters.

3.1 Data generation

We will learn the FitzHugh-Nagumo system (FHN) [6, 5] from data. FHN is a
simplification of the Hodgkin-Huxley model [4]. The motivation for the FHN
model was to isolate conceptually the essential mathematical properties of
excitation and propagation from the electrochemical properties of sodium and
potassium ion flow. The model consists of a voltage-like variable $\it{u}$
having cubic nonlinearity that allows regenerative self-excitation via a
positive feedback, and a recovery variable $\it{v}$ having a linear dynamics
that provides a slower negative feedback.

The partial differential equation version of FHN is given below:

$\begin{split}\frac{\partial u}{\partial
t}&=\frac{1}{c}(u-\frac{1}{3}u^{3}-v)+\triangle u\\\ \frac{\partial
v}{\partial t}&=c(u-av+b)\end{split}$ (1)

The parameters in (1) satisfies the following inequalities, $0<a<1$, $b>0$,
$c>0$. An appropriate initial condition is chosen by finding an homogeneous
solution of (1), which satisfies the ordinary differential equations below
(2):

$\begin{split}\frac{du}{dt}&=\frac{1}{c}(u-\frac{1}{3}u^{3}-v)\\\
\frac{dv}{dt}&=c(u-av+b)\end{split}$ (2)

The equilibrium solution of (2) is $(u_{e},v_{e})$ as shown in Figure (1).
Where the nullclines of the FHN equations are a line and a cubic that
intersect in a single rest point.

Figure 1. Diagram of the nullcline and critical manifold for homgeneous
FitzHugh-Nagumo model.

We consider the following initial conditions $u(x,0)$ and $v(x,0)$ for system
(1) given as follows:

$u(x,0)=u_{0}(x),v(x,0)=v_{0}(x)$

where,

$u_{0}(x)=\begin{cases}u_{e}&\text{if $x\geq\frac{L}{2}$}\\\ -u_{e}&\text{if
$x<\frac{L}{2}$}\end{cases}$

$v_{0}(x)=-v_{e}$

and $x\in[0,L]$.

Training and testing dataset is generated by solving (2) using ode solvers
such as scipy’s odeint function.

3.2 A discrete physics loss function based Neural Network

We consider the general form of a multi-dimensional, nonlinear, coupled PDE
system of the form

$\mathbf{u}_{t}(x,t)=\mathcal{F}[\mathbf{u},\mathbf{u}^{2},\ldots,\mathbf{u}_{x},\mathbf{u}_{xx},\ldots;\bf{\lambda}]$
(3)

where $\mathbf{u}(x,t)$ is the physical quantity we want to study. The
temporal domain is $t\in[0,T]$ and the spatial domain is
$x\in\Omega\subset\mathbf{R}$. $\mathbf{u}_{t}(x,t)$ is the first-order time
derivative term and $\mathbf{u}_{x}$ denotes the gradient operator with
respect to the spatial variable $x$. $\mathcal{F}[\cdot]$ is the nonlinear
functional parameterized by $\bf{\lambda}$. Eq (3) satisfies appropriate
initial and boundary conditions.

We will learn the FitzHugh-Nagumo system(FHN) [6, 5] from data. We consider
(2) in 1 spatial dimensions, and $t>0$.Here, we assume only a temporal domain
$t\in[0,T]$, and $\mathbf{u}$ denotes the variables $(u,v)$. In this approach,
we do not assume the form of the right hand side of (2). Suppose
$\mathbf{\hat{u}}$ is the target variable for which there is a known
measurements, we can formulate a neural network that take as input;
$\mathbf{\hat{u}}$, and the temporal domain $t\in[0,T]$. The output of this
neural network is the learned right hand side of (2) denoted as $\mathcal{F}$.
This neural network is a feedforward neural network, see Figure (2).

Figure 2. Schematic of ForwardEulerNet

The loss function can then be formulated as follows

$\mathcal{R}(t;\theta)=\mathbf{\hat{u}}_{t}-\mathcal{F}[t,\mathbf{\hat{u}},\ldots;\bf{\lambda},\bf{\theta}]$
(4)

A discrete loss is used in training the neural network, where the discrete
$\mathcal{F}$ is given by $\mathcal{N}$. The discrete loss function is
inspired by the forward Euler method to connect the temporal state
$\mathbf{\hat{u}}_{i}$ and temporal state $\mathbf{\hat{u}}_{i+1}$.

$\mathbf{\hat{u}}_{i+1}=\mathbf{\hat{u}}_{i}+\partial
t\cdot\mathcal{N}[\mathbf{\hat{u}};\theta],\quad\partial
t=\frac{T}{K},K\in\mathbf{R}^{+}$ (5)

The trainable parameters $\theta$ are optimized during training of the
forward-Euler based neural network that minimizes the mean squared error loss
function (MSE) given as follows

$MSE=\sum_{i=1}^{K}|-\mathbf{\hat{u}}_{i+1}+\mathbf{\hat{u}}_{i}+\partial
t\cdot\mathcal{N}[\mathbf{\hat{u}};\theta]|^{2},\quad\partial
t=\frac{T}{K},K\in\mathbf{R}^{+}$ (6)

When training is done, we obtain $\mathcal{N}$, which is the learned dynamics
of an ODE system such as eq (2). Now to obtain learned $\mathbf{u}$ denoted as
$\mathbf{u}^{pred}$, we use an ODE solver such. Using a solver such as the
scipy odeint. We observe that the learned model parameters $\lambda=[a,b,c]$
in eq (2) are embedded into the dynamics $\mathcal{N}$.

$\mathbf{u}^{pred}=odeint(\mathcal{N},(u(0),v(0)),tspan)$

In Figures (3)–(8), we obtained $\mathbf{u}^{pred}$ and compared with the
target $\mathbf{\hat{u}}$. In Figure (8), we observe that when we add $5\%$
noise to the training data, we are no longer able to learn an accurate
dynamics of (2) from data.

Figure 3. No noise, $\triangle t$ = 0.01, 64 neurons  Figure 4. No noise,
$\triangle t$ = 0.025, 64 neurons  Figure 5. No noise, $\triangle t$ = 0.05,
64 neurons  Figure 6. $2\%$ noise, $\triangle t$ = 0.01, 64 neurons  Figure 7.
$2\%$ noise, $\triangle t$ = 0.01, 128 neurons  Figure 8. $5\%$ noise,
$\triangle t$ = 0.01, 64 neurons

4\. CONCLUDING REMARKS

We present results and performance of the forward-Euler neural network
approach that learns the FitzHugh-Nagumo ODE (2) from data. In Table 1, we
observe that increasing the number of hidden layers, the MSE error of the
variable $u$ produces slightly more accurate results when the width (number of
neurons per layer) of the neural network is 32, there are not much difference
in the MSE error when we increase the hidden layers and the width is 64, 128,
and 256. However, the benefit of increasing the width is more noticeable in
the case of one hidden layer.

Width
---
Hidden layers |  32 |  64 |  128 |  256
1 |  7.6e-03 |  4.3e-03 |  1.3e-03 |  2.3e-03
2 |  4.8e-03 |  3.0e-03 |  2.3e-03 |  3.53-03
3 |  4.9e-03 |  4.4e-03 |  2.8e-03 |  3.3e-03
Table 1. Table showing MSE error of the variable $u$, where we assumed a noise
free training data and set the time step to $0.01$.

REFERENCES

## References

  * [1] K. Olumoyin, A.Q.M. Khaliq, and K.M. Furati, Data-Driven Deep-Learning Algorithm for Asymptomatic COVID-19 Model with Varying Mitigation Measures and Transmission Rate, Epidemiologia 2 471-489, 2021.
  * [2] J. Fourier, Theorie Analytique de la Chaleur Firmin Didot (reissued by Cambridge University Press),1822.
  * [3] S. Hawkings, On the Shoulders of Giants: The Great Works of Physics and Astronomy, Running Press, 2002.
  * [4] A.L. Hodgekin and A.F. Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve, J Physiol. 117 500 - 544, 1952.
  * [5] J. Nagumo, An active pulse transmission line simulating nerve axon Proc. IRE. 50 2061 - 2070, 1962.
  * [6] R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane Biophysical. 1 445 - 466, 1961.
  * [7] S. Brunton, J. Proctor, and J. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems Proceedings Of The National Academy Of Sciences 113 3932 - 3937, 2016.
  * [8] P. Ren, C. Rao, Y. Liu, J. Wang, and H. Sun, PhyCRNet: Physics-informed Convolutional-Recurrent Network for solving Spatiotemporal PDEs ArXiv:2106.14103v1 [cs.LG], 2021.
  * [9] N. Geneva and N. Zabaras, Modeling the dynamics of PDE systems with physics-constrained deep auto-regressive networks Journal Of Computational Physics 403 109056, 2020.
  * [10] M. Raissi, P. Perdikaris, and G. Karniadakis, Physics informed deep learning: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations Journal Of Computational Physics 686-707, 2019.
  * [11] G. Pang, L. Lu, and G. Karniadakis, fPINNs: Fractional Physics-informed Neural Networks SIAM J. Sci. Comput. 41, A2603-A2626, 2019.
  * [12] G. Mamikon, M. Raissi, P. Perdikaris, and G. Karniadakis, Machine learning of space fractional differential equations SIAM J. Sci. Comput. 41, A2485-A2509, 2019.

DEPARTMENT OF MATHEMATICS, Middle Tennessee State University, Murfreesboro,
TN, USA

E-mail address<EMAIL_ADDRESS>