# Condition number estimates for matrices arising in NURBS based isogeometric
discretizations of elliptic partial differential equations††thanks: This
research was supported in part by the Austrian Sciences Fund (Project
P21516-N18), National Science Foundation grant EPS-1135483, Award No.
KUS-C1-016-04, made by King Abdullah University of Science & Technology
(KAUST), and the KAUST Numerical Porous Media Center.

Krishan P. S. Gahalaut Division of Computer, Electrical and Mathematical
Sciences and Engineering, King Abdullah University of Science and Technology,
Thuwal 23955-6900, Kingdom of Saudi Arabia. (). Questions, comments, or
corrections to this document may be directed to that email address.
<EMAIL_ADDRESS>Satyendra K. Tomar Dornacher Strasse 6/21, 4040
Linz, Austria<EMAIL_ADDRESS>Craig C. Douglas School of Energy
Resources and Mathematics Department, University of Wyoming, Laramie, WY
82071-3036, USA<EMAIL_ADDRESS>
###### Abstract

We derive bounds for the minimum and maximum eigenvalues and the spectral
condition number of matrices for isogeometric discretizations of elliptic
partial differential equations in an open, bounded, simply connected Lipschitz
domain $\Omega\subset\mathbb{R}^{d}$, $d\in\\{2,3\\}$. We consider refinements
based on mesh size $h$ and polynomial degree $p$ with maximum regularity of
spline basis functions. For the $h$-refinement, the condition number of the
stiffness matrix is bounded above by a constant times $h^{-2}$ and the
condition number of the mass matrix is uniformly bounded. For the
$p$-refinement, the condition number grows exponentially and is bounded above
by $p^{2d+2}4^{pd}$ and $p^{2d}4^{pd}$ for the stiffness and mass matrices,
respectively. Rigorous theoretical proofs of these estimates and supporting
numerical results are provided.

###### keywords:

Elliptic PDEs, Galerkin formulation, B-Splines, NURBS, Isogeometric method,
Stiffness matrix, Mass matrix, $h$-refinement, $p$-refinement, Eigenvalues,
Condition number

sinumxxxxxxxx–x

##  1 Introduction

Isogeometric analysis is a term introduced by Hughes et al. in 2005 [28]. Most
of the research activity in isogeometric analysis has focused on using Non-
Uniform Rational B-Spline (NURBS) as basis functions, e.g., [1, 5, 9, 28].
Isogeometric analysis is not restricted to NURBS basis functions. Other types
of basis functions are used by researchers, e.g., T-Splines, hierarchical
B-Splines, and subdivision schemes. Use of splines as finite element basis
functions dates back to the 1970’s, however [41, 46, 50].

In isogeometric analysis the computational geometry (e.g., a circle) is
represented exactly from the information and the basis functions given by
Computer Aided Design (CAD). It holds an advantage over classical finite
element methods (FEM), where the basis functions are defined using piecewise
polynomials and the computational geometry (i.e., a mesh) is defined on
polygonal elements. It has been argued in [9] that NURBS based isogeometric
method leads to qualitatively more accurate results than a standard piecewise
polynomial based finite element method. Typically, the solution computed by an
isogeometric method has a higher continuity than the one computed in a
classical finite element method. It is a difficult and cumbersome task to
achieve even $C^{1}$ inter-element continuity in the piecewise polynomial
based finite element method, whereas isogeometric method offers up to
$C^{p-m}$ continuity, where $p$ denotes the degree of the basis functions and
$m$ denotes the knot-multiplicity. Finally, isogeometric analysis provides a
powerful tool to compute highly continuous numerical solution of PDEs arising
in engineering sciences.

Since the introduction of isogeometric analysis, most of its progress has been
focused on the applications and discretization properties. Nevertheless, when
dealing with large problems, the cost of solving the linear system of
equations arising from the isogeometric discretization becomes an important
issue. Clearly, the discretization matrix $A$ gets denser by increasing the
polynomial degree $p$. Therefore, the cost of a direct solver, particularly
for large problems, becomes prohibitively expensive. The most practical way to
solve them is to resort to an iterative method. Since the convergence rate of
such methods is strongly affected by the condition number of the system matrix
$A$, it is important to assess this quantity as a function of the mesh size
$h$ for the $h$-refinement, or as a function of the degree $p$ for the
$p$-refinement. Note that in the $p$-refinement, improved approximate
solutions are sought by increasing $p$ while the mesh of the domain, and thus
the maximum quadrilateral diameter $h$, is held fixed, whereas in the
$h$-refinement, improved approximations are obtained by refining the mesh, and
thus reducing $h$, while $p$ is held fixed. In this paper we consider both the
cases: $h$\- and $p$-refinements. Similar efforts are made in [22] on the
spectrum of stiffness matrices and in [40] on bounding the influence of the
domain parameterization and knot spacing. However, these papers primarily
derive bounds with respect to the mesh size $h$. To the best of our knowledge
there is no study that discusses the bounds on condition number estimates of
isogeometric matrices with respect to $p$-refinement. Our main results provide
upper bounds for the condition number of the stiffness matrix and the mass
matrix for both the $h$\- and $p$-refinements.

For $h$-refinement applied to second order elliptic problems on a regular
mesh, the condition number of the finite element stiffness matrix scales as
$h^{-2}$ and the condition number of the mass matrix is bounded uniformly,
independent of $h$ [3]. This is true for a great variety of elements and
independent of the dimension of the problem domain. Our results here are in
agreement [2] and are useful in theoretical analysis that relates to
$h$-refinement. For example, in convergence analysis of multigrid methods,
these results are one of the key elements in deriving convergence factors, for
finite element analysis [7, 23, 43] and for isogeometric analysis [20].

The order of the approximation error of the numerical solution depends on the
choice of the finite dimensional subspace, not on the choice of its basis [8].
Therefore, when working with a finite element method or an isogeometric method
for elliptic problems, we only consider function spaces rather than the choice
of particular basis functions. Nevertheless, the choice of the basis functions
affects the condition number of the stiffness and the mass matrices, which
influences the performance of iterative solvers. There is no general theory to
characterize the extremal eigenvalues or the condition number based on a set
of general polynomial basis functions [4, 33, 34, 35]. Unlike the
$h$-refinement case, the condition number heavily depends on the choice of
basis functions for the $p$-refinement.

For different choices of basis functions the condition number may grow
algebraically or exponentially. Olsen and Douglas [37] estimated the condition
number bounds of finite element matrices for tensor product elements with two
choices of basis functions. For Lagrange elements, it is proved that the
condition number grows exponentially in $p$. For hierarchical basis functions
based on Chebychev polynomials, the condition number grows rapidly but only
algebraically in $p$. Similar results on the condition number bounds can be
found in [18, 27, 32].

Due to the larger support of NURBS basis functions, the band of the stiffness
matrix corresponding to the NURBS-based isogeometric method is less sparse
than the one arising from piecewise polynomial finite element procedures.
Therefore, a larger condition number is expected. Our results for the
$p$-refinement case show that the condition number of system matrices in an
isogeometric method grows exponentially.

Throughout this paper we deal with the maximum regularity $C^{p-1}$ of a
B-spline unless otherwise specified. The generic constant $C$, which will be
used often takes different values at different occasions, and is independent
of $h$ and $p$ in the analysis with respect to $h$-refinement and
$p$-refinement, respectively. Moreover, in our numerical studies the coarsest
and finest meshes use $h=1$ and $h=1/128$, respectively.

The remainder of the paper is organized as follows. In Section 2, we describe
the model problem and its discretization. In Section 3, we define B-Splines
and NURBS and a basic notation. We recall bounds for the condition number of a
B-Spline basis function. In Section 4, we derive bounds for the eigenvalues
and the condition number of the stiffness and mass matrices arising in
isogeometric discretizations for the $h$\- and $p$-refinement cases. In
Section 5, we provide numerical experiments that support the theoretical
estimates. In Section 6, we draw some conclusions and discuss future work.

##  2 Model problem and its discretization

Let $\Omega\subset\mathbb{R}^{d}$, $d\in\\{2,3\\}$, be an open, bounded, and
simply connected Lipschitz domain with Dirichlet boundary $\partial\Omega$. We
consider the Poisson equation,

(1a) $\displaystyle\Delta u$ $\displaystyle=-f\quad\mathrm{in~{}}\quad\Omega,$
(1b) $\displaystyle u$
$\displaystyle=0\quad\mathrm{on~{}}\quad\partial\Omega,$

where $f:\Omega\rightarrow\mathbb{R}$ is given. The aim is to find
$u:(\Omega\cup\partial\Omega)\rightarrow\mathbb{R}$ that satisfies (1). We
consider Galerkin’s formulation of the problem, which is commonly used in
isogeometric analysis. Since we are interested in the study of the condition
number, therefore we shall not go into the details of the solution properties,
and restrict ourselves to the study of the condition number of the resulting
system matrices.

Isogeometric analysis has the same theoretical foundation as finite element
analysis, namely the variational form of a partial differential equation. We
define the function space $\mathcal{S}$ as all the functions that have square
integrable derivatives and also satisfy $u|_{\partial\Omega}=0$,

(2) $\mathcal{S}=\\{u:u\in H^{1}(\Omega),u|_{\partial\Omega}=0\\},$

where $H^{1}(\Omega)=\\{u:D^{\alpha}u\in L^{2}(\Omega),|\alpha|\leq 1\\}$ is a
Sobolev space, $\alpha\in\mathbb{N}^{d}$ is a multi-index,
$D^{\alpha}=D_{1}^{\alpha_{1}}D_{2}^{\alpha_{2}}\ldots D_{d}^{\alpha_{d}}$,
and $D_{i}^{j}=\displaystyle\frac{\partial^{j}}{\partial x_{i}^{j}}$.

We write the variational formulation of the model problem by multiplying it by
an arbitrary function $v\in\mathcal{S}$ and integrating by parts. For a given
$f$: find $u\in\mathcal{S}$ such that for all $v\in\mathcal{S}$,

$\int_{\Omega}\nabla u\cdot\nabla v\hskip
4.26773ptd\Omega=\int_{\Omega}fv\hskip 4.26773ptd\Omega.$

We rewrite the formulation as: find $u\in\mathcal{S}$ such that for all
$v\in\mathcal{S}$,

(3) $a(u,v)=L(v),\quad\forall v\in\mathcal{S},$

where

$a(u,v)=\int_{\Omega}\nabla u\cdot\nabla v\hskip 4.26773ptd\Omega,\quad\text{
and }\quad L(v)=\int_{\Omega}fv\hskip 4.26773ptd\Omega.$

Note that $a(\cdot,\cdot)$ is a bilinear form that is continuous and coercive
on $\mathcal{S}$. $L(\cdot)$ is a linear form associated with the original
equation.

Let $\mathcal{S}^{h}\subset\mathcal{S}$ be a finite dimensional approximation
of $\mathcal{S}$. The Galerkin form of the problem is: Find
$u^{h}\in\mathcal{S}^{h}$ such that for all $\text{ }v^{h}\in\mathcal{S}^{h}$,

(4) $a(u^{h},v^{h})=L(v^{h}),$

which is a well-posed problem with a unique solution [8].

By approximating $u_{h}$ and $v_{h}$ using spline (see Section 3) basis
functions $N_{i}$, $i=1,2,\ldots,n_{h}=\mathcal{O}(h^{-2})$, the variational
formulation (4) is transformed into a set of linear algebraic equations,

(5) $A\bf u=f.$

$A$ denotes the stiffness matrix obtained from the bilinear form
$a(\cdot,\cdot)$,

$A=(a_{i,j})=(a(N_{i},N_{j})),\quad i,j=1,2,3,\ldots.,n_{h}.$

$\bf u$ denotes the vector of unknown degrees of freedom and $\bf f$ denotes
the right hand side vector from the known data of the problem. $A$ is a real,
symmetric positive definite matrix.

##  3 Splines and their condition number bounds

Non-uniform rational B-Splines (NURBS) are commonly used in isogeometric
analysis and are built from B-Splines. In Section 3.1, we give a brief
description of B-Splines and NURBS and their properties. In Section 3.2, we
define the derivatives of B-Splines. In Section 3.3, we prove bounds on the
condition number of B-Spline basis functions.

###  3.1 B-Splines and NURBS

In this section, we define B-Spline and NURBS functions. We also define
surfaces and describe higher order objects based on both types of functions.

The Cox-de Boor reursion formula [13] is given by

######  Definition 1.

Let $\Xi_{1}=\\{\xi_{i}:i=1,\ldots,n+p+1\\}$ be a non-decreasing sequence of
real numbers called the $knot$ $vector$, where $\xi_{i}$ is the $i^{th}$ knot,
$p$ is the polynomial degree, and $n$ is the number of basis function. With a
knot vector in hand, the B-Spline basis functions denoted by $N^{p}_{i}(\xi)$
are (recursively) defined starting with a piecewise constant $(p=0)$:

(6a) $\displaystyle N^{0}_{i}(\xi)$ $\displaystyle=\begin{cases}1&\text{if
$\xi\in[\xi_{i},\xi_{i+1})$,}\\\ 0&\text{otherwise,}\end{cases}$ (6b)
$\displaystyle\quad N^{p}_{i}(\xi)$
$\displaystyle=\frac{\xi-\xi_{i}}{\xi_{i+p}-\xi_{i}}N^{p-1}_{i}(\xi)+\frac{\xi_{i+p+1}-\xi}{\xi_{i+p+1}-\xi_{i+1}}N^{p-1}_{i+1}(\xi),$

where $0\leq i\leq n,p\geq 1$ and $\displaystyle\frac{0}{0}$ is considered as
zero.

For a B-Spline basis function of degree $p$, an interior knot can be repeated
at most $p$ times, and the boundary knots can be repeated at most $p+1$ times.
A knot vector for which the two boundary knots are repeated $p+1$ times is
said to be open. In this case, the basis functions are interpolatory at the
first and the last knot. Important properties of the B-Spline basis functions
include nonnegativity, partition of unity, local support and
$C^{p-k}$-continuity.

Higher dimensional B-Spline objects are defined using tensor products.

######  Definition 2.

A B-Spline curve $C(\xi)$ is defined by

(7) $C(\xi)=\sum^{n}_{i=1}P_{i}N^{p}_{i}(\xi),$

where $\\{P_{i}:i=1,\ldots,n\\}$ are the control points and $N^{p}_{i}$ are
B-Spline basis functions defined in (6).

######  Definition 3.

A B-Spline surface $S(\xi,\eta)$ is defined by

(8)
$S(\xi,\eta)=\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}N_{i,j}^{p_{1},p_{2}}(\xi,\eta)P_{i,j},$

where $P_{i,j}$, $i=1,2,\ldots,n_{1}$, $j=1,2,\ldots,n_{2}$, denote the
control points, $N_{i,j}^{p_{1},p_{2}}$ is the tensor product of B-Spline
basis functions $N_{i}^{p_{1}}$ and $N_{j}^{p_{2}}$, and
$\Xi_{1}=\\{\xi_{1},\xi_{2},\ldots,\xi_{n_{1}+p_{1}+1}\\}$ and
$\Xi_{2}=\\{\eta_{1},\eta_{2},\ldots,\eta_{n_{2}+p_{2}+1}\\}$ are the
corresponding knot vectors.

Similarly three dimensional B-Spline solids can be defined using two tensor
products.

Polynomials cannot exactly describe frequently encountered shapes in
engineering, particularly the conic family (e.g., a circle). While B-Splines
are flexible and have many nice properties for curve design, they are also
incapable of representing such curves exactly. Such limitations are overcome
by NURBS functions that can exactly represent a wide array of objects.

Rational representation of conics originates from projective geometry. The
“coordinates” in the additional dimension are called weights, which we shall
denote by $w$. Furthermore, let $\\{P^{w}_{i}\\}$ be a set of control points
for a projective B-Spline curve in $\mathbb{R}^{3}$. For the desired NURBS
curve in $\mathbb{R}^{2}$, the weights and the control points are derived by
the relations

(9) $w_{i}=(P^{w}_{i})_{3},\qquad(P_{i})_{d}=(P^{w}_{i})/w_{i},\quad d=1,2,$

where $w_{i}$ is called the $i^{\mathrm{th}}$ weight and $(P_{i})_{d}$ is the
$d^{\mathrm{th}}$-dimension component of the vector $P_{i}$. The weight
function $w(\xi)$ is defined as

(10) $w(\xi)=\sum_{i=1}^{n}N_{i}^{p}(\xi)w_{i}.$

We now formally define NURBS objects.

######  Definition 4.

The NURBS basis functions and curve are defined by

(11) $R_{i}^{p}(\xi)=\frac{N_{i}^{p}(\xi)w_{i}}{w(\xi)}\quad{\rm and}\quad
C(\xi)=\sum_{i=1}^{n}R_{i}^{p}(\xi)P_{i}.$

######  Definition 5.

The NURBS surfaces are defined by

(12)
$S(\xi,\eta)=\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}R_{i,j}^{p_{1},p_{2}}(\xi,\eta)P_{i,j},$

where $R_{i,j}^{p_{1},p_{2}}$ is the tensor product of NURBS basis functions
$R_{i}^{p_{1}}$ and $R_{j}^{p_{2}}$.

NURBS functions also satisfy the properties of B-Spline functions [39, 42,
47].

###  3.2 Derivatives of B-Splines

Derivatives of B-Splines [19] and their conditioning are very important for
the estimation of the condition number of the stiffness matrix. The recursive
definition of B-Spline functions allow us to seek the relationship between the
derivative of a given B-Spline basis function and lower degree basis function.

######  Definition 6.

The derivative of the $i^{th}$ B-Spline basis function defined in (6) is given
by

(13)
$\frac{d}{d\xi}N^{p}_{i}(\xi)=\frac{p}{\xi_{i+p}-\xi_{i}}N^{p-1}_{i}(\xi)-\frac{p}{\xi_{i+p+1}-\xi_{i+1}}N^{p-1}_{i+1}(\xi).$

By repeated differentiation of (13) we get the general formula for any order
derivative. Since we are only interested in the first derivative, we ignore
further details [13].

The derivatives of rational functions will clearly depend on the derivatives
of their non-rational counterpart. Definition 6 can be generalized for NURBS.

######  Definition 7.

The derivative of the $i^{th}$ NURBS basis function is given by

(14)
$\displaystyle\frac{d}{d\xi}R^{p}_{i}(\xi)=w_{i}\displaystyle\frac{w(\xi)\displaystyle\frac{d}{d\xi}N^{p}_{i}(\xi)-\displaystyle\frac{d}{d\xi}w(\xi)N^{p}_{i}(\xi)}{(w(\xi)^{2})}.$

where $w_{i}$ and $w(\xi)$ are defined in (9) and (10), respectively.

###  3.3 Condition number of B-Splines

In this section, we recall bounds for the condition number of B-Splines.

We need to know the bounds on B-Spline basis functions in some $L_{s}$-norm,
where $s\in[1,\infty]$. We estimate the size of the coefficients of a
polynomial of degree $p$ in two dimensions when it is represented using the
tensor product structure of B-Spline basis functions. The condition number of
a basis can be defined as follows.

######  Definition 8.

A basis $\\{N_{i}\\}$ of a normed linear space is said to be stable with
respect to a vector norm if there are constants $K_{1}$ and $K_{2}$ such that
for all coefficients $\\{v_{i}\\}$ the following relation holds:

(15) $K_{1}^{-1}\|\\{v_{i}\\}\|\leq\Big{\|}\sum_{i}v_{i}N_{i}\Big{\|}\leq
K_{2}\|\\{v_{i}\\}\|.$

The number $\kappa=K_{1}K_{2}$, with $K_{1}$ and $K_{2}$ as small as possible,
is called the condition number of $\\{N_{i}\\}$ with respect to $\|\cdot\|$.
Note that we use the symbols $\|\cdot\|$ and $\|\\{\cdot\\}\|$ for the norms
in the vector space and the vector norm, respectively.

Such condition numbers give an upper bound for magnification of the error in
the coefficients to the function values. Indeed, if $\displaystyle
f=\sum_{i}f_{i}N_{i}\neq 0$ and $\displaystyle g=\sum_{i}g_{i}N_{i}$, then it
follows immediately from (15) that

$\displaystyle\frac{\|f-g\|}{\|f\|}\leq\kappa\frac{\|\\{f_{i}-g_{i}\\}\|}{\|\\{f_{i}\\}\|}.$

More details on the approximation properties and the stability of B-Splines
can be found in [24, 25, 26, 30, 31, 36, 38]. We use these estimates of
$\kappa$ to estimate the bounds of the condition number of the stiffness
matrix and the mass matrix.

It is of central importance for working with B-Spline basis functions that its
condition number is bounded independently of the underlying knot sequence.
That is, the condition number of B-Splines does not depend on the multiplicity
of the knots of knot vector [10, 11, 12, 15]. In [11] is a direct estimate
that the worst condition number of a B-Spline of degree $p$ with respect to
any $L_{s}$-norm is bounded above by $p9^{p}$. It is also conjectured that the
real value of $\kappa$ grows like $2^{p}$, which is superior to the direct
estimate:

(16a) $\displaystyle\kappa$ $\displaystyle<p9^{p}$
$\displaystyle\quad\text{(direct estimate)},$ (16b) $\displaystyle\kappa$
$\displaystyle\sim 2^{p}$ $\displaystyle\quad\text{(de Boor's conjecture)}.$

In [14], the exact condition number of a B-Spline basis is shown to be
difficult to determine.

Scherer and Shadrin [44] proved that the upper bound of the condition number
$kappa$ of a B-Spline of degree $p$ with respect to $L_{s}$-norm is bounded by

(17) $\kappa<\displaystyle p^{\frac{1}{2}}4^{p},$

which is closer to the conjecture in (16b). Scherer and Shadrin [45] proved
the following result.

######  Lemma 9.

For all $p$ and all $s\in[1,\infty]$,

(18) $\kappa<p2^{p}.$

Lemma 9 confirms the conjecture (16b) up to a polynomial factor. Possible
approaches to eliminate the polynomial factor are also discussed in [45].
Lemma 9 can be easily generalized to $d$-dimensions.

######  Lemma 10.

Using a tensor product B-Spline basis of degree $p$ in $d$-dimensions and
(18), the following is immediate:

(19) $\kappa<(p2^{p})^{d}.$

##  4 Estimates of condition number

In this section, we give estimates for the condition number of the stiffness
matrix (in Section 4.1) and the mass matrix (in Section 4.2) obtained from
isogeometric discretization. In each case, we have bounds on the condition
number with respect to both $h$\- and $p$-refinements. For $h$-refinement,
upper bounds for the maximum eigenvalues, a lower bound for the minimum
eigenvalue, and an upper bound for the condition number are given. For
$p$-refinement, we prove upper and lower bounds for the maximum eigenvalue,
lower bounds for the minimum eigenvalue, and upper bounds for the condition
number.

###  4.1 Stiffness matrix

In this section, we give estimates for the condition number of the stiffness
matrix with estimates for $h$-refinement in Section 4.1.1 and for
$p$-refinement in Section 4.1.2.

####  4.1.1 $h$-refinement

Without loss of generality, we begin with a two-dimensional open parametric
domain $\Omega=(0,1)^{2}$ that we refer to as a patch. Given two open knot
vectors $\Xi_{1}=\\{0=\xi_{1},\xi_{2},\xi_{3},\ldots,\xi_{m_{1}}=1\\}$ and
$\Xi_{2}=\\{0=\eta_{1},\eta_{2},\eta_{3},\ldots,\eta_{m_{2}}=1\\}$, we
partition the patch $\Omega$ into a mesh

$\mathcal{Q}_{h}=\\{Q=(\xi_{i},\xi_{i+1})\otimes(\eta_{j},\eta_{j+1}),i=p_{1}+1,2,\ldots,m_{1}-p_{1}-1,j=p_{2}+1,2,\ldots,m_{2}-p_{2}-1\\},$

where $Q$ is a two-dimensional open knot-span whose diameter is denoted by
$h_{Q}$. We consider a family of quasi-uniform meshes
$\\{\mathcal{Q}_{h}\\}_{h}$ on $\Omega$, where
$h=\max\\{h_{Q}|Q\in\mathcal{Q}_{h}\\}$ denotes the family index [5]. Let
$\mathcal{S}_{h}$ denote the B-spline space associated with the mesh
$\mathcal{Q}_{h}$. Given two adjacent elements $Q_{1}$ and $Q_{2}$, we denote
by $m_{Q_{1}Q_{2}}$ the number of continuous derivatives across their common
face $\partial Q_{1}\cap\partial Q_{2}$. In the analysis, we will use the
usual Sobolev space of order $m\in\mathbb{N}$,

(20) $\displaystyle\mathcal{H}^{m}(\Omega)=$ $\displaystyle\Big{\\{}v\in
L^{2}(\Omega)\mathrm{~{}such~{}that~{}}v|_{Q}\in H^{m}(Q),\forall
Q\in\mathcal{Q}_{h},\mathrm{~{}and}$
$\displaystyle\nabla^{i}(v|_{Q_{1}})=\nabla^{i}(v|_{Q_{2}})\mathrm{~{}on~{}}\partial
Q_{1}\cap\partial Q_{2},$ $\displaystyle\forall
i\in\mathbb{N}\mathrm{~{}with~{}}0\leq
i\leq\min\\{m_{Q_{1}Q_{2}},m-1\\},\forall
Q_{1},Q_{2}\mathrm{~{}with~{}}\partial Q_{1}\cap\partial
Q_{2}\neq\emptyset\Big{\\}},$

where $\nabla^{i}$ has the usual meaning of $i^{\mathrm{th}}$-order partial
derivative. The space $\mathcal{H}^{m}$ is equipped with the following semi-
norms and norm

$|v|_{\mathcal{H}^{i}(\Omega)}^{2}=\sum_{Q\in\mathcal{Q}_{h}}|v|_{H^{i}(Q)}^{2},~{}~{}0\leq
i\leq m,\quad{\rm
and}\quad\|v\|_{\mathcal{H}^{m}(\Omega)}^{2}=\sum_{i=0}^{m}|v|_{\mathcal{H}^{i}(\Omega)}^{2}.$

On a regular mesh of size $h$, the condition number of the finite element
equations for a second-order elliptic boundary value problem can be obtained
using inverse estimates [2, 7, 8]. Similar inverse estimates are of interest
for the isogeometric framework using NURBS basis functions.

To keep the article self-contained, we recall some results from [5, 47].

######  Theorem 11.

Let $\mathcal{S}_{h}$ be the spline space consisting of piecewise polynomials
of degree $p$ associated with uniform partitions. Then there exists a constant
$C=C(shape)$, such that for all $0\leq l\leq m$,

(21) $\|v\|_{\mathcal{H}^{m}(\Omega)}\leq
Ch^{l-m}\|v\|_{\mathcal{H}^{l}(\Omega)},\quad\forall v\in\mathcal{S}_{h}.$

The proof of the above theorem, for a particular case $m=2$ and $l=1$, is
given in [5]. More general inverse inequalities can be easily derived
following the same approach. By taking $m=1$ and $l=0$, the following can be
easily derived from (21)

(22) $a(v,v)=\int_{\Omega}|\nabla v|^{2}\leq Ch^{-2}\|v\|^{2}.$

Under suitable conditions the condition number related to elliptic problems in
finite element analysis scales as $h^{-2}$ [18, 29, 48]. We prove the similar
result for the stiffness matrix arising in isogeometric discretization.

We first prove

######  Lemma 12.

There exist constants $C_{1}$ and $C_{2}$ independent of $h$ (but may depend
on $p$), such that for all
$v=\displaystyle\sum_{i=1}^{n_{h}}v_{i}N_{i}\in\mathcal{S}_{h}$,

(23)
$C_{1}h^{2}\|\\{v_{i}\\}\|^{2}\leq\Big{\|}\sum_{i=1}^{n_{h}}v_{i}N_{i}\Big{\|}^{2}\leq
C_{2}h^{2}\|\\{v_{i}\\}\|^{2}.$

###### Proof.

We only consider the non-trivial case: there exists some $i$ for which
$v_{i}\neq 0$. For any $Q\in\mathcal{Q}_{h}$, there are $(p+1)^{2}$ basis
functions with non-zero support. Let
$\mathcal{I}_{h}^{Q}\equiv\\{i^{Q}_{1},i^{Q}_{2},\ldots,i^{Q}_{p+1}\\}\times\\{j^{Q}_{1},j^{Q}_{2},\ldots,j^{Q}_{p+1}\\}\subset\\{1,2,\ldots,n_{h}\\}$
denote the index set for the basis functions that have non-zero support in
$Q$. Also, let
$\bar{v}_{q}={\displaystyle\max_{i\in\mathcal{I}_{h}^{Q}}|v_{i}|}$ and
$\bar{v}={\displaystyle\max_{i=1,2,\ldots,n_{h}}|v_{i}|}$. Now using
positivity and partition of unity properties of basis functions, the right
hand side inequality can be proved as follows:

$\displaystyle\|v\|^{2}=$
$\displaystyle\sum_{Q\in\mathcal{Q}_{h}}\int_{Q}v^{2}=\sum_{Q\in\mathcal{Q}_{h}}\int_{Q}\Bigg{(}\sum_{i\in\mathcal{I}_{h}^{Q}}v_{i}N_{i}\Bigg{)}^{2}\leq\sum_{Q\in\mathcal{Q}_{h}}\int_{Q}\Bigg{(}\bar{v}_{q}\sum_{i\in\mathcal{I}_{h}^{Q}}N_{i}\Bigg{)}^{2}$
$\displaystyle\leq$
$\displaystyle\sum_{Q\in\mathcal{Q}_{h}}\int_{Q}\bar{v}_{q}^{2}\leq\sum_{Q\in\mathcal{Q}_{h}}h_{Q}^{2}\bar{v}^{2}_{q}\leq\sum_{Q\in\mathcal{Q}_{h}}h_{Q}^{2}\sum_{i\in\mathcal{I}_{h}^{Q}}v_{i}^{2}$
$\displaystyle\leq$ $\displaystyle
h^{2}\sum_{Q\in\mathcal{Q}_{h}}\sum_{i\in\mathcal{I}_{h}^{Q}}v_{i}^{2}\leq
C_{2}h^{2}\sum_{i=1}^{n_{h}}v_{i}^{2}=C_{2}h^{2}\|\\{v_{i}\\}\|^{2}.$

For the left hand side inequality,

$\displaystyle h^{2}\|\\{v_{i}\\}\|^{2}=$ $\displaystyle
h^{2}\sum_{i=1}^{n_{h}}v_{i}^{2}\leq
h^{2}\sum_{i=1}^{n_{h}}\bar{v}^{2}=h^{2}{n_{h}}\bar{v}^{2}\leq
h^{2}\bigg{(}\frac{C}{h}\bigg{)}^{2}\bar{v}^{2}=C^{2}\bar{v}^{2}$
$\displaystyle=$ $\displaystyle C^{2}\|\\{v_{i}\\}\|_{L_{\infty}}^{2}\leq
C^{2}K_{1}^{2}\|v\|_{L_{\infty}}^{2}\text{ $\Big{(}$using
(\ref{eq:BasisCondP}), $K_{1}^{-1}\|\\{v_{i}\\}\|_{L_{\infty}}\leq\big{\|}\sum
v_{i}N_{i}\big{\|}_{L_{\infty}}\Big{)}$ }$ $\displaystyle\leq$ $\displaystyle
C^{2}K_{1}^{2}\|v\|^{2}.$

The result then follows by taking $\displaystyle
C_{1}=\left(\displaystyle\frac{1}{C^{2}K_{1}^{2}}\right).$ ∎

We now turn to the problem of obtaining bounds on the extremal eigenvalues and
the condition number.

######  Theorem 13.

Let $A$ be the stiffness matrix $A=(a_{ij})$, where
$a_{ij}=a(N_{i},N_{j})=\displaystyle\int_{\Omega}\nabla N_{i}\cdot\nabla
N_{j}$. Then the upper bound on $\lambda_{\text{max}}$ and lower bound on
$\lambda_{\text{min}}$ are given by

$\lambda_{\text{max}}\leq c_{1}$ and $\lambda_{\text{min}}\geq c_{2}h^{2}$,

where $c_{1},c_{2}$ are constants independent of $h$. The bound on $\kappa(A)$
is given by

$\kappa(A)\leq Ch^{-2}$,

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

###### Proof.

Let $v=\displaystyle\sum_{i=1}^{n_{h}}v_{i}N_{i}$. Then
$a(v,v)=\\{v_{i}\\}\cdot A\\{v_{i}\\}$, where
$\\{v_{i}\\}=\\{v_{1},v_{2},\ldots,v_{n_{h}}\\}$. Using the inverse estimate
(22),

$\displaystyle\frac{\\{v_{i}\\}\cdot A\\{v_{i}\\}}{\|\\{v_{i}\\}\|^{2}}=$
$\displaystyle\frac{a(v,v)}{\|\\{v_{i}\\}\|^{2}}\leq\frac{Ch^{-2}\|v\|^{2}}{\|\\{v_{i}\\}\|^{2}}.$

Using (23),

$\displaystyle\displaystyle\frac{Ch^{-2}\|v\|^{2}}{\|\\{v_{i}\\}\|^{2}}\leq\frac{Ch^{-2}C_{2}h^{2}\|\\{v_{i}\\}\|^{2}}{\|\\{v_{i}\\}\|^{2}}=CC_{2}=c_{1}.$

Hence,

(24) $\displaystyle\lambda_{\text{max}}=\sup_{\begin{subarray}{c}v\neq
0\end{subarray}}\frac{\\{v_{i}\\}\cdot A\\{v_{i}\\}}{\|\\{v_{i}\\}\|^{2}}\leq
c_{1}.$

On the other hand, for the bounds on $\lambda_{\text{min}}$, by using
coercivity of bilinear form $a(v,v)$,

$\displaystyle\frac{\\{v_{i}\\}\cdot A\\{v_{i}\\}}{\|\\{v_{i}\\}\|^{2}}=$
$\displaystyle\frac{a(v,v)}{\|\\{v_{i}\\}\|^{2}}\geq\frac{\alpha\|v\|_{H^{1}}^{2}}{\|\\{v_{i}\\}\|^{2}}\geq\frac{\alpha\|v\|^{2}}{\|\\{v_{i}\\}\|^{2}}.$

Using (23) again,

$\displaystyle\displaystyle\frac{\alpha\|v\|^{2}}{\|\\{v_{i}\\}\|^{2}}\geq\frac{\alpha_{1}C_{1}h^{2}\|\\{v_{i}\\}\|^{2}}{\|\\{v_{i}\\}\|^{2}}=$
$\displaystyle\alpha_{1}C_{1}h^{2}=c_{2}h^{2}.$

Hence,

(25) $\displaystyle\lambda_{\text{min}}=\inf_{\begin{subarray}{c}v\neq
0\end{subarray}}\frac{\\{v_{i}\\}\cdot A\\{v_{i}\\}}{\|\\{v_{i}\\}\|^{2}}\geq
c_{2}h^{2}.$

The condition number of the stiffness matrix is given by

$\kappa(A)=\frac{\lambda_{\text{max}}}{\lambda_{\text{min}}},\text{where}~{}\lambda_{\text{max}}=\max_{\begin{subarray}{c}v\neq
0\end{subarray}}\frac{\\{v_{i}\\}\cdot
A\\{v_{i}\\}}{\|\\{v_{i}\\}\|^{2}},~{}\text{and}~{}\lambda_{\text{min}}=\min_{\begin{subarray}{c}v\neq
0\end{subarray}}\frac{\\{v_{i}\\}\cdot A\\{v_{i}\\}}{\|\\{v_{i}\\}\|^{2}}.$

From (24) and (25),

(26) $\kappa(A)\leq Ch^{-2}.$

∎

####  4.1.2 $p$-refinement

In this section for $p$-refinement, we prove upper and lower bounds for the
maximum eigenvalue, lower bounds for the minimum eigenvalue, and upper bounds
for the condition number.

Let $\mathcal{S}_{p}$ be the tensor product space of spline functions of
degree $p$.

The following lemma is well known generalization of a theorem of Markov due to
Hill, Szechuan and Tamarkin [6, 37].

######  Lemma 14 (Schmidt’s inequality).

There exists a constant $C$ independent of $p$ such that for any polynomial
$f(x)$ of degree $p$,

(27) $\int_{-1}^{1}(f^{\prime}(x))^{2}dx\leq Cp^{4}\int_{-1}^{1}(f(x))^{2}dx.$

Note: No such constant $C$ exists so that (27) holds for all $f(x)$ with the
exponent smaller than 4.

Let $I=(-1,1)$. Using (27),

(28) $\int_{I}\left(\frac{dN_{p}(\xi)}{d\xi}\right)^{2}d\xi\leq
Cp^{4}\int_{I}(N_{p}(\xi))^{2}d\xi.$

Using (28),

(29) $\begin{split}\int_{\Omega}\nabla N_{p}(\xi,\eta)\cdot\nabla
N_{p}(\xi,\eta)d\xi
d\eta&=\displaystyle\int_{I}\int_{I}\left[\left(\frac{\partial
N_{p}(\xi,\eta)}{\partial\xi}\right)^{2}+\left(\frac{\partial
N_{p}(\xi,\eta)}{\partial\eta}\right)^{2}\right]d\xi d\eta\\\
&\leq\displaystyle Cp^{4}\int_{I\times I}(N_{p}(\xi,\eta))^{2}d\xi
d\eta.\end{split}$

Moreover, the following estimate directly follows from Schmidt’s inequality
and (29):

(30) $a(v,v)=\int_{\Omega}|\nabla v|^{2}\leq Cp^{4}\|v\|^{2}.$

We now have a similar result like Lemma 12 for the $p$-refinement.

######  Lemma 15.

There exist constants $C_{1}$ and $C_{2}$ independent of $p$ such that for all
$\displaystyle v=\sum_{i=1}^{n_{p}}v_{i}N_{i}\in\mathcal{S}_{p}$,

(31)
$\displaystyle\frac{C_{1}}{(p^{2}4^{p})^{2}}\|\\{v_{i}\\}\|^{2}\leq\displaystyle\Big{\|}\sum_{i=1}^{n_{p}}v_{i}N_{i}\Big{\|}^{2}\leq
C_{2}\|\\{v_{i}\\}\|^{2}.$

###### Proof.

From the stability of B-Splines there exists a constant $\gamma$ that depends
on the degree $p$ such that

(32)
$\displaystyle\Big{\|}\sum_{i=1}^{n_{p}}v_{i}N_{i}\Big{\|}\leq\|\\{v_{i}\\}\|\leq\gamma\Big{\|}\sum_{i=1}^{n_{p}}v_{i}N_{i}\Big{\|},$

From (19), $\gamma=p^{2}4^{p}$. In the estimate (31), the right hand side
inequality follows easily from nonnegativity and the partition of unity
properties of basis functions. The left hand side inequality follows from
(32). ∎

For the $p$-refinement of isogeometric discretization, the analog to Theorem
13 is

######  Theorem 16.

Let $\\{N_{i}\\}$ be a set of basis functions of $\mathcal{S}_{p}$ on a unit
square. Then

$\kappa(A)\leq Cp^{8}16^{p}$.

###### Proof.

We prove this theorem following the same approach as for the $h$-refinement
estimates. Let $v=\displaystyle\sum_{i=1}^{n_{p}}v_{i}N_{i}$, where
$\\{v_{i}\\}=\\{v_{1},v_{2},\ldots,v_{n_{p}}\\}$. Now using (30) and (31),

$\begin{split}\frac{\\{v_{i}\\}\cdot
A\\{v_{i}\\}}{\|\\{v_{i}\\}\|^{2}}=\frac{a(v,v)}{\|\\{v_{i}\\}\|^{2}}\leq\frac{Cp^{4}\|v\|^{2}}{\|\\{v_{i}\\}\|^{2}}\leq\frac{Cp^{4}C_{2}\|\\{v_{i}\\}\|^{2}}{\|\\{v_{i}\\}\|^{2}}=CC_{2}p^{4}=Cp^{4}.\end{split}$

Hence,

(33) $\lambda_{\text{max}}=\max_{\begin{subarray}{c}v\neq
0\end{subarray}}\frac{\\{v_{i}\\}\cdot A\\{v_{i}\\}}{\|\\{v_{i}\\}\|^{2}}\leq
Cp^{4}.$

To prove the lower bound for $\lambda_{\text{min}}$ we use (31) and coercivity
of bilinear form,

$\begin{split}\frac{\\{v_{i}\\}\cdot
A\\{v_{i}\\}}{\|\\{v_{i}\\}\|^{2}}&=\frac{a(v,v)}{\|\\{v_{i}\\}\|^{2}}\geq\frac{\alpha\|v\|_{H^{1}}^{2}}{\|\\{v_{i}\\}\|^{2}}\geq\frac{\alpha\|v\|^{2}}{\|\\{v_{i}\\}\|^{2}}\\\
&\geq\frac{\alpha\displaystyle\frac{C_{1}}{(p^{2}4^{p})^{2}}\|\\{v_{i}\\}\|^{2}}{\|\\{v_{i}\\}\|^{2}}={\frac{\alpha
C_{1}}{(p^{2}4^{p})^{2}}}={\frac{C}{(p^{4}16^{p})}}.\end{split}$

Hence,

(34) $\lambda_{\text{min}}=\min_{\begin{subarray}{c}v\neq
0\end{subarray}}\frac{\\{v_{i}\\}\cdot
A\\{v_{i}\\}}{\|\\{v_{i}\\}\|^{2}}\geq\frac{C}{(p^{4}16^{p})}.$

From (33) and (34),

$\displaystyle\kappa(A)=\frac{\lambda_{\text{max}}}{\lambda_{\text{min}}}\leq\frac{Cp^{4}}{{\left(\displaystyle\frac{C}{(p^{4}16^{p})}\right)}}\leq
C(p^{8}16^{p}).$

∎

######  Remark 17.

Theorem 16 can be easily generalized for higher dimensions. The bound for the
condition number of the stiffness matrix for a $d$-dimensional problem is
given by $(p^{4+2d}4^{pd})$.

While we proved an upper bound on the maximum eigenvalue of the stiffness
matrix using B-Spline basis functions, Theorem 16 is independent of the choice
of the basis functions (it holds for all kind of basis functions, not just
spline functions). From numerical experiments using B-Spline basis functions
(see Table 1), we observe that $\lambda_{\text{max}}$ depends linearly on the
polynomial degree $p$, which motivates further investigations.

Table 1:  Maximum eigenvalue of the stiffness matrix $A$ $p$ $h^{-1}$ | 1 | 2 | 4 | 8 | 16 | 32 | 64 | 128
---|---|---|---|---|---|---|---|---
2 | 0.36 | 1.42 | 1.42 | 1.49 | 1.50 | 1.50 | 1.50 | 1.50
3 | 0.45 | 1.04 | 1.37 | 1.52 | 1.56 | 1.57 | 1.57 | 1.57
4 | 0.41 | 0.94 | 1.33 | 1.72 | 1.81 | 1.83 | 1.84 | 1.84
5 | 0.35 | 0.88 | 1.32 | 1.93 | 2.10 | 2.14 | 2.14 | 2.14
6 | 0.34 | 0.85 | 1.32 | 2.12 | 2.40 | 2.46 | 2.47 | 2.47
7 | 0.33 | 0.84 | 1.32 | 2.26 | 2.70 | 2.78 | 2.80 | 2.80
8 | 0.32 | 0.83 | 1.33 | 2.36 | 2.99 | 3.11 | 3.13 | 3.14
9 | 0.31 | 0.82 | 1.33 | 2.43 | 3.29 | 3.44 | 3.47 | 3.47
10 | 0.31 | 0.82 | 1.34 | 2.47 | 3.56 | 3.77 | 3.80 | 3.81
20 | 0.29 | 0.78 | 1.36 | 2.65 | 5.02 | 6.95 | 7.20 | 7.23
30 | 0.29 | 0.78 | 1.36 | 2.69 | 5.28 | 9.38 | 10.55 | 10.66

The lower bound on the minimum eigenvalue depends on the stability of the
B-Spline basis functions, which cannot be improved further (especially beyond
the de Boor’s conjecture). On the other hand, the upper bound on the maximum
eigenvalue directly depends on the upper bound of the bilinear form $a(v,v)$.
We can improve the bound for $a(v,v)$ given in (30). In the following theorem
we improve this bound and provide our main result.

######  Theorem 18.

For the two dimensional problem the improved upper bound for the condition
number of the stiffness matrix $A$ is given by

(35) $\kappa(A)\leq Cp^{2}(p^{2}4^{p})^{2}=Cp^{6}16^{p}.$

The bound for a $d$-dimensional problem is given by

(36) $\kappa(A)\leq Cp^{2d+2}4^{pd}.$

For the sake of clarity we will give the proof of Theorem 18 in parts in
Lemmas 19-23.

It is clear from Table 1 that the maximum eigenvalue of the stiffness matrix
is independent of $p$ for the coarsest mesh size $h=1$, and linearly dependent
of $p$ asymptotically. In the analysis, we consider two dimensional problem on
the coarsest mesh first and extend it to finer meshes later. On the coarsest
mesh we have B-Spline basis functions of degree $p$ in one variable $\xi$,

$N^{p}_{i,\xi}=(-1)^{i}{p\choose i}(\xi-1)^{p-i}\xi^{i},\quad
i=0,1,2,\ldots,p.$

Similarly in variable $\eta$,

$N^{p}_{j,\eta}=(-1)^{j}{p\choose j}(\eta-1)^{p-j}\eta^{j},\quad
j=0,1,2,\ldots,p.$

Two variable B-Spline basis functions on the coarsest mesh is given by the
tensor product

$N^{p,p}_{i,j,\xi,\eta}=(-1)^{i+j}{p\choose i}{p\choose
j}\xi^{i}\eta^{j}(\xi-1)^{p-i}(\eta-1)^{p-j},\quad i,j=0,1,2,\ldots,p.$

We construct an upper bound of the diagonal entries of the stiffness matrix on
the coarsest mess i.e. single element stiffness matrix $A^{e}$.

######  Lemma 19.

There exists a constant $C$ independent of $p$, such that

(37)
$A^{e}_{(i,j),(i,j)}=a(N^{p,p}_{i,j,\xi,\eta},N^{p,p}_{i,j,\xi,\eta})=\int_{0}^{1}\int_{0}^{1}\nabla
N^{p,p}_{i,j,\xi,\eta}\cdot\nabla N^{p,p}_{i,j,\xi,\eta}d\xi d\eta\leq C.$

###### Proof.

We provide the major points of the proof. Some of the details can be found in
the research report [21]. For all $i,j=0,1,2,\ldots,p$,

$\begin{split}&a(N^{p,p}_{i,j,\xi,\eta},N^{p,p}_{i,j,\xi,\eta})=\int_{0}^{1}\int_{0}^{1}\nabla
N^{p,p}_{i,j,\xi,\eta}\cdot\nabla N^{p,p}_{i,j,\xi,\eta}d\xi d\eta\\\
&={p\choose i}^{2}{p\choose
j}^{2}\int_{0}^{1}\int_{0}^{1}\big{\\{}i\xi^{i-1}\eta^{j}(\xi-1)^{p-i}(\eta-1)^{p-j}+\\\
&\hskip
142.26378pt(p-i)\xi^{i}\eta^{j}(\xi-1)^{p-i-1}(\eta-1)^{p-j}\big{\\}}^{2}d\xi
d\eta\\\ &+{p\choose i}^{2}{p\choose
j}^{2}\int_{0}^{1}\int_{0}^{1}\big{\\{}j\xi^{i}\eta^{j-1}(\xi-1)^{p-i}(\eta-1)^{p-j}+\\\
&\hskip
142.26378pt(p-j)\xi^{i}\eta^{j}(\xi-1)^{p-i}(\eta-1)^{p-j-1}\big{\\}}^{2}d\xi
d\eta\\\ &\equiv I+II.\end{split}$

Now,

$\begin{split}I&={p\choose i}^{2}{p\choose
j}^{2}\int_{0}^{1}\int_{0}^{1}\left(i^{2}\xi^{2(i-1)}\eta^{2j}(\xi-1)^{2(p-i)}(\eta-1)^{2(p-j)}\right)d\xi
d\eta+\\\ &\quad{p\choose i}^{2}{p\choose
j}^{2}\int_{0}^{1}\int_{0}^{1}\left((p-i)^{2}\xi^{2i}\eta^{2j}(\xi-1)^{2(p-i-1)}(\eta-1)^{2(p-j)}\right)d\xi
d\eta+\\\ &\quad{p\choose i}^{2}{p\choose
j}^{2}\int_{0}^{1}\int_{0}^{1}\left(2i(p-i)\xi^{2i-1}\eta^{2j}(\xi-1)^{2p-2i-1}(\eta-1)^{2(p-j)}\right)d\xi
d\eta\\\ &\equiv I_{1}+I_{2}+I_{3}.\end{split}$

After simplifying (using results on factorial functions), we get

$\begin{split}I_{1}&\leq\frac{1}{2}{p\choose i}^{2}{p\choose
j}^{2}\frac{(2i)!(2p-2i)!}{(2p-1)!}\frac{(2j)!(2p-2j)!}{(2p+1)!},\\\
I_{2}&\leq\frac{1}{2}{p\choose i}^{2}{p\choose
j}^{2}\frac{(2i)!(2p-2i)!}{(2p-1)!}\frac{(2j)!(2p-2j)!}{(2p+1)!},\\\
I_{3}&=-\frac{1}{2}{p\choose i}^{2}{p\choose
j}^{2}\frac{(2i)!(2p-2i)!}{(2p-1)!}\frac{(2j)!(2p-2j)!}{(2p+1)!}.\end{split}$

For all $i=0,1,2,\ldots,p$,

$\begin{split}I&=I_{1}+I_{2}+I_{3}\\\ &=\begin{cases}I_{2},&{\text{ if
}i=0,}\\\ I_{1}+I_{2}+I_{3},&{\text{ if }i=1,2,\ldots,p-1,}\\\ I_{1},&{\text{
if }i=p,}\end{cases}\\\ &\leq\left\\{{p\choose
i}^{2}\frac{(2i)!(2p-2i)!}{(2p)!}\right\\}\left\\{{p\choose
j}^{2}\frac{(2j)!(2p-2j)!}{(2p)!}\right\\}=I_{a}I_{b},\text{ where}\\\
I_{a}&={p\choose
i}^{2}\frac{(2i)!(2p-2i)!}{(2p)!}=\frac{p!p!}{i!i!(p-i)!(p-i)!}\frac{(2i)!(2p-2i)!}{(2p)!},\\\
I_{b}&={p\choose
j}^{2}\frac{(2j)!(2p-2j)!}{(2p)!}=\frac{p!p!}{j!j!(p-j)!(p-j)!}\frac{(2j)!(2p-2j)!}{(2p)!}.\end{split}$

We prove that $I_{a}\leq C$ by induction on $p$, where $C$ is a constant
independent of $p$. For $p=1$, we have $I_{a}=1\text{ for all }i=0,1$. Hence,
the result holds for the base case. Assume that the result holds for
$p=m\text{ and for all }i=0,1,2,\ldots,m$,

(38)
$\displaystyle\frac{m!m!}{i!i!(m-i)!(m-i)!}\frac{(2i)!(2m-2i)!}{(2m)!}\leq C.$

Now we show that the result holds for $p=m+1~{}\text{and for
all}~{}i=0,1,2,\ldots,m+1$. We have

$\begin{split}&\frac{(m+1)!(m+1)!}{i!i!(m+1-i)!(m+1-i)!}\frac{(2i)!(2(m+1)-2i)!}{(2(m+1))!}\\\
&=\begin{cases}\displaystyle\left(\frac{m^{2}+2m+1}{4m^{2}+6m+2}\right)\left(\frac{4(m-i)^{2}+6(m-i)+2}{(m-i)^{2}+2(m-i)+1}\right)\left\\{\frac{m!m!}{i!i!(m-i)!(m-i)!}\frac{(2i)!(2m-2i)!}{(2m)!}\right\\},\\\
\hfill{\text{if }i=0,1,2,\ldots,m,}\\\ 1,{\text{if
}i=m+1.}\end{cases}\end{split}$

Using (38) and since
$\displaystyle\left(\frac{m^{2}+2m+1}{4m^{2}+6m+2}\right)\left(\frac{4(m-i)^{2}+6(m-i)+2}{(m-i)^{2}+2(m-i)+1}\right)\leq
1$, we get for all $i=0,1,2,\ldots,m+1$,

$\frac{(m+1)!(m+1)!}{i!i!(m+1-i)!(m+1-i)!}\frac{(2i)!(2(m+1)-2i)!}{(2(m+1))!}\leq
C.$

We now have $I_{a}\leq C$, where $C$ is a constant independent of $p$.
Similarly we can obtain that $I_{b}\leq C$. Hence,

$I=I_{a}I_{b}\leq C$.

Proceeding in the same way for ${II}$, we can prove that

${II}\leq C.$

Finally,

$\displaystyle a(N^{p,p}_{i,j,\xi,\eta},N^{p,p}_{i,j,\xi,\eta})=I+{II}\leq C.$

∎

Thus, we have proved that $a(N^{p,p}_{i,j,\xi,\eta},N^{p,p}_{i,j,\xi,\eta})$
is bounded by a constant independent of $p$. Since the upper bound of the
diagonal entries is the upper bound of all the entries of the stiffness
matrix, the maximum entry of the stiffness matrix is bounded by a constant
independent of $p$,

(39) $\displaystyle a(N^{p,q}_{i,j,\xi,\eta},N^{p,q}_{k,l,\xi,\eta})\leq C.$

Similarly, we can prove for three dimensional problem that

(40) $\displaystyle
a(N^{p,q,r}_{i,j,k,\xi,\eta,\zeta},N^{p,q,r}_{l,m,n,\xi,\eta,\zeta})\leq C.$

Using (39) and (40) we have

######  Lemma 20.

The maximum eigenvalue of the element stiffness matrix $A^{e}$ can be bounded
below by a constant $C$ independent of $p$,

$\lambda_{\text{max}}(A^{e})\geq C.$

###### Proof.

We prove this by using the basics of matrix norms. The max-norm of a matrix is
the element-wise norm defined by

$\|A^{e}\|_{\text{max}}=\ \text{max}\\{|a_{ij}|\\}.$

From (37),

$\text{max}\\{|a_{ij}|\\}=C,$

where $C$ is independent of $p$. By the equivalence of norms we have

$\|A^{e}\|_{2}\geq\|A^{e}\|_{\text{max}}=C.$

Hence,

$\lambda_{\text{max}}(A^{e})\geq C.$

∎

To bound $\lambda_{\text{max}}$ from above we bound the spectral norm by the
$\ell_{1}$-norm in

######  Lemma 21.

For any fixed $k$ and $l$ such that $0\leq k,l\leq p$ and for any $0\leq
i,j\leq p$,

$\displaystyle\sum_{i=0}^{p}\sum_{j=0}^{p}\big{|}a(N^{p,p}_{k,l,\xi,\eta},N^{p,p}_{i,j,\xi,\eta})\big{|}<C$,

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

###### Proof.

We again provide the main steps and for details refer the reader to the
research report [21]. We have

$N^{p,p}_{0,0,\xi,\eta}=(1-\xi)^{p}(1-\eta)^{p}.$

We first prove

$\displaystyle\sum_{i=0}^{p}\sum_{j=0}^{p}\big{|}a(N^{p,p}_{0,0,\xi,\eta},N^{p,p}_{i,j,\xi,\eta})\big{|}<C,$

where $C$ is a constant independent of $p$. We have

$\begin{split}&a(N^{p,p}_{0,0,\xi,\eta},N^{p,p}_{i,j,\xi,\eta})=\int_{0}^{1}\int_{0}^{1}\nabla
N^{p,p}_{0,0,\xi,\eta}\cdot\nabla N^{p,p}_{i,j,\xi,\eta}d\xi d\eta\\\
&=\int_{0}^{1}\int_{0}^{1}\left(\frac{\partial}{\partial\xi}N^{p,p}_{0,0,\xi,\eta}\frac{\partial}{\partial\xi}N^{p,p}_{i,j,\xi,\eta}\right)d\xi
d\eta+\int_{0}^{1}\int_{0}^{1}\left(\frac{\partial}{\partial\eta}N^{p,p}_{0,0,\xi,\eta}\frac{\partial}{\partial\eta}N^{p,p}_{i,j,\xi,\eta}\right)d\xi
d\eta\\\ &=I+II.\end{split}$

Now,

$\begin{split}I&=-p{p\choose i}{p\choose
j}\left(\int_{0}^{1}\eta^{j}(1-\eta)^{2p-j}d\eta\right)\\\
&\quad\quad\quad\quad\left(\int_{0}^{1}i\xi^{i-1}(1-\xi)^{2p-i-1}d\xi-\int_{0}^{1}(p-i)\xi^{i}(1-\xi)^{2p-i-2}d\xi\right)\\\
&I=\begin{cases}\displaystyle{p\choose
j}\frac{p^{2}}{(4p^{2}-1)}\frac{(j)!(2p-j)!}{(2p)!},\hfill\text{if }i=0,\\\
\displaystyle-p{p\choose i}{p\choose
j}\frac{2p}{(2p+1)}\frac{(i)!(2p-i)!}{(2p)!}\frac{(j)!(2p-j)!}{(2p)!}\frac{1}{(2p-i)}\left(1-\frac{(p-i)}{(2p-i-1)}\right),\\\
\hfill\text{if }i=1,2,\ldots,p-1,\\\ \displaystyle-{p\choose
j}\frac{2p}{(2p+1)}\frac{(p)!(p)!}{(2p)!}\frac{(j)!(2p-j)!}{(2p)!},\hfill\text{if
}i=p.\end{cases}\end{split}$

A similar expression can be obtained for $II$. We want to calculate

$\displaystyle\sum_{i=0}^{p}\sum_{j=0}^{p}\big{|}a(N^{p,p}_{0,0,\xi,\eta},N^{p,p}_{i,j,\xi,\eta})\big{|}$.

For $i=0$,

$\begin{split}\displaystyle\sum_{j=0}^{p}|I|&=\sum_{j=0}^{p}{p\choose
j}\displaystyle\frac{p^{2}}{(4p^{2}-1)}\frac{(j)!(2p-j)!}{(2p)!}<\frac{1}{3}\sum_{j=0}^{p}{p\choose
j}\displaystyle\frac{(j)!(2p-j)!}{(2p)!}\\\
&=\frac{1}{3}\sum_{j=0}^{p}\displaystyle\frac{(p)!(2p-j)!}{(2p)!(p-j)!}<\frac{1}{3}\left(1+\frac{1}{2}+\frac{1}{4}+\frac{1}{8}+\ldots+\frac{p!p!}{(2p)!}\right)<1.\end{split}$

For $i=1,2,\ldots,p-1$,

$\begin{split}&\displaystyle\sum_{i=1}^{p-1}\sum_{j=0}^{p}|I|\\\
&=\sum_{i=1}^{p-1}\sum_{j=0}^{p}{p\choose i}{p\choose
j}\displaystyle\frac{2p^{2}}{(2p+1)}\frac{(i)!(2p-i)!}{(2p)!}\frac{(j)!(2p-j)!}{(2p)!}\frac{1}{(2p-i)}\left(1-\frac{(p-i)}{(2p-i-1)}\right)\\\
&<\frac{1}{2}\displaystyle\sum_{i=1}^{p-1}\frac{p!(2p-i-2)!}{(2p-2)!(p-i)!}<\frac{1}{2}\left(\frac{1}{2}+\frac{1}{4}+\frac{1}{8}+\ldots+\frac{p!(p-1)!}{(2p-2)!}\right)<1.\end{split}$

For $i=p$,

$\begin{split}\displaystyle\sum_{j=0}^{p}|I|&=\sum_{j=0}^{p}{p\choose
j}\displaystyle\frac{2p}{(2p+1)}\frac{(p)!(p)!}{(2p)!}\frac{(j)!(2p-j)!}{(2p)!}<\sum_{j=0}^{p}\displaystyle\frac{p!p!p!(2p-j)!}{(2p)!(2p)!(p-j)!}\\\
&<\left(\frac{1}{2}+\frac{1}{4}+\frac{1}{8}+\ldots+\frac{(p!)^{4}}{((2p)!)^{2}}\right)<1.\end{split}$

Hence,

(41) $\displaystyle\sum_{i=0}^{p}\sum_{j=0}^{p}|I|<C,$

where $C$ is independent of $p$. Similarly, we have

(42) $\displaystyle\sum_{i=0}^{p}\sum_{j=0}^{p}|II|<C.$

Therefore, from (41) and (42),

(43)
$\displaystyle\sum_{i=0}^{p}\sum_{j=0}^{p}\big{|}a(N^{p,p}_{0,0,\xi,\eta},N^{p,p}_{i,j,\xi,\eta})\big{|}<C,$

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

We have bounded by a constant the absolute row sum for the first row of the
element stiffness matrix. Since on a uniform mesh the absolute row sum for all
rows of the element stiffness matrix are of the same order upto a constant, we
get the desired result:

$\displaystyle\sum_{i=0}^{p}\sum_{j=0}^{p}\big{|}a(N^{p,p}_{k,l,\xi,\eta},N^{p,p}_{i,j,\xi,\eta})\big{|}<C$,

for any fixed $k$ and $l$ such that $0\leq k,l\leq p$ and for any $0\leq
i,j\leq p$, where $C$ is a constant independent of $p$. ∎

Similar results can be obtained for higher dimensions. The next lemma is a
direct consequence of the Lemma 21 and gives an upper bound for the maximum
eigenvalue.

######  Lemma 22.

The maximum eigenvalue of the element stiffness matrix $A^{e}$ can be bounded
above by a constant $C$ that is independent of $p$:

$\lambda_{\text{max}}(A^{e})\leq C.$

###### Proof.

We have

$\displaystyle\sum_{i=0}^{p}\sum_{j=0}^{p}\big{|}a(N^{p,p}_{k,l,\xi,\eta},N^{p,p}_{i,j,\xi,\eta})\big{|}<C$,

where $C$ is a constant independent of $p$, which implies $\|A^{e}\|_{1}\leq
C$. Since $A^{e}$ is a symmetric matrix, we have
$\|A^{e}\|_{1}=\|A^{e}\|_{\infty}$. Therefore, using
$\|A^{e}\|_{2}\leq\sqrt{\|A^{e}\|_{1}\|A^{e}\|_{\infty}}$, we get

$\|A^{e}\|_{2}\leq|A^{e}\|_{1}\leq C$.

Thus,

$\lambda_{\text{max}}(A^{e})\leq C.$

∎

From Lemmas 20 and 22, for the element stiffness matrix we have

######  Lemma 23.

$\lambda_{\text{max}}(A^{e})=C$, where $C$ is a constant independent of $p$.

The results in Lemma 20 and Lemma 22 are proved for an element stiffness
matrix on a single element mesh. Obviously these results hold for all element
stiffness matrices on finer meshes. Therefore, Lemma 23 holds for all element
stiffness matrices on refined meshes.

Now we bound the spectral norm (maximum eigenvalue) of the global stiffness
matrix by its $\ell_{1}$-norm, i.e., the maximum of the row sum over all the
rows of the global stiffness matrix. The $\ell_{1}$-norm of global stiffness
matrix will depend on the $\ell_{1}$ norm of the element stiffness matrices
and on their assembly. Therefore, the bounds for maximum eigenvalue of the
global stiffness matrix can be expressed in terms of the maximum eigenvalues
of the corresponding element stiffness matrices and the maximum number of
overlaps within the rows and columns of element stiffness matrices (which is
the number of element stiffness matrices that contributes at a particular
nonzero position in the global matrix).

In the process of assembling the global stiffness matrix, the overlaps within
the element stiffness matrices depend on the regularity of the basis functions
used in the discretization. For $C^{0}$\- and $C^{p-1}$-continuous basis
fucntion the overlaps within the elements will be minimum and maximum,
respectively. It is easy to see that for $C^{p-1}$-continuous basis functions
the overlaps will be in $(p+1)^{2}$ knot spans (e.g., see Fig. 1).

Using the bound for maximum eigenvalue of element stiffness matrices we state
the following lemma for the bound for maximum eigenvalue of global stiffness
matrix.

######  Lemma 24.

The maximum eigenvalue of the global stiffness matrix
$\lambda_{\text{max}}(A)=Cp^{2}$, where $C$ is a constant independent of $p$.

###### Proof.

In the assembly of the element stiffness matrices, the maximum number of
overlaps for a particular nonzero position in the global matrix is
$(p+1)^{2}$. We have

$\begin{split}\|A\|_{1}&\leq(\text{maximum number of
overlaps})\times\|A^{e}\|_{1},\\\ &\leq((p+1)^{2})\times C.\end{split}$

The inequality $\|A\|_{2}\leq\sqrt{\|A\|_{1}\|A\|_{\infty}}$ and symmetry of
$A$ imply

$\|A\|_{2}\leq\|A\|_{1}\leq C(p+1)^{2}$.

Hence,

$\lambda_{\text{max}}(A)\leq Cp^{2}.$

∎

Now, using the bound for maximum eigenvalue given in Lemma 24, the proof of
Theorem 18 follows directly.

######  Remark 25.

The estimate for maximum eigenvalue given in Lemma 24 is not sharp. In reality
this estimate is not quadratic in $p$. This can be explained by the following
observation. In the overlapping all of the elements of element stiffness
matrices are not summed in absolute value. Some of the negative entries
overlap with positive entries which in result reduces the row sum of the
global stiffness matrix. From our numerical experiments, we conjecture the
following.

The maximum eigenvalue of the global stiffness matrix
$\lambda_{\text{max}}(A)=Cp$, where $C$ is a constant independent of $p$.

Fig. 1: Global stiffness matrix: Assembly of element stiffness matrices for
$p=4$ on $8\times 8$ spans

######  Remark 26.

We used the condition number of B-Splines $\kappa\sim p2^{p}$ and
$\lambda_{max}\sim p^{2}$ in reaching the above estimates. If we use the de
Boor’s conjecture (the condition number of B-Splines $\kappa\sim 2^{p}$) and
$\lambda_{max}\sim p$ (see Remark 25) instead, then the upper bound of the
stiffness matrix can be further improved and given by

(44) $\kappa(A)\leq Cp4^{pd}.$

###  4.2 Mass matrix

In this section, we give estimates for the condition number of the mass matrix
with estimates for $h$-refinement in Section 4.2.1 and for $p$-refinement in
Section 4.2.2.

####  4.2.1 $h$-refinement

Let $M=(m_{ij})$ be the mass matrix, where

$m_{ij}=(N_{i},N_{j})=\displaystyle\int_{\Omega}N_{i}N_{j}\quad\quad
i,j=1,2,\ldots,n_{h}$.

The following lemma gives estimates for the maximum and minimum eigenvalues of
the mass matrix with respect to $h$.

######  Lemma 27.

For the extremal eigenvalues of the mass matrix $M=(m_{ij})=(N_{i},N_{j})$,

$C_{1}h^{2}\leq\lambda_{\text{min}}\leq\lambda_{\text{max}}\leq C_{2}h^{2},$

where $C_{1},C_{2}$ are constants independent of $h$. Furthermore,

$c_{1}\leq\kappa(M)\leq c_{2},$

where $c_{1},c_{2}$ are constants independent of $h$.

###### Proof.

Using (23), we bound both the extremal eigenvalues of the mass matrix. For the
minimum eigenvalue,

$\frac{\\{v_{i}\\}\cdot
M\\{v_{i}\\}}{\|\\{v_{i}\\}\|^{2}}=\frac{(v,v)}{\|\\{v_{i}\\}\|^{2}}\geq\frac{C_{1}h^{2}\|\\{v\\}\|^{2}}{\|\\{v_{i}\\}\|^{2}}=C_{1}h^{2}.$

For the maximum eigenvalue,

$\frac{\\{v_{i}\\}\cdot
M\\{v_{i}\\}}{\|\\{v_{i}\\}\|^{2}}=\frac{(v,v)}{\|\\{v_{i}\\}\|^{2}}\leq\frac{C_{2}h^{2}\|\\{v\\}\|^{2}}{\|\\{v_{i}\\}\|^{2}}=C_{2}h^{2}.$

So,

$C_{1}h^{2}\leq\lambda_{\text{min}}\leq\lambda_{\text{max}}\leq C_{2}h^{2}.$

Hence,

$c_{1}\leq\kappa(M)\leq c_{2}.$

∎

####  4.2.2 $p$-refinement

In this section, we estimate the bounds on the extremal eigenvalues and the
condition number of the mass matrices for $p$-refinement.

######  Lemma 28.

The element mass matrix is a positive matrix and all of the entries of the
element mass matrix are bounded above by $\displaystyle\frac{C}{(2p+1)^{2}}$,
where $C$ is a constant independent of $p$.

###### Proof.

We have

$\begin{split}M^{e}_{(i,j),(k,l)}&=(N^{p,p}_{i,j,\xi,\eta},N^{p,p}_{k,l,\xi,\eta})=\int_{0}^{1}\int_{0}^{1}N^{p,p}_{i,j,\xi,\eta}\cdot
N^{p,p}_{k,l,\xi,\eta}d\xi d\eta\\\
&=\int_{0}^{1}\int_{0}^{1}\left((-1)^{i+j}{p\choose i}{p\choose
j}\xi^{i}\eta^{j}(\xi-1)^{p-i}(\eta-1)^{p-j}\right)\\\
&\quad\quad\quad\left((-1)^{k+l}{p\choose k}{p\choose
l}\xi^{k}\eta^{l}(\xi-1)^{p-k}(\eta-1)^{p-l}\right)d\xi d\eta\\\
&=\left(I\right)\left(II\right),\end{split}$

where

$\begin{split}I&={p\choose i}{p\choose
k}\left(\int_{0}^{1}\xi^{(i+k+1)-1}(1-\xi)^{(2p-i-k+1)-1}d\xi d\eta\right)\\\
&=\frac{p!p!}{i!k!(p-i)!(p-k)!}\frac{(i+k)!(2p-i-k)!}{(2p+1)!}\\\
&=\frac{1}{2p+1}{\left\\{\frac{p!p!}{i!k!(p-i)!(p-k)!}\frac{(i+k)!(2p-i-k)!}{(2p)!}\right\\}}=\frac{1}{2p+1}I_{1},\end{split}$

and

$\begin{split}II&={p\choose j}{p\choose
l}\left(\int_{0}^{1}\eta^{(j+l+1)-1}(1-\eta)^{(2p-j-l+1)-1}d\xi
d\eta\right)\\\
&=\frac{p!p!}{j!l!(p-j)!(p-l)!}\frac{(j+l)!(2p-j-l)!}{(2p+1)!}\\\
&=\frac{1}{2p+1}{\left\\{\frac{p!p!}{j!l!(p-j)!(p-l)!}\frac{(j+l)!(2p-j-l)!}{(2p)!}\right\\}}=\frac{1}{2p+1}II_{1}.\end{split}$

By induction on $p$ we easily obtain that (as we proved in Lemma 19),

$\begin{split}I_{1}=\left\\{\frac{p!p!}{i!k!(p-i)!(p-k)!}\frac{(i+k)!(2p-i-k)!}{(2p)!}\right\\}\leq
C.\end{split}$

Similarly, $II_{1}\leq C$. Therefore

(45) $M^{e}_{(i,j),(k,l)}\leq\frac{C}{(2p+1)^{2}}.$

It is also clear that for all $p\geq 1$ and $i,k=0,1,2,\ldots,p$, $I_{1}>0$,
and $II_{1}>0$. Hence, the mass matrix $M^{e}_{(i,j),(k,l)}$ is a positive
matrix. ∎

######  Lemma 29.

The maximum eigenvalue of the element mass matrix $M^{e}$ can be bounded below
by

$\lambda_{\text{max}}(M^{e})\geq\frac{C}{(2p+1)^{2}}.$

###### Proof.

Following the proof of Lemma 20 and (45) we get the desired result. ∎

To bound $\lambda_{\text{max}}$ from above we bound the spectral norm by the
$\ell_{1}$-norm of the mass matrix. In the following lemma we first compute
the $\ell_{1}$-norm of the mass matrix.

######  Lemma 30.

For the mass matrix $M^{e}$ on the coarsest mesh,

$\|M^{e}\|_{1}=\frac{1}{(p+1)^{2}}.$

###### Proof.

We have

$\begin{split}\displaystyle\|M\|_{1}&=\max_{\begin{subarray}{c}i,j\end{subarray}}\sum_{k,l}\int_{0}^{1}\int_{0}^{1}N_{i,j,\xi,\eta}^{p,p}\cdot
N_{k,l,\xi,\eta}^{p,p}d\xi
d\eta=\max_{\begin{subarray}{c}i,j\end{subarray}}\sum_{k,l}(N_{i,j,\xi,\eta}^{p,p},N_{k,l,\xi,\eta}^{p,p})\\\
&=\max_{\begin{subarray}{c}i,j\end{subarray}}(N_{i,j,\xi,\eta}^{p,p},\sum_{k,l}N_{k,l,\xi,\eta}^{p,p})=\max_{\begin{subarray}{c}i,j\end{subarray}}(N_{i,j,\xi,\eta}^{p,p},1)\quad\quad\Big{(}\text{since
}\sum_{k,l}N_{k,l,\xi,\eta}^{p,p}=1\Big{)}\\\
&=\max_{\begin{subarray}{c}i,j\end{subarray}}\int_{0}^{1}\int_{0}^{1}N_{i,j,\xi,\eta}^{p,p}d\xi
d\eta.\end{split}$

Now,

$\begin{split}&\int_{0}^{1}\int_{0}^{1}N_{i,j,\xi,\eta}^{p,p}d\xi
d\eta=\int_{0}^{1}\int_{0}^{1}(-1)^{i+j}{p\choose i}{p\choose
j}\xi^{i}\eta^{j}(\xi-1)^{p-i}(\eta-1)^{p-j}d\xi d\eta\\\ &={p\choose
i}{p\choose
j}\left(\int_{0}^{1}\xi^{(i+1)-1}(1-\xi)^{(p-i+1)-1}d\xi\right)\left(\int_{0}^{1}\eta^{(j+1)-1}(1-\eta)^{(p-j+1)-1}d\eta\right)\\\
&={p\choose i}{p\choose
j}\left(\frac{\Gamma(i+1)\Gamma(p-i+1)}{\Gamma(p+2)}\right)\left(\frac{\Gamma(j+1)\Gamma(p-j+1)}{\Gamma(p+2)}\right)\\\
&=\frac{p!}{i!(p-i)!}\frac{p!}{j!(p-j)!}\frac{i!(p-i)!}{(p+1)!}\frac{j!(p-j)!}{(p+1)!}=\frac{1}{(p+1)^{2}}.\end{split}$

Hence,

$\max_{\begin{subarray}{c}i,j\end{subarray}}\int_{0}^{1}\int_{0}^{1}N_{i,j,\xi,\eta}^{p,p}d\xi
d\eta=\frac{1}{(p+1)^{2}}.$

∎

The symmetry of $M^{e}$ implies

(46) $\|M^{e}\|_{\infty}=\|M^{e}\|_{1}=\frac{1}{(p+1)^{2}}.$

######  Lemma 31.

The maximum eigenvalue of the mass matrix $M^{e}$ on the coarsest mesh can be
bounded above by

$\lambda_{\text{max}}(M^{e})\leq C\frac{1}{(p+1)^{2}}.$

###### Proof.

We have the following inequality for matrix norms

$\|M^{e}\|_{2}^{2}\leq\|M^{e}\|_{1}\|M^{e}\|_{\infty}.$

Using Lemma 30 and (46) we get the bound on the spectral norm of $M^{e}$,

$\|M^{e}\|_{2}\leq C\frac{1}{(p+1)^{2}}.$

∎

######  Remark 32.

In fact, for the coarsest mesh we get
$\lambda_{\text{max}}(M^{e})=\displaystyle\frac{1}{(p+1)^{2}}$ by Lemma 28 and
by [49, Lemma 2.5].

Using the same argument as in the stiffness matrix case we can give the
estimate for the maximum eigenvalue of the global mass matrix using the
estimates for the element mass matrices.

######  Lemma 33.

The maximum eigenvalue of the global mass matrix $M$ can be bounded above by

$\lambda_{\text{max}}(M)\leq C,$

where $C$ is a constant independent of $p$ (may depend on $h$).

###### Proof.

Following the proof of Lemma 24, we have

$\begin{split}\lambda_{\text{max}}(M)&\leq
C((p+1)^{2})\times\displaystyle\frac{1}{(p+1)^{2}}.\end{split}$

Hence,

$\lambda_{\text{max}}(M)\leq C.$

∎

######  Remark 34.

Unlike the stiffness matrix case, this estimate for the mass matrix case is
sharp. Since all of the entries of the mass matrix are positive, therefore in
the overlapping, the entries of element mass matrices are always added up
without any cancellations or reductions.

######  Lemma 35.

There exists a constant $C$ that is independent of $p$ such that the minimum
eigenvalue of the mass matrix $M$ can be bounded below by

$\lambda_{\text{min}}(M)\geq\frac{C}{p^{4}16^{p}}.$

###### Proof.

To bound the minimum eigenvalue from below we use the left hand side
inequality of (31):

$\displaystyle\frac{\\{v_{i}\\}\cdot M\\{v_{i}\\}}{\|\\{v_{i}\\}\|^{2}}=$
$\displaystyle\frac{(v,v)}{\|\\{v_{i}\\}\|^{2}}\geq\frac{\displaystyle\frac{C}{p^{4}16^{p}}\|\\{v_{i}\\}\|^{2}}{\|\\{v_{i}\\}\|^{2}}=\frac{C}{p^{4}16^{p}}.$

Therefore, $\displaystyle\lambda_{\text{min}}(M)\geq\frac{C}{p^{4}16^{p}},$
where $C$ is a constant that is independent of $p$. ∎

The following lemma gives us the upper bound for the condition number of the
mass matrix.

######  Lemma 36.

The condition number of the mass matrix $M$ is bounded above by

$\kappa(M)\leq Cp^{4}16^{p},$

where $C$ is a constant that is independent of $p$.

###### Proof.

From Lemma 31 and Lemma 35,

$\frac{C}{p^{4}16^{p}}\leq\lambda_{\text{min}}\leq\lambda_{\text{max}}\leq C.$

Hence,

$\kappa(M)\leq Cp^{4}16^{p}.$

∎

######  Remark 37.

The above bound can be easily generalized for a $d$-dimensional problem:

(47) $\kappa(M)\leq p^{2d}4^{pd}.$

Following Remark 26 and using de Boor’s conjecture (16b), the upper bound for
the condition number of the mass matrix can be further improved and given by

(48) $\kappa(M)\leq\displaystyle 4^{pd}.$

######  Remark 38.

We have done all the analysis on the parametric domain $(0,1)^{2}$. To get the
results for the physical domain we can define an invertible NURBS geometrical
map from the parametric domain to the physical domain. With suitable
transformations we get the results for the physical domain. For details, see
[5].

##  5 Numerical results

In this section, we provide the numerical results for $h$-refined (in Section
5.1) and $p$-refined (in Section 5.2) stiffness and mass matrices. The
numerical discretizations are performed using the Matlab toolbox GeoPDEs [16,
17].

###  5.1 $h$-refinement

For $h$-refinement, the condition number of the stiffness matrix is shown in
Table 2. Numerical results are provided from $p=2$ to $p=5$. In the classical
finite element method, the condition number of the stiffness matrix is of
order $h^{-2}$ even for a coarse mesh-size. However, in isogeometric
discretizations, for higher $p$ on coarse mesh, the condition number is highly
influenced by the stability constant of B-Splines. The condition number of
B-Splines heavily depends on the polynomial degree (see Section 3) and scales
as $(p2^{p})^{d}$. The factor $(p2^{p})^{d}$ dominates the factor $h^{-2}$ for
coarse meshes. Nevertheless, the numerical results support the theoretical
findings asymptotically (for reasonably refined meshes) for any polynomial
degree.

Table 2:  Condition number of the stiffness matrix $A$ $p$ $h^{-1}$ | 2 | 4 | 8 | 16 | 32 | 64 | 128
---|---|---|---|---|---|---|---
2 | 4.00 | 4.00 | 5.22 | 19.77 | 78.14 | 311.58 | 1245.36
3 | 30.93 | 29.51 | 29.19 | 28.56 | 82.10 | 327.21 | 1307.67
4 | 339.92 | 269.23 | 240.03 | 222.55 | 215.00 | 381.73 | 1525.40
5 | 4177.20 | 3220.60 | 2148.25 | 1812.58 | 1700.63 | 1688.11 | 1781.51

In Table 3, we present the condition number of the mass matrix. We see that
the condition number is bounded uniformly by a constant independent of $h$,
which confirms the theoretical estimates.

Table 3:  Condition number of the mass matrix $M$ $p$ $h^{-1}$ | 2 | 4 | 8 | 16 | 32 | 64 | 128
---|---|---|---|---|---|---|---
2 | 89.679 | 109.68 | 108.51 | 109.85 | 111.29 | 111.69 | 111.79
3 | 915.558 | 799.941 | 737.379 | 708.010 | 715.89 | 719.45 | 720.33
4 | 11773.17 | 6795.46 | 5381.96 | 4762.53 | 4750.07 | 4779.41 | 4786.90
5 | 163371.70 | 77448.11 | 42580.04 | 33560.40 | 32587.27 | 32808.69 | 32871.70

###  5.2 $p$-refinement

We perform numerical experiments for $p$-refinement to obtain the maximum and
minimum eigenvalues, and the condition number of the stiffness matrix and the
mass matrix. The eigenvalues and the condition number are obtained on the
coarsest mesh and the finest mesh. For higher $p$ ($p>10$) roundoff errors
start contaminating the results and we stop reporting with $10$.

In Tables 4 and 5, we present the extremal eigenvalues and the condition
number of the stiffness matrix for $p=2$ to $p=10$. We observe that the
maximum eigenvalue scales as a constant independent of $p$ for the coarsest
mesh and linearly dependent on $p$ for refined meshes, and that the minimum
eigenvalue is bounded from below by the bound given in Theorem 16.

The extremal eigenvalues and the condition number of the mass matrix for $p=2$
to $p=10$ are presented in Table 6 and Table 7. Numerical results confirm the
theoretical estimates given in Lemma 31, Lemma 33, Lemma 35, and Lemma 36.

Table 4:  $\lambda_{\text{max}}$, $\lambda_{\text{min}}$, and $\kappa(A)$ on the coarsest mesh $p$ | $\lambda_{\text{max}}$ | $\lambda_{\text{min}}$ | $\kappa(A)$
---|---|---|---
2 | 0.35 | 3.5e-01 | 1.0e+00
3 | 0.45 | 3.8e-02 | 1.1e+01
4 | 0.41 | 2.9e-03 | 1.3e+02
5 | 0.35 | 2.1e-04 | 1.6e+03
6 | 0.33 | 1.5e-05 | 2.1e+04
7 | 0.33 | 1.1e-06 | 2.9e+05
8 | 0.31 | 7.8e-08 | 4.0e+06
9 | 0.30 | 5.4e-09 | 5.6e+07
10 | 0.30 | 3.7e-10 | 8.1e+08
Table 5:  $\lambda_{\text{max}}$, $\lambda_{\text{min}}$, and $\kappa(A)$ on the finest mesh $p$ | $\lambda_{\text{max}}$ | $\lambda_{\text{min}}$ | $\kappa(A)$
---|---|---|---
2 | 1.50 | 1.2e-03 | 1.2e+03
3 | 1.58 | 1.2e-03 | 1.3e+03
4 | 1.84 | 1.2-03 | 1.5e+03
5 | 2.14 | 1.2e-03 | 1.8e+03
6 | 2.47 | 1.8e-04 | 1.4e+04
7 | 2.80 | 2.6e-05 | 1.1e+05
8 | 3.14 | 3.6e-06 | 8.8e+05
9 | 3.47 | 4.9e-07 | 7.1e+06
10 | 3.81 | 6.6e-08 | 5.7e+07
Table 6:  $\lambda_{\text{max}}$, $\lambda_{\text{min}}$, and $\kappa(M)$ on the coarsest mesh $p$ | $\lambda_{\text{max}}$ | $\lambda_{\text{min}}$ | $\kappa(M)$
---|---|---|---
2 | 1.1e-01 | 1.1e-03 | 1.0e+02
3 | 6.2e-02 | 5.1e-05 | 1.2e+03
4 | 4.0e-02 | 2.5e-06 | 1.5e+04
5 | 2.7e-02 | 1.3e-07 | 2.1e+05
6 | 2.0e-02 | 6.9e-09 | 2.9e+06
7 | 1.5e-02 | 3.7e-10 | 4.1e+07
8 | 1.2e-02 | 2.0e-11 | 5.9e+08
9 | 1.0e-02 | 1.1e-12 | 8.5e+09
10 | 8.2e-03 | 6.6e-14 | 1.2e+11
Table 7:  $\lambda_{\text{max}}$, $\lambda_{\text{min}}$, and $\kappa(M)$ on the finest mesh $p$ | $\lambda_{\text{max}}$ | $\lambda_{\text{min}}$ | $\kappa(M)$
---|---|---|---
2 | 6.1e-05 | 5.5e-07 | 1.1e+02
3 | 6.1e-05 | 8.5e-08 | 7.2e+02
4 | 6.1e-05 | 1.3e-08 | 4.8e+03
5 | 6.1e-05 | 1.9e-09 | 3.3e+04
6 | 6.1e-05 | 2.6e-10 | 2.3e+05
7 | 6.1e-05 | 3.7e-11 | 1.7e+06
8 | 6.1e-05 | 5.1e-12 | 1.2e+07
9 | 6.1e-05 | 6.9e-13 | 8.9e+07
10 | 6.1e-05 | 9.2e-14 | 6.6e+08

##  6 Conclusions

We have provided the bounds for the minimum eigenvalue, maximum eigenvalue,
and the condition numbers of the stiffness and mass matrices for the Laplace
operator with $h$\- and $p$-refinements of the isogeometric discretizations
that are based on B-Spline (NURBS) basis functions. We proved that in the
$h$-refinement case, like the classical finite element method, the condition
number of the stiffness matrix scales as $h^{-2}$. For the mass matrix, it
scales as constant independent of $h$. For the $p$-refinement case, we proved
that the condition number of the stiffness and mass matrices grow
exponentially in $p$.

The estimates for the minimum eigenvalues of the stiffness and mass matrices
depend on the stability constant of B-Splines. In reaching these estimates we
have used the stability constant of B-Splines as $p2^{p}$. Using the de Boor’s
conjecture (the stability constant of B-Splines given by $2^{p}$, which is the
best known bound), these estimates can be further improved according to
Remarks 26 and 37.

Unfortunately, a sharp estimate for the stability constant is unknown.
Therefore, a sharp estimate for the minimum eigenvalue cannot be determined at
this time and will be the subject of future research by us and others. It is a
very difficult problem.

## Acknowledgments

The authors would like to thank Prof. U. Langer (Johannes Kepler University
Linz, Austria) and Prof. L. Zikatanov (Pennsylvania Stae University, USA) for
helpful suggestions on the topic of this paper.

## References

  * [1] F. Auricchio, L. Beirao da Veiga, A. Buffa, C. Lovadina, A. Reali and G. Sangalli. A fully “locking-free” isogeometric approach for plane linear elasticity problems: A stream function formulation. Comput. Methods Appl. Mech. Engrg. 197, 160-172, 2007.
  * [2] O. Axelsson and V. A. Barker. Finite Element Solution of Boundary Value Problems: Theory and Computation. Society for Industrial and Applied Mathematics Philadelphia, PA, USA, 2001.
  * [3] I. Babuska, J.E. Osborn. Finite element-Galerkin approximation of the eigenvalues and eigenvectors of selfadjoint problems Mathematics of computation, Volume 52, No. 186. 275-297, 1989.
  * [4] I. Babuska, B.A. Szabo, and I.N. Katz. The $p$-version of the Finite Element Method. SIAM Journal on Numerical Analysis, Volume 18, No. 3. 515-545, 1981.
  * [5] Y. Bazilevs, L. Beirao Da Veiga, J.A. Cottrell, T.J.R. Hughes and G. Sangalli. Isogeometric analysis: Approximation, Stability and error estimates for $h$\- refined meshes. Math. Models Methods Appl. Sci. 16(7), 1031-1090, 2006.
  * [6] R. Bellman. A note on an inequality of E. Schmidt. Amer. Math. Soc., 50, 734-736, 1944.
  * [7] D. Braess. Finite Elements: Theory, Fast Solvers and Applications in Solid Mechanics. Cambridge University Press, 2007.
  * [8] P.G. Ciarlet. The Finite Element Method for Elliptic Problems. North Holland Publishing Company, 1978.
  * [9] J.A. Cottrell, T.J.R. Hughes and Y. Bazilevs. Isogeometric Analysis: Toward Integration of CAD and FEA. Wiley, 2009.
  * [10] C. de Boor. On Calculating with B-Splines. Journal of Approximation Theory 6, pp. 50-62, 1972.
  * [11] C. de Boor. On local linear functionals which vanish at all $B$-splines but one. Theory of Approximation with Applications, A. G. Law and N. B. Sahney (eds.), Academic Press (New York). 120-145, 1976.
  * [12] C. de Boor. Splines as linear combinations of B-Splines, a survey. In Approximation Theory II, G. G. Lorentz, C. K. Chui and L. L. Schumaker, (Eds.), Academic Press (New York), pp. 1-47, 1976.
  * [13] C. de Boor. A Practical Guide to Splines. Springer-Verlag, New York, 1978.
  * [14] C. de Boor. The exact condition of the B-Spline basis may be hard to determine. Journal of Approximation theory. 60, 344-359, 1990.
  * [15] C. de Boor and J.W. Daniel. Splines with Nonnegative B-Spline Coefficients. Mathematics of Computation 28(126), pp. 565-568, 1974.
  * [16] C. de Falco, A. Reali and R. Vazquez. GeoPDEs: A research tool for Isogeometric Analysis of PDEs. Adv. Eng. Softw. 42, 1020-1034, 2011.
  * [17] C. de Falco, A. Reali and R. Vazquez. GeoPDEs
webpage. http://geopdes.sourceforge.net

  * [18] A. Ern, J.-L. Guermond. Evaluation of the condition number in linear systems arising in finite element approximations. ESAIM Math. Mod. Numer. Anal. 40(1), 29-48, 2006.
  * [19] M.S. Floater. Evaluation and properties of the derivative of a NURBS curve. Mathematical Methods in CAGD, T. Lyche and L.L. Schumaker (eds.), Academic Press, Boston, 2, pp. 261-274, 1992.
  * [20] K.P.S. Gahalaut, J.K. Kraus and S.K. Tomar. Multigrid Methods for Isogeometric Discretization. Comput. Methods Appl. Mech. Engrg., 253, pp. 413–425, 2013.
  * [21] K.P.S. Gahalaut and S.K. Tomar. Condition number estimates for matrices arising in the isogeometric discretizations. RICAM report, 23–2012.
  * [22] C. Garoni, C. Manni, F. Pelosi, S. Serra-Capizzano and H. Speleers. On the spectrum of stiffness matrices arising from isogeometric analysis. Numerische Mathematik. DOI 10.1007/s00211-013-0600-2, 2013.
  * [23] W. Hackbusch. Iterative Solution of Large Sparse Systems of Equations. Springer, 1994.
  * [24] K. Hoellig. Multivariate Splines. SIAM J. Numer. Anal. 19(5), pp. 1013-1031, 1982.
  * [25] K. Hoellig. Stability of B-Spline basis via knot insertion. Computer Added Geometric Design. 17, pp. 447-450, 2000.
  * [26] K. Hoellig, U. Reif and J. Wipper. Weighted Extended B-Spline Approximation of Dirichlet Problems. SIAM J. Numer. Anal. 39(2), pp. 442-462, 2002.
  * [27] N. Hu, X-Z Guo, and I.N. Katz. Bounds for eigenvalues and condition numbers in the $p$-version of the finite element method. Mathematics of computation. 67/224, 1423-1450, 1998.
  * [28] T.J.R. Hughes, J.A. Cottrell and Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput. Methods Appl. Mech. Engrg. 194, 4135-4195, 2005.
  * [29] C. Johnson. Numerical solution of partial differential equations by the finite element method. Cambridge University Press, Cambridge, New York, 1987.
  * [30] T. Lyche, K. Scherer. On the sup-norm condition number of the multivariate triangular Bernstein basis. Multivariate Approximation and Splines, G. Nuernberger, J. W. Schmidt, and G. Walz (eds.), ISNM.125, Birkhäuser Verlag, Basel. 141-151, 1997.
  * [31] T. Lyche, K. Scherer. On the $p$-norm condition number of the multivariate triangular Bernstein basis Journal of Computational and Applied Mathematics. 119, 259–273, 2000.
  * [32] J.F. Maitre and O. Pourquier. About the conditioning of matrices in the $p$-version of the finite element method for second order elliptic problems. Journal of Computational and Applied Mathematics 63, 341-348, 1995.
  * [33] J.F. Maitre and O. Pourquier. Condition number and diagonal preconditioning: Comparison of the $p$-version and the spectral element methods. Numer. Math. 74, pp. 69-84, 1996.
  * [34] J.M. Melenk. On condition numbers in hp-FEM with Gauss–Lobatto-based shape functions. Journal of Computational and Applied Mathematics. Volume 139, Issue 1. 21-48, 2002.
  * [35] J.M. Melenk and I. Babuska. The Partition of Unity Finite Element Method: Basic Theory and Applications. Comput. Methods Appl. Mech. Engrg.139, pp. 289-314, 1996.
  * [36] B. Mössner and U. Reif. Stability of tensor product B-Splines on domains. J. Approx. Theory 154, 1-19, 2008.
  * [37] E.T. Olsen, J. Douglas, Jr. Bounds on spectral condition numbers of matrices arising in the $p$-version of the finite element method. Numerische Mathematik. 69, 333-352, 1995.
  * [38] J.M. Pe$\tilde{\text{n}}$a. B-Splines and Optimal Stability. Mathematics of Computation 66, pp. 1555-1560, 1997.
  * [39] L. Piegl and W. Tiller. The NURBS Book (Monographs in Visual Communication), Second ed., Springer-Verlag, 1997.
  * [40] E. Pilgerstorfer, B. Jüttler. Bounding the influence of domain parameterization and knot spacing on numerical stability in Isogeometric Analysis. Computer Methods in Applied Mechanics and Engineering. 268, pp. 589–613, 2014.
  * [41] P.M. Prenter. Splines and Variational Methods. John Wiley & Sons, New York, 1975.
  * [42] D.F. Rogers. An Introduction to NURBS With Historical Perspective. Academic Press, 2001.
  * [43] Y. Saad. Iterative Methods for Sparse Linear Systems. PWS Publishing Company, Boston, 1996.
  * [44] K. Scherer, A.Y. Shadrin. New upper bound for the B-Spline basis condition number, I. East J. Approx. 2, 331-342, 1996.
  * [45] K. Scherer, A.Y. Shadrin. New upper bound for the B-Spline basis condition number, II. A proof of de Boor’s $2^{k}$-conjecture. Journal of Approximation theory. 99, 217-229, 1999.
  * [46] M.H. Schultz. Spline Analysis. Prentice-Hall, Englewood Cliffs, 1973.
  * [47] L.L. Schumaker. Spline Functions: Basic Theory. Cambridge University Press, 2007.
  * [48] G. Strang and G.J. Fix. An Analysis of the Finite Element Method. Prentice Hall, 1973.
  * [49] R.S. Varga. Matrix Iterative Analysis. Prentice Hall, New Jersey, 1965.
  * [50] A. Weiser, S.C. Eisenstat, and M.H. Schultz. On solving elliptic equations to moderate accuracy. SIAM J. Numer. Anal. 17, 908-929, 1980.