mydate[],

# Multigrid solvers for isogeometric discretizations of the second biharmonic
problem

Jarle Sogn1 and Stefan Takacs2
1 Johann Radon Institute for Computational and Applied Mathematics,
Austrian Academy of Sciences,
Altenberger Str. 69, 4040 Linz, Austria
2 Institute of Computational Mathematics,
Johannes Kepler University Linz,
Altenberger Str. 69, 4040 Linz, Austria

###### Abstract

We develop a multigrid solver for the second biharmonic problem in the context
of Isogeometric Analysis (IgA), where we also allow a zero-order term. In a
previous paper, the authors have developed an analysis for the first
biharmonic problem based on Hackbusch’s framework. This analysis can only be
extended to the second biharmonic problem if one assumes uniform grids. In
this paper, we prove a multigrid convergence estimate using Bramble’s
framework for multigrid analysis without regularity assumptions. We show that
the bound for the convergence rate is independent of the scaling of the zero-
order term and the spline degree. It only depends linearly on the number of
levels, thus logarithmically on the grid size. Numerical experiments are
provided which illustrate the convergence theory and the efficiency of the
proposed multigrid approaches.

##  1 Introduction

We consider multigrid methods for biharmonic problems discretized by
Isogeometric Analysis (IgA). In particular, we consider the following model
problem: Given a bounded domain $\Omega\subset\mathbb{R}^{d}$,
$d\in\\{2,3\\}$, with Lipschitz boundary $\partial\Omega$, a parameter
$\beta\geq 0$ and sufficiently smooth functions $f$, $g_{1}$, and $g_{2}$,
find a function $u$ such that

$\displaystyle\begin{split}\beta u+\Delta^{2}u&=f\quad\text{in}\quad\Omega,\\\
u&=g_{1}\quad\text{on}\quad\partial\Omega,\\\ \Delta
u&=g_{2}\quad\text{on}\quad\partial\Omega\end{split}$ (1.1)

holds in a variational sense. For $\beta=0$, this problem is known as the
second biharmonic problem, which is of interest for plate theory (cf. [6]) and
Stokes streamline equations (cf. [9]). Problems with $\beta>0$ are of
particular interest in the context of optimal control problems, where the
constraint is a second order elliptic operator. The optimality systems
associated to these optimal control problems can be preconditioned robustly
using preconditioners that rely on solving (1.1), see [18, 26, 1, 19]. The
problem (1.1) is obtained when considering the full observation; if one
considers an optimal control problem with limited observation, one would
obtain a similar problem, where the mass term $\beta u$ is multiplied with the
characteristic function for the observation domain.

We derive a standard variational formulation of the model problem, which lives
in the Sobolev space $H^{2}(\Omega)$. For the discretization, we use
Isogeometric Analysis (IgA) since it easily allows for $H^{2}$-conforming
discretizations. Particularly, we consider a discretization based on tensor
product B-splines of some degree $p>1$ and maximum smoothness, i.e., $p-1$
times continuously differentiable. For the derivation of the multigrid solver,
we set up a hierarchy of grids as obtained by uniform refinement. Since we
keep spline degree and spline smoothness fixed, we obtain nested spaces.

Concerning the choice of the smoother, there are many possibilities. We are
interested in a smoother that yields a $p$-robust multigrid method. The first
$p$-robust multigrid solvers were based on the boundary corrected mass
smoother [15] and the subspace corrected mass smoother [14]. Both have been
formulated for the Poisson problem. Since the subspace corrected mass smoother
is more flexible and has proven itself more efficient in practice, we restrict
ourselves to that smoother. The multigrid solvers with subspace corrected mass
smoother have been extended to the first biharmonic problem in [25] and to the
second and third biharmonic problem in the thesis [24]. The convergence
estimates are shown using the standard splitting of the analysis into
approximation property and smoothing property, as proposed by Hackbusch, cf.
[12].

The theory in all of these papers requires that the grids are uniform since
they have been based on the $p$-robust approximation error estimates from
[28], which are valid only in this case. Since then, newer $p$-robust
approximation error estimates, see [23, 22], have been proposed, which do not
require uniform grids. Using these new estimates, it is straightforward to
relax this assumption and to show analogous results for the Poisson problem as
well as the first biharmonic problem for quasi uniform grids. However, this is
not straightforward for the second biharmonic problem, since the proof
requires a certain commutativity property (cf. [24, Lemma 9.2]), which is only
valid in case of uniform grids.

In this paper, we go another way. We base the analysis on the framework
introduced by Bramble et al., cf. [3, 2]. This allows us to drop the
requirement that the grids are uniform. While this analysis could also be
performed for other kinds of boundary conditions, like the first biharmonic
problem, we restrict ourselves to the second biharmonic problem since it has
previously turned out to be the more challenging one. For this setting, we
prove a multigrid convergence estimate which is robust with respect to the
spline degree $p$ and which only depends logarithmically on the grid size $h$.

Moreover, we show that the convergence is robust in the parameter $\beta\geq
0$. This analysis is motivated by the mentioned optimal control problem. Such
parameter-robust multigrid solvers are also known for the Poisson problem, see
[20] for an analysis based on Hackbusch’s framework. There, the authors also
provide a regularity result for the corresponding partial differential
equation (PDE), which is based on standard results for the Poisson problem. In
our case, we do not need to do that since Bramble’s analysis is not based on
any regularity assumptions.

In the numerical experiments, one can observe that the convergence of a
multigrid solver with subspace corrected mass smoother degrades if the
geometry gets distorted. While this is also true for the Poisson problem, this
dependence is significantly amplified for the biharmonic problem. The reason
for the geometry dependence of the convergence rates is that the subspace
corrected mass smoother is based on the tensor product structure of the spline
space. This tensor product structure is distorted by the geometry mapping. So,
the contributions of the geometry function are ignored when setting up the
smoother. We aim to overcome this problem by considering a hybrid smoother
that combines the proposed smoother with Gauss-Seidel sweeps, see also [25,
24].

Alternative smoothers based on overlapping multiplicative Schwarz techniques
have been considered in [7, 19]. Both approaches give good numerical results
for the biharmonic problem. However, there is no rigorous, $p$-robust
convergence theory available for these methods. It is worth mentioning that,
as an alternative for solving biharmonic problems on the primal form, various
kinds of mixed or non-conforming formulations have been developed, cf. [4, 29,
13, 21, 5].

The remainder of the paper is organized as follows. We introduce IgA, the
biharmonic model problem in its variational form and its discretization in
Section 2. In Section 3, the multigrid method is introduced and we state
sufficient conditions for its convergence. We develop the approximation error
estimates needed for the convergence estimates in Section 4. The choice of the
smoother, the smoothing properties and the resulting multigrid convergence
results are addressed in Section 5. Finally, we provide numerical results in
Section 6.

##  2 Model problem and its discretization

###  2.1 The biharmonic model problem

Following the usual design principles of IgA, we assume that the computational
domain $\Omega\subset\mathbb{R}^{d}$ has a Lipschitz boundary $\partial\Omega$
and that it is parameterized by a geometry function

$\textbf{G}:\widehat{\Omega}=(0,1)^{d}\rightarrow\Omega=\textbf{G}(\widehat{\Omega}),$

whose third weak derivatives are almost everywhere uniformly bounded. The
parameterization has the property

$\|\nabla^{r}\mathbf{G}\|_{L^{\infty}(\widehat{\Omega})}\leq
c_{1}\quad\text{and}\quad\|\left(\nabla^{r}\mathbf{G}\right)^{-1}\|_{L^{\infty}(\widehat{\Omega})}\leq
c_{2},\quad\text{for}\quad r=1,2,3,$ (2.1)

for some constants $c_{1}$ and $c_{2}$.

After homogenization, the variational formulation of the model problem (1.1)
reads as follows. Given $f\in L^{2}(\Omega)$ and $\beta\in\mathbb{R}$ with
$\beta\geq 0$, find $u\in V:=H^{2}(\Omega)\cap H^{1}_{0}(\Omega)$ such that

$\beta(u,v)_{L^{2}(\Omega)}+(\Delta u,\Delta
v)_{L^{2}(\Omega)}=(f,v)_{L^{2}(\Omega)}\quad\forall\,v\in V.$ (2.2)

Here and in what follows, $L^{2}(\Omega)$ and $H^{r}(\Omega)$ denote the
standard Lebesgue and Sobolev spaces with standard inner products
$(\cdot,\cdot)_{L^{2}(\Omega)}$, $(\cdot,\cdot)_{H^{r}(\Omega)}$ and norms
$\|\cdot\|_{L^{2}(\Omega)}$, $\|\cdot\|_{H^{r}(\Omega)}$. $H^{1}_{0}(\Omega)$
is the standard subspace of $H^{1}(\Omega)$ containing the functions with
vanishing trace. On $V$, we define the bilinear form
$\left(\cdot,\cdot\right)_{\mathcal{B}}$ via

$\left(u,v\right)_{\mathcal{B}}:=(\Delta u,\Delta
v)_{L^{2}(\Omega)}\quad\forall\,u,v\in V,$

which is an inner product since we have the Poincaré like inequality

$\|u\|_{H^{2}(\Omega)}\leq c_{\Omega}\|\Delta
u\|_{L^{2}(\Omega)}=c_{\Omega}\|u\|_{\mathcal{B}}\quad\forall\,u\in V,$ (2.3)

where $c_{\Omega}$ is a constant that depends only on the shape of $\Omega$,
cf. [19].

Using the substitution rule for integration and the chain rule for
differentiation, (2.2) can be expressed in terms of integrals on the parameter
domain $\widehat{\Omega}$. In IgA, this is usually done in order to simplify
the evaluation of the integrals using quadrature rules. Besides these inner
products, there are also standard inner products for the parameter domain,
like $(\cdot,\cdot)_{L^{2}(\widehat{\Omega})}$ and
$(\cdot,\cdot)_{\widehat{\mathcal{B}}}$, where the latter is given by

$\left(\widehat{u},\widehat{v}\right)_{\widehat{\mathcal{B}}}:=(\Delta\widehat{u},\Delta\widehat{v})_{L^{2}(\widehat{\Omega})}\quad\forall\,\widehat{u},\widehat{v}\in\widehat{V}:=H^{2}(\widehat{\Omega})\cap
H^{1}_{0}(\widehat{\Omega}).$

Also for the parameter domain $\widehat{\Omega}$, the result (2.3) holds. So,
we know

$\|u\|_{H^{2}(\widehat{\Omega})}\leq c_{\widehat{\Omega}}\|\Delta
u\|_{L^{2}(\widehat{\Omega})}=c_{\widehat{\Omega}}\|u\|_{\widehat{\mathcal{B}}}\quad\forall\,u\in\widehat{V}.$

We know (cf. [24]) that there exist constants $\underline{c}_{M}$,
$\overline{c}_{M}$, $\underline{c}_{B}$ and $\overline{c}_{B}$ only depending
on the constants $c_{1}$, $c_{2}$ and the shape of $\Omega$ such that

$\begin{split}\underline{c}_{M}\,(u,u)_{L^{2}(\Omega)}&\leq(\widehat{u},\widehat{u})_{L^{2}(\widehat{\Omega})}\leq\overline{c}_{M}\,(u,u)_{L^{2}(\Omega)}\quad\text{and}\\\
\underline{c}_{B}\,(u,u)_{\mathcal{B}}&\leq(\widehat{u},\widehat{u})_{\widehat{\mathcal{B}}}\leq\overline{c}_{B}\,(u,u)_{\mathcal{B}}\end{split}$
(2.4)

for all $u\in V$ with $\widehat{u}=u\circ\bm{G}\in\widehat{V}$. We define a
simplified bilinear form $\left(\cdot,\cdot\right)_{\bar{\mathcal{B}}}$ as the
inner product obtained by removing the cross terms from the inner product
$\left(\cdot,\cdot\right)_{\widehat{\mathcal{B}}}$, that is,

$\left(\widehat{u},\widehat{v}\right)_{\bar{\mathcal{B}}}:=\sum^{d}_{k=1}\left(\partial_{x_{k}x_{k}}\widehat{u},\partial_{x_{k}x_{k}}\widehat{v}\right)_{L^{2}(\widehat{\Omega})}\quad\forall\,\widehat{u},\widehat{v}\in\widehat{V}.$

Here and in what follows, $\partial_{x}:=\frac{\partial}{\partial x}$ and
$\partial_{xy}:=\partial_{x}\partial_{y}$ and
$\partial_{x}^{r}:=\frac{\partial^{r}}{\partial x^{r}}$ denote partial
derivatives. The original bilinear form and the simplified bilinear form are
spectrally equivalent, which implies that also the simplified bilinear form is
an inner product.

######  Lemma 2.1.

The inner products $(\cdot,\cdot)_{\widehat{\mathcal{B}}}$ and
$(\cdot,\cdot)_{\bar{\mathcal{B}}}$ are spectrally equivalent, that is,

$\left(\widehat{u},\widehat{u}\right)_{\bar{\mathcal{B}}}\leq\left(\widehat{u},\widehat{u}\right)_{\widehat{\mathcal{B}}}\leq
d\left(\widehat{u},\widehat{u}\right)_{\bar{\mathcal{B}}}\quad\forall\,\widehat{u}\in\widehat{V}.$

###### Proof.

From [11, 10], it follows that
$\|\Delta\widehat{u}\|_{L^{2}(\widehat{\Omega})}=\|\nabla^{2}\widehat{v}\|_{L^{2}(\widehat{\Omega})}$
for $\widehat{u},\widehat{v}\in\widehat{V}$. Using this, we obtain

$\displaystyle\|\widehat{u}\|_{\widehat{\mathcal{B}}}^{2}$
$\displaystyle=\|\Delta\widehat{u}\|_{L^{2}(\widehat{\Omega})}^{2}=\|\nabla^{2}\widehat{u}\|_{L^{2}(\widehat{\Omega})}^{2}=\underbrace{\sum_{k=1}^{d}\|\partial_{x_{k}x_{k}}\widehat{u}\|_{L^{2}(\widehat{\Omega})}^{2}}_{\displaystyle=\|\widehat{u}\|_{\bar{\mathcal{B}}}^{2}}+\underbrace{\sum_{k=1}^{d}\sum_{l\in\\{1,\ldots,d\\}\backslash\\{k\\}}\|\partial_{x_{k}x_{l}}\widehat{u}\|_{L^{2}(\widehat{\Omega})}^{2}}_{\displaystyle\geq
0},$

which shows the first side of the inequality. Using the Cauchy-Schwarz
inequality and $ab\leq\tfrac{1}{2}(a^{2}+b^{2})$, we obtain

$\displaystyle\|\widehat{u}\|_{\widehat{\mathcal{B}}}^{2}$
$\displaystyle=\sum_{k=1}^{d}\sum_{l=1}^{d}\left(\partial^{2}_{x_{k}}\widehat{u},\partial^{2}_{x_{l}}\widehat{u}\right)_{L^{2}(\widehat{\Omega})}\leq\frac{1}{2}\sum_{k=1}^{d}\sum_{l=1}^{d}\left(\left\|\partial_{x_{k}}^{2}\widehat{u}\right\|_{L^{2}(\widehat{\Omega})}^{2}+\left\|\partial_{x_{l}}^{2}\widehat{u}\right\|_{L^{2}(\widehat{\Omega})}^{2}\right)=d\|\widehat{u}\|_{\bar{\mathcal{B}}}^{2},$

which shows second side of the inequality. ∎

######  Remark 2.1.

A analogous result holds for the domain $\Omega$, which satisfies condition
(2.1). In this case, the constants also depend on the shape of $\Omega$.

###  2.2 Discretization

We consider a discretization using tensor product B-splines in the context of
IgA. We start by defining these splines on the parameter domain
$\widehat{\Omega}$. Let $C^{k}(0,1)$ denote the space of all continuous
functions mapping $(0,1)\rightarrow\mathbb{R}$ that are $k$ times continuously
differentiable and let $\mathcal{P}_{p}$ be the space of polynomials of degree
at most $p$. For any sequence of grid points
$\bm{\tau}:=(\tau_{0},\ldots,\tau_{N+1})$ with

$0=\tau_{0}<\tau_{1}<\cdots<\tau_{N}<\tau_{N+1}=1,$

we define the space $S_{p,\bm{\tau}}$ of splines of degree $p$ with maximum
smoothness by

$S_{p,\bm{\tau}}:=\left\\{v\in
C^{p-1}(0,1):v|_{(\tau_{j},\tau_{j+1})}\in\mathcal{P}_{p},\;j=0,1,\ldots,N\right\\}.$

The size of the largest and the smallest interval are denoted by

$h_{\bm{\tau}}:=\max_{j=0,\ldots,N}(\tau_{j+1}-\tau_{j})\quad\text{and}\quad
h_{\bm{\tau},\mathrm{min}}:=\min_{j=0,\ldots,N}(\tau_{j+1}-\tau_{j}),$

respectively. For the parameter domain, we define a spline space by
tensorization, which we transfer to the physical domain using the pull-back
principle, thus we define for given sequences of grid points
$\bm{\tau}_{\ell,1},\ldots,\bm{\tau}_{\ell,d}$ the spaces

$\widehat{V}_{\ell}:=\left(\bigotimes^{d}_{i=1}S_{p,\bm{\tau}_{\ell,i}}\right)\cap
H^{1}_{0}(\widehat{\Omega})\subset\widehat{V}\quad\text{and}\quad
V_{\ell}:=\\{f\circ\textbf{G}^{-1}:f\in\widehat{V}_{\ell}\\}\subset V.$

Here and in what follows, the tensor product space
$\bigotimes^{d}_{i=1}S_{p,\bm{\tau}_{\ell,i}}$ is the space of all linear
combinations of functions of the form
$v(x_{1},\ldots,x_{d})=v_{1}(x_{1})\cdots v_{d}(x_{d})$ with $v_{i}\in
S_{p,\bm{\tau}_{\ell,i}}$. The spline degree $p$ could be different for each
of the spacial directions. For notational convenience, we restrict ourselves
to a uniform choice of the degree.

The corresponding minimum and maximum grid size are denoted by

$h_{\ell}:=\max_{i=1,\ldots,d}h_{\bm{\tau}_{\ell,i}}\quad\text{and}\quad
h_{\ell,\mathrm{min}}:=\min_{i=1,\ldots,d}h_{\bm{\tau}_{\ell,i},\mathrm{min}}.$

For the multigrid methods we, set up a sequence nested spline spaces

$V_{0}\subset V_{1}\subset\cdots\subset V_{L}\subset V\quad\text{with}\quad
h_{0}>h_{1}>\cdots>h_{L}>0$

based on a sequence of nested grids.

We assume that all grids are quasi uniform, that is, there is a constant
$c_{q}$ such that

$h_{\ell}\leq
c_{q}\,h_{\ell,\mathrm{min}}\quad\text{for}\quad\ell=0,1,\ldots,L.$ (2.5)

We also assume that the ratio of the grid sizes of any two consecutive grids
is bounded, that is, there is a constant $c_{r}$ such that

$h_{\ell}\leq c_{r}\,h_{\ell-1}\quad\text{for}\quad\ell=1,\ldots,L.$ (2.6)

If the grids are obtained by uniform refinements of the coarsest grid, then
this condition is naturally satisfied with $c_{r}=2$.

By applying a Galerkin discretization, we obtain the following discrete
problem: Find $u_{\ell}\in V_{\ell}$ such that

$\beta(u_{\ell},v_{\ell})_{L^{2}(\Omega)}+(u_{\ell},v_{\ell})_{\mathcal{B}}=(f,v_{\ell})_{L^{2}(\Omega)}\quad\forall\,v_{\ell}\in
V_{\ell}.$ (2.7)

By fixing a basis for the space $V_{\ell}$, we can rewrite (2.7) in matrix-
vector notation as

$(\beta\mathcal{M}_{\ell}+\mathcal{B}_{\ell})\underline{u}_{\ell}=\underline{f}_{\ell},$
(2.8)

where $\mathcal{B}_{\ell}$ is the biharmonic stiffness matrix,
$\mathcal{M}_{\ell}$ is the mass matrix, $\underline{u}_{\ell}$ is the vector
representation of the corresponding function $u_{\ell}$ with respect to the
chosen basis and the vector $\underline{f}_{\ell}$ is obtained by testing the
right-hand side functional $(f,\cdot)_{L^{2}(\Omega)}$ with the basis
functions.

######  Notation 2.1.

Throughout this paper, $c$ is a generic positive constant that is independent
of $h$ and $p$, but may depend on $d$, the constants $c_{1}$, $c_{2}$,
$c_{q}$, and $c_{r}$ and the shape of $\Omega$.

For any two square matrices $A,B\in\mathbb{R}^{n\times n}$, $A\leq B$ means
that

$\underline{x}^{T}A\underline{x}\leq\underline{x}^{T}B\underline{x}\quad\forall\,\underline{x}\in\mathbb{R}^{n}.$

##  3 The multigrid solver

In this section, we present an abstract multigrid method and give a
convergence theorem that is based on the analysis by Bramble et al., see [3,
Theorem 1].

###  3.1 The multigrid framework

Let us assume that we have nested spaces $V_{0}\subset
V_{1}\subset\cdots\subset V_{L}\subset V$. Let $I^{\ell}_{\ell-1}$ be the
matrix representation of the canonical embedding from $V_{\ell-1}$ into
$V_{\ell}$ and let the restriction matrix $I^{\ell-1}_{\ell}$ be its
transpose, this is $I^{\ell-1}_{\ell}:=(I^{\ell}_{\ell-1})^{T}$.

On each grid level, $\ell=0,\ldots,L$, we have a linear system

$\mathcal{A}_{\ell}\,\underline{u}_{\ell}=\underline{f}_{\ell},$

which is obtained by discretizing a symmetric, bounded and coercive bilinear
form $a(\cdot,\cdot)$ in the space $V_{\ell}$ using the Galerkin principle.
The matrix induces a norm via
$\|\underline{u}_{\ell}\|_{\mathcal{A}_{\ell}}:=(\mathcal{A}_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})^{1/2}=\|\mathcal{A}_{\ell}^{1/2}\underline{u}_{\ell}\|$.
Here and in what follows, $(\cdot,\cdot)$ and $\|\cdot\|$ are the Euclidean
scalar product and norm, respectively. In the continuous setting, the matrix
can be represented by an operator

$\mathcal{A}:V\rightarrow
V^{\prime}\quad\text{with}\quad\mathcal{A}u=a(u,\cdot).$

We have
$\|u_{\ell}\|_{\mathcal{A}}=\|\underline{u}_{\ell}\|_{\mathcal{A}_{\ell}}$ for
all functions $u_{\ell}\in V_{\ell}$ with coefficient representation
$\underline{u}_{\ell}$.

For the analysis, we can additionally choose symmetric positive definite
matrices $X_{\ell}$ for all grid levels $\ell=0,1,\ldots,L$, which induce
norms via
$\|\underline{u}_{\ell}\|_{X_{\ell}}=(X_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})^{1/2}=\|X^{1/2}_{\ell}\underline{u}_{\ell}\|$.
The norm $\|u_{\ell}\|_{X_{\ell}}$ of a function $u_{\ell}\in V_{\ell}$ is
interpreted as $\|\underline{u}_{\ell}\|_{X_{\ell}}$, where
$\underline{u}_{\ell}$ is the coefficient representation of $u_{\ell}$.

For the abstract framework, we assume to have a symmetric and positive
definite matrix $\tau_{\ell}L^{-1}_{\ell}$ for every grid level
$\ell=1,\ldots,L$, representing the smoother.

Later, for the model problem, the bilinear form $a(\cdot,\cdot)$, the matrices
$\mathcal{A}_{\ell}$, $\ell=0,\ldots,L$ and our choice of ${X}_{\ell}$ will be

$a(u,v)=\beta(u,v)_{L^{2}(\Omega)}+(u,v)_{\mathcal{B}},\quad\mathcal{A}_{\ell}=\beta\mathcal{M}_{\ell}+\mathcal{B}_{\ell}\quad\text{and}\quad
X_{\ell}=(\beta+h_{\ell}^{-4})\mathcal{M}_{\ell}+\mathcal{B}_{\ell}.$

As smoothers, we will choose a subspace corrected mass smoother, a symmetric
Gauss-Seidel smoother and a hybrid smoother in Section 5.

Based on these choices, the overall algorithm reads as follows.

######  Algorithm 3.1.

One multigrid cycle, applied to some iterate $\underline{u}_{\ell}^{(0)}$ and
a right-hand side $\underline{f}_{\ell}$ consists of the following steps:

  * •

Apply $\nu_{\ell}$ pre-smoothing steps, i.e., compute

$\underline{u}_{\ell}^{(i)}=\underline{u}_{\ell}^{(i-1)}+\tau_{\ell}L^{-1}_{\ell}(\underline{f}_{\ell}-\mathcal{A}_{\ell}\underline{u}_{\ell}^{(i-1)})\quad\mbox{for}\quad
i=1,\ldots,\nu_{\ell}.$ (3.1)

  * •

Apply recursive coarse-grid correction, i.e., apply the following steps.
Compute the residual and restrict it to the next coarser grid level:

$\underline{r}_{\ell-1}=I_{\ell}^{\ell-1}(\underline{f}_{\ell}-\mathcal{A}_{\ell}\underline{u}_{\ell}^{(\nu_{\ell})}).$

If $\ell-1=0$, compute the update
$\underline{q}_{0}:=A_{0}^{-1}\underline{r}_{0}$ using a direct solver.
Otherwise, compute the update $\underline{q}_{\ell-1}$ by applying the
algorithm $r$ ($r\in\mathbb{N}:=\\{1,2,\ldots\\}$) times recursively to the
right-hand side $\underline{r}_{\ell-1}$ and a zero vector as initial guess.
Then set

$\underline{u}^{(\nu_{\ell}+1)}_{\ell}=\underline{u}^{(\nu_{\ell})}_{\ell}+I_{\ell-1}^{\ell}\underline{q}_{\ell-1}.$

  * •

Apply $\nu_{\ell}$ post-smoothing steps, i.e., compute
$\underline{u}_{\ell}^{(i)}$ using (3.1) for
$i=\nu_{\ell}+2,\ldots,2\nu_{\ell}+1$ to obtain the next iterate
$\underline{u}_{\ell}^{(2\nu_{\ell}+1)}$.

This abstract algorithm coincides with the algorithm presented in [3]. Since
each multigrid cycle is linear, its application can be expressed by the matrix
$B_{\ell}^{s}$, which is recursively given by
$B_{0}^{s}:=\mathcal{A}_{0}^{-1}$ and

$B_{\ell}^{s}:=\big{(}I-(I-\tau_{\ell}L_{\ell}^{-1}\mathcal{A}_{\ell})^{\nu_{\ell}}(I-I_{\ell-1}^{\ell}B_{\ell-1}^{s}I_{\ell}^{\ell-1}\mathcal{A}_{\ell})^{r}(I-\tau_{\ell}L_{\ell}^{-1}\mathcal{A}_{\ell})^{\nu_{\ell}}\big{)}\mathcal{A}_{\ell}^{-1},\quad\ell=1,\ldots,L.$

The iteration matrix corresponding to one multigrid cycle is given by

$I-B_{\ell}^{s}\mathcal{A}_{\ell}=(I-\tau_{\ell}L_{\ell}^{-1}\mathcal{A}_{\ell})^{\nu_{\ell}}(I-I_{\ell-1}^{\ell}B_{\ell-1}^{s}I_{\ell}^{\ell-1}\mathcal{A}_{\ell})^{r}(I-\tau_{\ell}L_{\ell}^{-1}\mathcal{A}_{\ell})^{\nu_{\ell}},\quad\ell=1,\ldots,L.$

######  Remark 3.1.

The integer $r$ represents the recursively of the algorithm, where $r=1$
corresponds to the $V$-cycle and $r=2$ corresponds to the $W$-cycle.

###  3.2 Abstract convergence framework

The assumptions used to show convergence can be split into two groups:
approximation properties and smoother properties.

######  Theorem 3.1.

Let $\lambda_{\ell}$ be the largest eigenvalue of
$X^{-1}_{\ell}\mathcal{A}_{\ell}$. Assume that the following estimates hold:

  * •

_Approximation properties._ There are constants $C_{1}$ and $C_{2}$,
independent of $\ell$, and linear operators $Q_{\ell}:V_{L}\rightarrow
V_{\ell}$ for $\ell=0,1,\ldots,L$ with $Q_{L}=I$ such that

$\displaystyle\|(Q_{\ell}-Q_{\ell-1})u_{L}\|^{2}_{X_{\ell}}$
$\displaystyle\leq C_{1}\lambda_{\ell}^{-1}(u_{L},u_{L})_{\mathcal{A}}\quad$
for $\displaystyle\ell=1,\ldots,L,$ (3.2)
$\displaystyle(Q_{\ell}u_{L},Q_{\ell}u_{L})_{\mathcal{A}}$ $\displaystyle\leq
C_{2}(u_{L},u_{L})_{\mathcal{A}}\quad$ for $\displaystyle\ell=0,\ldots,L-1,$
(3.3)

for all $u_{L}\in V_{L}$.

  * •

_Smother properties._ We assume there exist a constant $C_{S}$ independent of
$\ell$ such that

$\frac{\|\underline{u}_{\ell}\|^{2}_{X_{\ell}}}{\lambda_{\ell}}\leq
C_{S}(\tau_{\ell}L^{-1}_{\ell}X_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})_{X_{\ell}}\quad\forall\,\underline{u}_{\ell}\in\mathbb{R}^{\dim
V_{\ell}}$ (3.4)

and

$(\tau_{\ell}L^{-1}_{\ell}\mathcal{A}_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})_{\mathcal{A}_{\ell}}\leq(\underline{u}_{\ell},\underline{u}_{\ell})_{\mathcal{A}_{\ell}}\quad\forall\,\underline{u}_{\ell}\in\mathbb{R}^{\dim
V_{\ell}}$ (3.5)

holds for $\ell=1,\ldots,L$.

Then, the estimate

$\left((I-B^{s}_{L}\mathcal{A}_{L})\underline{u}_{L},\underline{u}_{L}\right)_{\mathcal{A}_{L}}\leq\left(1-\frac{1}{CL}\right)\left(\underline{u}_{L},\underline{u}_{L}\right)_{\mathcal{A}_{L}},$

holds for all $\underline{u}_{L}\in\mathbb{R}^{\dim V_{L}}$, where
$C=[1+C_{2}^{1/2}+(C_{S}C_{1})^{1/2}]^{2}$.

For a proof, see [3, Theorem 1].

######  Remark 3.2.

Condition (3.4) is only required for functions $u_{\ell}$ in the range of
$Q_{\ell}-Q_{\ell-1}$. However, since we do not exploit this, we have stated
the stronger condition.

Now, we provide conditions that guarantee (3.4) and (3.5), which fit our needs
better than the original conditions.

######  Lemma 3.1.

If there exists a constant $C_{S}$, independent of $\ell$, which satisfies

$(\mathcal{A}_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})\leq\frac{1}{\tau_{\ell}}(L_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})\leq\lambda_{\ell}C_{S}(X_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})\quad\forall\,\underline{u}_{\ell}\in\mathbb{R}^{\dim
V_{\ell}}$ (3.6)

for each $\ell=1,\ldots,L$. Then, the assumptions (3.4) and (3.5) hold for the
same $C_{S}$.

###### Proof.

We start by showing that the first inequality implies (3.5), i.e., the
smoothing operator $I-\tau_{\ell}L^{-1}_{\ell}\mathcal{A}_{\ell}$ is
nonnegative in $\mathcal{A}_{\ell}$. Let
$\underline{w}_{\ell}\in\mathbb{R}^{\dim V_{\ell}}$ be an arbitrary vector.
Using the Cauchy-Schwarz inequality and the first inequality in (3.6), we
obtain

$\displaystyle\tau_{\ell}(L^{-1}_{\ell}\underline{w}_{\ell},\underline{w}_{\ell})$
$\displaystyle=\tau_{\ell}(\mathcal{A}^{1/2}_{\ell}L^{-1}_{\ell}\underline{w}_{\ell},\mathcal{A}^{-1/2}_{\ell}\underline{w}_{\ell})$
$\displaystyle\leq\tau_{\ell}(\mathcal{A}_{\ell}L^{-1}_{\ell}\underline{w}_{\ell},L^{-1}_{\ell}\underline{w}_{\ell})^{1/2}(\mathcal{A}^{-1}_{\ell}\underline{w}_{\ell},\underline{w}_{\ell})^{1/2}$
$\displaystyle\leq\tau^{1/2}_{\ell}(L^{-1}_{\ell}\underline{w}_{\ell},\underline{w}_{\ell})^{1/2}(\mathcal{A}^{-1}_{\ell}\underline{w}_{\ell},\underline{w}_{\ell})^{1/2}$

It follows that

$\tau_{\ell}(L^{-1}_{\ell}\underline{w}_{\ell},\underline{w}_{\ell})\leq(\mathcal{A}^{-1}_{\ell}\underline{w}_{\ell},\underline{w}_{\ell})\quad\forall\,\underline{w}_{\ell}\in\mathbb{R}^{\dim
V_{\ell}}.$

By substituting $\underline{w}_{\ell}$ with
$\mathcal{A}_{\ell}\underline{u}_{\ell}$, we get (3.5). Next, we use the
Cauchy-Schwarz inequality and the second inequality in (3.6) to show (3.4).
Let $\underline{w}_{\ell}\in\mathbb{R}^{\dim V_{\ell}}$, we have

$\displaystyle(X^{-1}_{\ell}\underline{w}_{\ell},\underline{w}_{\ell})$
$\displaystyle=(L_{\ell}^{1/2}X^{-1}_{\ell}\underline{w}_{\ell},L^{-1/2}_{\ell}\underline{w}_{\ell})\leq(L_{\ell}X^{-1}_{\ell}\underline{w}_{\ell},X^{-1}_{\ell}\underline{w}_{\ell})^{1/2}(L^{-1}_{\ell}\underline{w}_{\ell},\underline{w}_{\ell})^{1/2}$
$\displaystyle\leq\tau^{1/2}_{\ell}\lambda_{\ell}^{1/2}{C}_{S}^{1/2}(X^{-1}_{\ell}\underline{w}_{\ell},\underline{w}_{\ell})^{1/2}(L^{-1}_{\ell}\underline{w}_{\ell},\underline{w}_{\ell})^{1/2}.$

By squaring the inequality, we get

$(X_{\ell}^{-1}\underline{w}_{\ell},\underline{w}_{\ell})=\tau_{\ell}\lambda_{\ell}C_{S}(L_{\ell}^{-1}\underline{w}_{\ell},\underline{w}_{\ell})\quad\forall\,\underline{w}_{\ell}\in\mathbb{R}^{\dim
V_{\ell}}.$

By substituting $\underline{w}_{\ell}$ with $X_{\ell}\underline{u}_{\ell}$, we
get (3.4). ∎

##  4 Approximation error estimates

In this section, we prove some approximation error estimates and provide a
projector which will be used to prove (3.2) and (3.3).

###  4.1 Error and stability estimates for the univariate case

We start by introducing a periodic spline space. For any given sequence of
grid points $\bm{\tau}=(0,\tau_{1},\ldots,\tau_{N},1)$, we define

$\bm{\tau}^{per}:=(-1,-\tau_{N},\cdots,-\tau_{1},0,\tau_{1},\cdots,\tau_{N},1).$

For each $p\in\mathbb{N}$, we define the periodic spline space

$S_{p,\bm{\tau}}^{per}:=\left\\{v\in
S_{p,\bm{\tau}^{per}}\,:\,\partial^{l}v\left(-1\right)=\partial^{l}v\left(1\right)\quad\forall\,l\in\mathbb{N}_{0}\mbox{
with }l<p\right\\}$

and a spline space with vanishing even derivatives on the boundary

$S^{0}_{p,\bm{\tau}}:=\left\\{v\in
S_{p,\bm{\tau}}\,:\,\partial^{2l}v\left(0\right)=\partial^{2l}v\left(0\right)=0\quad\forall\,l\in\mathbb{N}_{0}\mbox{
with }2l<p\right\\}.$ (4.1)

We also define the periodic Sobolev space

$H^{q}_{per}(-1,1):=\left\\{v\in
H^{q}(-1,1)\,:\,\partial^{l}v\left(-1\right)=\partial^{l}v\left(1\right),\quad\forall\,l\in\mathbb{N}_{0}\mbox{
with }l<q\right\\}$

for each $q\in\mathbb{N}$. Let
$\Pi^{per}_{p,\bm{\tau}}:H^{2}_{per}(-1,1)\rightarrow S_{p,\bm{\tau}}^{per}$
be the $H^{2}$-orthogonal projector satisfying

$\displaystyle\begin{split}\left(\partial^{2}\Pi^{per}_{p,\bm{\tau}}u,\partial^{2}v\right)_{L^{2}(-1,1)}&=\left(\partial^{2}u,\partial^{2}v\right)_{L^{2}(-1,1)}\quad\forall\,v\in
S_{p,\bm{\tau}}^{per},\\\
\left(\Pi^{per}_{p,\bm{\tau}}u,1\right)_{L^{2}(-1,1)}&=\left(u,1\right)_{L^{2}(-1,1)}.\end{split}$
(4.2)

We use the following approximation error estimate for spline spaces which does
not require uniform knot spans.

######  Theorem 4.1.

For any $p\geq 3$, we have

$\|\partial^{2}(u-\Pi^{per}_{p,\bm{\tau}}u)\|_{L^{2}(-1,1)}\leq\frac{h_{\bm{\tau}}^{2}}{\pi^{2}}\|\partial^{4}u\|_{L^{2}(-1,1)}\quad\forall\,u\in
H^{4}_{per}(-1,1).$

For a proof, see [22, Theorem 4].

Using the $H^{2}$–$H^{4}$ result above and an Aubin-Nitsche duality trick, we
obtain the following $L^{2}$–$H^{2}$ result.

######  Theorem 4.2.

For any $p\geq 3$, we have

$\|u-\Pi^{per}_{p,\bm{\tau}}u\|_{L^{2}(-1,1)}\leq\frac{h_{\bm{\tau}}^{2}}{\pi^{2}}\|\partial^{2}u\|_{L^{2}(-1,1)}\quad\forall\,u\in
H^{2}_{per}(-1,1).$

###### Proof.

Let $u\in H^{2}_{per}(-1,1)$ be arbitrary but fixed. Let $w\in H^{4}(-1,1)\cap
H^{3}_{per}(-1,1)$ be such that $\partial^{4}w=u-\Pi^{per}_{p,\bm{\tau}}u$.
Note that (4.2) gives
$0=(u-\Pi^{per}_{p,\bm{\tau}}u,1)_{L^{2}(-1,1)}=(\partial^{4}w,1)_{L^{2}(-1,1)}=\partial^{3}w(1)-\partial^{3}w(-1)$.
So, we know that $w\in H^{4}_{per}(-1,1)$.

Using integration by parts (which does not introduce boundary terms since
$u-\Pi^{per}_{p,\bm{\tau}}u\in H^{2}_{per}(-1,1)$ and $w\in
H^{4}_{per}(-1,1)$) and using Theorem 4.1, we obtain

$\displaystyle\|u-\Pi^{per}_{p,\bm{\tau}}u\|^{2}_{L^{2}}$
$\displaystyle=\frac{(u-\Pi^{per}_{p,\bm{\tau}}u,u-\Pi^{per}_{p,\bm{\tau}}u)_{L^{2}}}{\|u-\Pi^{per}_{p,\bm{\tau}}u\|_{L^{2}}}=\frac{(u-\Pi^{per}_{p,\bm{\tau}}u,\partial^{4}w)_{L^{2}}}{\|\partial^{4}w\|_{L^{2}}}$
$\displaystyle=\frac{(\partial^{2}(u-\Pi^{per}_{p,\bm{\tau}}u),\partial^{2}w)_{L^{2}}}{\|\partial^{4}w\|_{L^{2}}}\leq\frac{h_{\bm{\tau}}^{2}}{\pi^{2}}\frac{(\partial^{2}(u-\Pi^{per}_{p,\bm{\tau}}u),\partial^{2}w)_{L^{2}}}{\|\partial^{2}(w-\Pi^{per}_{p,\bm{\tau}}w)\|_{L^{2}}}.$

From the definition of $\Pi^{per}_{p,\bm{\tau}}$, see (4.2), we have
$(\partial^{2}(u-\Pi^{per}_{p,\bm{\tau}}u),\partial^{2}\Pi^{per}_{p,\bm{\tau}}w)_{L^{2}}=0$.
This, together with the Cauchy-Schwarz inequality and the $H^{2}$-stability of
$\Pi^{per}_{p,\bm{\tau}}$, gives

$\displaystyle\|u-\Pi^{per}_{p,\bm{\tau}}u\|^{2}_{L^{2}}$
$\displaystyle\leq\frac{h_{\bm{\tau}}^{2}}{\pi^{2}}\frac{(\partial^{2}(u-\Pi^{per}_{p,\bm{\tau}}u),\partial^{2}(w-\Pi^{per}_{p,\bm{\tau}}w))_{L^{2}}}{\|\partial^{2}(w-\Pi^{per}_{p,\bm{\tau}}w)\|_{L^{2}}}$
$\displaystyle\leq\frac{h_{\bm{\tau}}^{2}}{\pi^{2}}\|\partial^{2}(u-\Pi^{per}_{p,\bm{\tau}}u)\|^{2}_{L^{2}}\leq\frac{h_{\bm{\tau}}^{2}}{\pi^{2}}\|\partial^{2}u\|^{2}_{L^{2}},$

which completes the proof. ∎

Let $\Pi^{0}_{p,\bm{\tau}}:H^{2}(0,1)\cap H^{1}_{0}(0,1)\rightarrow
S_{p,\bm{\tau}}^{0}$ be the $H^{2}$-orthogonal projector satisfying

$\displaystyle\left(\partial^{2}\Pi^{0}_{p,\bm{\tau}}u,\partial^{2}v\right)_{L^{2}(0,1)}$
$\displaystyle=\left(\partial^{2}u,\partial^{2}v\right)_{L^{2}(0,1)}\quad\forall\,v\in
S_{p,\bm{\tau}}^{0}.$

######  Theorem 4.3.

For any $p\geq 3$, we have

$\|u-\Pi^{0}_{p,\bm{\tau}}u\|_{L^{2}(0,1)}\leq\frac{h_{\bm{\tau}}^{2}}{\pi^{2}}\|\partial^{2}u\|_{L^{2}(0,1)}\quad\forall\,u\in
H^{2}(0,1)\cap H^{1}_{0}(0,1).$

###### Proof.

Let $u\in H^{2}(0,1)\cap H^{1}_{0}(0,1)$ be arbitrary but fixed. Define $w$ on
$(-1,1)$ to be

$w(x):=\mbox{sign}(x)\;u(|x|).$

Observe that we obtain $w\in H^{2}_{per}(-1,1)$. From Theorem 4.2, we have

$\|(I-\Pi^{per}_{p,\bm{\tau}})w\|_{L^{2}(-1,1)}\leq
ch_{\bm{\tau}}^{2}\|\partial^{2}w\|_{L^{2}(-1,1)}.$

Observe that
$\|\partial^{2}w\|_{L^{2}(-1,1)}=2^{1/2}\|\partial^{2}u\|_{L^{2}(0,1)}$.
Define $w_{\bm{\tau}}:=\Pi^{per}_{p,\bm{\tau}}w$ and let $u_{\bm{\tau}}$ be
the restriction of $w_{\bm{\tau}}$ to $(0,1)$. Observe that $w_{\bm{\tau}}$ is
anti-symmetric, which implies that $u_{\bm{\tau}}\in S^{0}_{p,\bm{\tau}}$. It
follows that
$\|w-w_{\bm{\tau}}\|_{L^{2}(-1,1)}=2^{1/2}\|u-u_{\bm{\tau}}\|_{L^{2}(0,1)}$.
Using this, we obtain

$\|u-u_{\bm{\tau}}\|_{L^{2}(0,1)}\leq
ch_{\bm{\tau}}^{2}\|\partial^{2}u\|_{L^{2}(0,1)}.$

It remains to show that $u_{\bm{\tau}}$ coincides with
$\Pi^{per}_{p,\bm{\tau}}u$, i.e., to show that $u-u_{\bm{\tau}}$ is
$H^{2}$-orthogonal to $S^{0}_{p,\bm{\tau}}$. By definition, this means that we
have to show

$(\partial^{2}(u-u_{\bm{\tau}}),\partial^{2}\tilde{u}_{\bm{\tau}})_{L^{2}(0,1)}=0\quad\forall\,\tilde{u}_{\bm{\tau}}\in
S^{0}_{p,\bm{\tau}}.$

Let $\tilde{w}_{\bm{\tau}}\in S_{p,\bm{\tau}}^{per}$ be
$\tilde{w}_{\bm{\tau}}:=\mbox{sign}(x)\,\tilde{u}_{\bm{\tau}}(|x|)$ and
observe that
$2(\partial^{2}(u-u_{\bm{\tau}}),\partial^{2}\tilde{u}_{\bm{\tau}})_{L^{2}(0,1)}=(\partial^{2}(w-w_{\bm{\tau}}),\partial^{2}\tilde{w}_{\bm{\tau}})_{L^{2}(0,1)}$,
since $u$, $u_{\bm{\tau}}$ and $\tilde{u}_{\bm{\tau}}$ are restrictions of
$w$, $w_{\bm{\tau}}$ and $\tilde{w}_{\bm{\tau}}$, respectively. Furthermore,
$(\partial^{2}(w-w_{\bm{\tau}}),\partial^{2}\tilde{w}_{\bm{\tau}})_{L^{2}(-1,1)}=0$
by construction, since $w_{\bm{\tau}}:=\Pi^{per}_{p,\bm{\tau}}w$, which
completes the proof. ∎

Let $Q^{0}_{p,\bm{\tau}}:H^{2}(0,1)\cap H^{1}_{0}(0,1)\rightarrow
S_{p,\bm{\tau}}^{0}$ be the $L^{2}$-orthogonal projector satisfying

$\displaystyle\left(Q^{0}_{p,\bm{\tau}}u,v\right)_{L^{2}(0,1)}$
$\displaystyle=\left(u,v\right)_{L^{2}(0,1)}\quad\forall\,v\in
S_{p,\bm{\tau}}^{0}.$

Since the $L^{2}$-orthogonal projector minimizes the error in the
$L^{2}$-norm, Theorem 4.3 immediately implies the following statement.

######  Theorem 4.4.

For any $p\geq 3$, we have

$\|u-Q^{0}_{p,\bm{\tau}}u\|_{L^{2}(0,1)}\leq\frac{h_{\bm{\tau}}^{2}}{\pi^{2}}\|\partial^{2}u\|_{L^{2}(0,1)}\quad\forall\,u\in
H^{2}(0,1)\cap H^{1}_{0}(0,1).$

Next, we show the stability of $Q_{p,\bm{\tau}}^{0}$ with respect to the
$H^{2}$-seminorm. Such a proof is possible since the space
$S_{p,\bm{\tau}}^{0}$ satisfies the following $p$-robust inverse inequality,
while the space $S_{p,\bm{\tau}}\cap H^{1}_{0}(0,1)$ does not satisfy such an
inverse inequality, cf. [28].

######  Theorem 4.5.

Let $p\in\mathbb{N}$ with $p\geq 2$. We have

$\|\partial^{2}u_{\bm{\tau}}\|_{L^{2}(0,1)}\leq
12h^{-2}_{\bm{\tau},\mathrm{min}}\|u_{\bm{\tau}}\|_{L^{2}(0,1)}\quad\forall
u_{\bm{\tau}}\in S^{0}_{p,\bm{\tau}}.$

A proof can be found in [25, Theorem 12].

######  Theorem 4.6.

Let $p\in\mathbb{N}$ with $p\geq 3$. Then there exists a constant $c>0$ such
that

$\displaystyle\|\partial^{2}(Q_{p,\bm{\tau}}^{0}u)\|^{2}_{L^{2}(0,1)}\leq
c\|\partial^{2}u\|^{2}_{L^{2}(0,1)}\quad\forall u\in H^{2}(0,1)\cap
H^{1}_{0}(0,1).$

###### Proof.

The proof is analogous to that of [25, Theorem 14], however it is given here
for completeness. Using the triangle inequality and the inverse inequality, we
obtain

$\displaystyle\|\partial^{2}Q_{p,\bm{\tau}}^{0}u\|^{2}_{L^{2}}$
$\displaystyle\leq
2\|\partial^{2}\Pi_{p,\bm{\tau}}^{0}u\|^{2}_{L^{2}}+2\|\partial^{2}(Q_{p,\bm{\tau}}^{0}u-\Pi_{p,\bm{\tau}}^{0}u)\|^{2}_{L^{2}}$
$\displaystyle\leq
2\|\partial^{2}\Pi_{p,\bm{\tau}}^{0}u\|^{2}_{L^{2}}+ch_{\bm{\tau},\mathrm{min}}^{-2}\|Q_{p,\bm{\tau}}^{0}u-\Pi_{p,\bm{\tau}}^{0}u\|^{2}_{L^{2}}$
$\displaystyle\leq
2\|\partial^{2}\Pi_{p,\bm{\tau}}^{0}u\|^{2}_{L^{2}}+ch_{\bm{\tau},\mathrm{min}}^{-2}\|u-\Pi_{p,\bm{\tau}}^{0}u\|^{2}_{L^{2}}+ch_{\bm{\tau},\mathrm{min}}^{-2}\|u-Q_{p,\bm{\tau}}^{0}u\|^{2}_{L^{2}}.$

The Theorems 4.3 and 4.4 and Assumption (2.5) give the desired result. ∎

###  4.2 Proof of the approximation properties

In this subsection, we consider the discretization framework from Section 2.
We choose

$X_{\ell}:=\mathcal{B}_{\ell}+(\beta+h_{\ell}^{-4})\mathcal{M}_{\ell},$

which corresponds to the norm $\|\cdot\|_{X_{\ell}}$ that satisfies

$\|u\|^{2}_{X_{\ell}}=\|u\|^{2}_{\mathcal{B}}+(\beta+h_{\ell}^{-4})\|u\|^{2}_{L^{2}(\Omega)}\quad\forall\,u\in
V.$

Now, we give a bound for the eigenvalues of $X^{-1}_{\ell}\mathcal{A}_{\ell}$.

######  Lemma 4.1.

Let $\lambda_{\ell}$ with $\ell\geq 1$ be the largest eigenvalue of
$X_{\ell}^{-1}\mathcal{A}_{\ell}$. For $p\geq 3$, we have
$\lambda_{\ell}\in(\frac{1}{1+c},1)$ for some positive constant $c$.

###### Proof.

Since $\mathcal{M}_{\ell}$ is symmetric positive definite and
$h_{\ell}^{-4}>0$, we have $\mathcal{A}_{\ell}<X_{\ell}$, which implies
$\lambda_{\ell}<1$.

For the lower bound, we use $V_{\ell-1}\subsetneqq V_{\ell}$, which implies
that there is some $w_{\ell}\in V_{\ell}$ that is $L_{2}$-orthogonal to
$V_{\ell-1}$, that is $(w_{\ell},u_{\ell-1})_{L^{2}(\Omega)}=0$ for all
$u_{\ell-1}\in V_{\ell-1}$. By combining Theorem 4.7 and (2.4), we obtain

$\|w_{\ell}\|_{L^{2}(\Omega)}=\sup_{u_{\ell-1}\in
V_{\ell-1}}\|w_{\ell}-u_{\ell-1}\|_{L^{2}(\Omega)}\leq
c\,h_{\ell-1}^{2}\|w_{\ell}\|_{\mathcal{B}}.$

In matrix-vector notation, this reads as

$\underline{w}_{\ell}^{T}\mathcal{M}_{\ell}\underline{w}_{\ell}\leq
c\,h_{\ell-1}^{4}\,\underline{w}_{\ell}^{T}\mathcal{B}_{\ell}\underline{w}_{\ell}.$

Using (2.6), we know that there is a constant $c>0$ such that

$\displaystyle\underline{w}_{\ell}^{\top}X_{\ell}\underline{w}_{\ell}$
$\displaystyle=\underline{w}_{\ell}^{\top}\mathcal{A}_{\ell}\underline{w}_{\ell}+h_{\ell}^{-4}\underline{w}_{\ell}^{\top}\mathcal{M}_{\ell}\underline{w}_{\ell}<(1+c)\underline{w}_{\ell}^{\top}\mathcal{A}_{\ell}\underline{w}_{\ell},$

which shows $\lambda_{\ell}>1/(1+c)$. ∎

Next, we prove (3.2) and (3.3). This requires that we choose the projectors
$\mathbf{Q}^{0}_{p,\ell}$, which have to map into the space $V_{\ell}$. We
first define a projector that maps from $\widehat{V}$ into
$\widehat{V}_{\ell}$ by tensorization of the univariate projectors:

$\widehat{\mathbf{Q}}^{0}_{p,\ell}:=Q^{0}_{p,\bm{\tau}_{\ell,1}}\otimes\cdots\otimes
Q^{0}_{p,\bm{\tau}_{\ell,d}},$

where the tensor product is to be understood as in [27, Section 3.2]. The next
two theorems follow from Theorems 4.4 and 4.6 by standard arguments.

######  Theorem 4.7.

Let $p\in\mathbb{N}$ with $p\geq 3$. Then there exists a constant $c$ such
that

$\|(I-\widehat{\mathbf{Q}}_{p,\ell}^{0})\widehat{u}\|_{L^{2}(\widehat{\Omega})}\leq
ch_{\ell}^{2}\|\widehat{u}\|_{\bar{\mathcal{B}}}\quad\forall\,\widehat{u}\in
H^{2}(\widehat{\Omega})\cap H^{1}_{0}(\widehat{\Omega}).$

###### Proof.

The proof is given for the two-dimensional case. We have by definition and
using the triangle inequality

$\displaystyle\|(I-\widehat{\mathbf{Q}}_{p,\ell}^{0})\widehat{u}\|_{L^{2}(\widehat{\Omega})}$
$\displaystyle=\|(I-Q_{p,\bm{\tau}_{\ell,1}}^{0}\otimes
Q_{p,\bm{\tau}_{\ell,2}}^{0})\widehat{u}\|_{L^{2}(\widehat{\Omega})}$
$\displaystyle\leq\|(I-Q_{p,\bm{\tau}_{\ell,1}}^{0}\otimes
I)\widehat{u}\|_{L^{2}(\widehat{\Omega})}+\|(Q_{p,\bm{\tau}_{\ell,1}}^{0}\otimes
I)(I-I\otimes
Q_{p,\bm{\tau}_{\ell,2}}^{0})\widehat{u}\|_{L^{2}(\widehat{\Omega})}.$

Using the $L^{2}$-stability of the $L^{2}$-projectors, we further obtain

$\displaystyle\|(I-\widehat{\mathbf{Q}}_{p,\ell}^{0})\widehat{u}\|_{L^{2}(\widehat{\Omega})}$
$\displaystyle\leq\|(I-Q_{p,\bm{\tau}_{\ell,1}}^{0}\otimes
I)\widehat{u}\|_{L^{2}(\widehat{\Omega})}+\|(I-I\otimes
Q_{p,\bm{\tau}_{\ell,2}}^{0})\widehat{u}\|_{L^{2}(\widehat{\Omega})}.$

The desired result immediately follows from Theorem 4.4. The extension to more
dimensions is obvious. ∎

######  Theorem 4.8.

Let $p\in\mathbb{N}$ with $p\geq 3$. Then there exists a constant $c>0$ such
that

$\displaystyle\|\widehat{\bm{Q}}_{p,\ell}^{0}\widehat{u}\|^{2}_{\bar{\mathcal{B}}}\leq
c\|\widehat{u}\|^{2}_{\bar{\mathcal{B}}}\quad\forall\widehat{u}\in
H^{2}(\widehat{\Omega})\cap H^{1}_{0}(\widehat{\Omega}).$

###### Proof.

The proof is given for the two-dimensional case. We have by definition and
using the triangle inequality

$\displaystyle\|\widehat{\mathbf{Q}}_{p,\ell}^{0}\widehat{u}\|_{\bar{\mathcal{B}}}^{2}$
$\displaystyle=\|\partial_{x_{1}}^{2}(Q_{p,\bm{\tau}_{\ell,1}}^{0}\otimes
Q_{p,\bm{\tau}_{\ell,2}}^{0})\widehat{u}\|_{L^{2}(\widehat{\Omega})}^{2}+\|\partial_{x_{2}}^{2}(Q_{p,\bm{\tau}_{\ell,1}}^{0}\otimes
Q_{p,\bm{\tau}_{\ell,2}}^{0})\widehat{u}\|_{L^{2}(\widehat{\Omega})}^{2}$
$\displaystyle=\|\partial_{x_{1}}^{2}(I\otimes
Q_{p,\bm{\tau}_{\ell,2}}^{0})(Q_{p,\bm{\tau}_{\ell,1}}^{0}\otimes
I)\widehat{u}\|_{L^{2}(\widehat{\Omega})}^{2}+\|\partial_{x_{2}}^{2}(Q_{p,\bm{\tau}_{\ell,1}}^{0}\otimes
I)(I\otimes
Q_{p,\bm{\tau}_{\ell,2}}^{0})\widehat{u}\|_{L^{2}(\widehat{\Omega})}^{2}.$

Using the $L^{2}$-stability of the $L^{2}$-projector, we obtain

$\displaystyle\|\widehat{\mathbf{Q}}_{p,\ell}^{0}\widehat{u}\|_{\bar{\mathcal{B}}}^{2}$
$\displaystyle\leq\|\partial_{x_{1}}^{2}(I\otimes
Q_{p,\bm{\tau}_{\ell,1}}^{0})\widehat{u}\|_{L^{2}(\widehat{\Omega})}^{2}+\|\partial_{x_{2}}^{2}(Q_{p,\bm{\tau}_{\ell,2}}^{0}\otimes
I)\widehat{u}\|_{L^{2}(\widehat{\Omega})}^{2}.$

Using Theorem 4.6, we further obtain

$\displaystyle\|\widehat{\mathbf{Q}}_{p,\ell}^{0}\widehat{u}\|_{\bar{\mathcal{B}}}^{2}$
$\displaystyle\leq
c\|\partial_{x_{1}}^{2}\widehat{u}\|_{L^{2}(\widehat{\Omega})}^{2}+c\|\partial_{x_{2}}^{2}\widehat{u}\|_{L^{2}(\widehat{\Omega})}^{2}=c\|\widehat{u}\|_{\bar{\mathcal{B}}}^{2},$

which finishes the proof. The extension to more dimensions is obvious. ∎

The projectors $\mathbf{Q}^{0}_{p,\ell}$ are now defined via the pull-back
principle, such that

$\mathbf{Q}^{0}_{p,\ell}u:=(\widehat{\mathbf{Q}}^{0}_{p,\ell}(u\circ\bm{G}))\circ\bm{G}^{-1}\quad\forall\,u\in
V.$ (4.3)

Note that, by construction, $\mathbf{Q}^{0}_{p,\ell}$ maps into a subspace of
$V_{\ell}$, where all even outer normal derivatives on the boundary vanish.

######  Theorem 4.9.

Let $d\in\mathbb{N}$ and $p\in\mathbb{N}$ with $p\geq 3$. For each level
$\ell=0,1,\ldots,L-1$, let $\mathbf{Q}^{0}_{p,\ell}:H^{2}(\Omega)\cap
H^{1}_{0}(\Omega)\rightarrow V_{\ell}$ be the projectors defined in (4.3).
There exists a constants $C_{1}$ and $C_{2}$ such that

$\displaystyle\|(\mathbf{Q}^{0}_{p,\ell}-\mathbf{Q}^{0}_{p,\ell-1})u_{L}\|^{2}_{X_{\ell}}$
$\displaystyle\leq C_{1}\lambda^{-1}_{\ell}(u_{L},u_{L})_{\mathcal{A}}\quad$
for $\displaystyle\ell=1,\ldots,L,$ (4.4)
$\displaystyle(\mathbf{Q}^{0}_{p,\ell}\,u_{L},\mathbf{Q}^{0}_{p,\ell}\,u_{L})_{\mathcal{A}}$
$\displaystyle\leq C_{2}(u_{L},u_{L})_{\mathcal{A}}\quad$ for
$\displaystyle\ell=0,\ldots,L-1,$ (4.5)

for all $u_{L}\in V_{L}$.

###### Proof.

Let $u_{L}\in V_{L}$ arbitrary but fixed and let
$\widehat{u}_{L}:=u_{L}\circ\bm{G}\in\widehat{V}_{L}$. Using (2.4), Lemma 2.1
and Theorem 4.8 and the $L^{2}$-stability of
$\widehat{\mathbf{Q}}^{0}_{p,\ell}$, we obtain

$\displaystyle(\mathbf{Q}^{0}_{p,\ell}\,u_{L},\mathbf{Q}^{0}_{p,\ell}\,u_{L})_{\mathcal{A}}$
$\displaystyle\leq
c(\widehat{\mathbf{Q}}^{0}_{p,\ell}\,\widehat{u}_{L},\widehat{\mathbf{Q}}^{0}_{p,\ell}\,\widehat{u}_{L})_{\widehat{\mathcal{A}}}=c\beta\|\widehat{\mathbf{Q}}^{0}_{p,\ell}\widehat{u}_{L}\|^{2}_{L^{2}(\widehat{\Omega})}+c\|\widehat{\mathbf{Q}}^{0}_{p,\ell}\widehat{u}_{L}\|^{2}_{\widehat{\mathcal{B}}}$
$\displaystyle\leq
c\beta\|\widehat{u}_{L}\|^{2}_{L^{2}(\widehat{\Omega})}+c\|\widehat{u}_{L}\|^{2}_{\widehat{\mathcal{B}}}\leq
c(\widehat{u}_{L},\widehat{u}_{L})_{\widehat{\mathcal{A}}}\leq
C_{2}(u_{L},u_{L})_{\mathcal{A}},$

which shows (4.5). Next we prove the auxiliary result

$\|(I-\mathbf{Q}^{0}_{p,\ell-1})u_{L}\|^{2}_{X_{\ell}}\leq
c\lambda^{-1}_{\ell}(u_{L},u_{L})_{\mathcal{A}}\quad\text{for}\quad\ell=1,\ldots,L.$
(4.6)

Using (2.4), (2.1), Theorem 4.8, Theorem 4.7 and the $L^{2}$-stability of
$\mathbf{Q}^{0}_{p,\ell-1}$, we get

$\displaystyle\|(I-\mathbf{Q}^{0}_{p,\ell-1})u_{L}\|^{2}_{X_{\ell}}$
$\displaystyle=\|(I-\mathbf{Q}^{0}_{p,\ell-1})u_{L}\|^{2}_{\mathcal{B}}+(\beta+h_{\ell}^{-4})\|(I-\mathbf{Q}^{0}_{p,\ell-1})u_{L}\|^{2}_{L^{2}(\Omega)}$
$\displaystyle\leq
c\|(I-\widehat{\mathbf{Q}}^{0}_{p,\ell-1})\widehat{u}_{L}\|^{2}_{\bar{\mathcal{B}}}+c(\beta+h_{\ell}^{-4})\|(I-\widehat{\mathbf{Q}}^{0}_{p,\ell-1})\widehat{u}_{L}\|^{2}_{L^{2}(\widehat{\Omega})}$
$\displaystyle\leq
c\|\widehat{u}_{L}\|^{2}_{\bar{\mathcal{B}}}+ch_{\ell}^{-4}h_{\ell-1}^{4}\|\widehat{u}_{L}\|^{2}_{\widehat{\mathcal{B}}}+c\beta\|\widehat{u}_{L}\|^{2}_{L^{2}(\widehat{\Omega})}$
$\displaystyle\leq
c(1+h_{\ell}^{-4}h_{\ell-1}^{4})\|u_{L}\|^{2}_{\mathcal{B}}+c\beta\|u_{L}\|^{2}_{L^{2}(\Omega)}.$

We use assumption (2.6) and Lemma 4.1 to get (4.6). To complete the proof, we
use the fact that
$\mathbf{Q}^{0}_{p,\ell-1}\mathbf{Q}^{0}_{p,\ell}=\mathbf{Q}^{0}_{p,\ell-1}$,
(4.6) and (4.5), to obtain

$\displaystyle\|(\mathbf{Q}^{0}_{p,\ell}-\mathbf{Q}^{0}_{p,\ell-1})u_{L}\|^{2}_{X_{\ell}}$
$\displaystyle=\|(I-\mathbf{Q}^{0}_{p,\ell-1})\mathbf{Q}^{0}_{p,\ell}u_{L}\|^{2}_{X_{\ell}}\leq
c\lambda^{-1}_{\ell}(\mathbf{Q}^{0}_{p,\ell}u_{L},\mathbf{Q}^{0}_{p,\ell}u_{L})_{\mathcal{A}}$
$\displaystyle\leq C_{1}\lambda^{-1}_{\ell}(u_{L},u_{L})_{\mathcal{A}}.$

This shows (4.4) and finishes the proof. ∎

######  Remark 4.1.

In [24, Lemma 9.2], a similar result to Theorem 4.7 is shown. There, the
$\mathcal{B}_{\ell}$-orthogonal projector is considered. That proof only holds
for uniform grids. By using an $L^{2}$-orthogonal projector, we avoid these
difficulties. Since the convergence theory by Hackbusch [12] requires the
error estimates for the $\mathcal{B}_{\ell}$-orthogonal projector, this
motivated us to use the convergence theory by Bramble [2], where this is not
the case.

##  5 The smoothers and the overall convergence results

###  5.1 Subspace corrected mass smoother

We consider the subspace corrected mass smoother, which was originally
proposed in [14] for a second order problem and was one of the first smoothers
to produce a multigrid method for IgA which is robust in both the grid size
and the spline degree. In [25, 24] this smoother was extended to biharmonic
problems. The smoother is based around the inverse inequality in Theorem 4.5,
which is independent of the spline degree.

First, we introduce a splitting for the one dimensional case as follows:

$S_{p,\bm{\tau}}\cap H^{1}_{0}(0,1)=S^{0}_{p,\bm{\tau}}\oplus
S^{1}_{p,\bm{\tau}},$

where $S^{0}_{p,\bm{\tau}}$ is as defined in (4.1) and $S^{1}_{p,\bm{\tau}}$
is its $L^{2}$-orthogonal complement in $S_{p,\bm{\tau}}\cap H^{1}_{0}(0,1)$.
For each of these spaces, we define the corresponding $L^{2}$-orthogonal
projection

$\displaystyle Q_{p,\bm{\tau}}^{0}:H^{2}(0,1)\cap H^{1}_{0}(0,1)\rightarrow
S^{0}_{p,\bm{\tau}},$ $\displaystyle Q_{p,\bm{\tau}}^{1}:H^{2}(0,1)\cap
H^{1}_{0}(0,1)\rightarrow S^{1}_{p,\bm{\tau}}.$

The next step, is to extend the splitting to the multivariate case. Let
$\alpha:=(\alpha_{1},\ldots,\alpha_{d})\in\\{0,1\\}^{d}$ be a multiindex. The
tensor product B-spline space $\widehat{V}_{\ell}=S_{p,\bm{\tau}_{\ell}}\cap
H^{1}_{0}(\widehat{\Omega})$ with
$\bm{\tau}_{\ell}=(\bm{\tau}_{\ell,1},\ldots,\bm{\tau}_{\ell,d})$ is split
into the direct sum of $2^{d}$ subspaces

$\widehat{V}_{\ell}=\bigoplus_{\alpha\in\\{0,1\\}^{d}}S^{\alpha}_{p,\bm{\tau}_{\ell}}\quad\text{where}\quad
S^{\alpha}_{p,\bm{\tau}_{\ell}}=S^{\alpha_{1}}_{p,\bm{\tau}_{\ell,1}}\otimes\cdots\otimes
S^{\alpha_{d}}_{p,\bm{\tau}_{\ell,d}}.$ (5.1)

Again, we define $L^{2}$-orthogonal projectors

$\widehat{\mathbf{Q}}_{p,\bm{\tau}_{\ell}}^{\alpha}:=Q_{p,\bm{\tau}_{\ell,1}}^{\alpha_{1}}\otimes\cdots\otimes
Q_{p,\bm{\tau}_{\ell,d}}^{\alpha_{d}}:\widehat{V}\rightarrow
S^{\alpha}_{p,\bm{\tau}_{\ell}}.$

The projector $\widehat{\mathbf{Q}}_{p,\bm{\tau}_{\ell}}^{0}$ from Section 4.2
is consistent with this definition, for the choice $\alpha=0$. Since the
splitting is $L^{2}$-orthogonal, we obviously have the following result.

$\widehat{u}_{\ell}=\sum_{\alpha\in\\{0,1\\}^{d}}\widehat{\mathbf{Q}}_{p,\bm{\tau}_{\ell}}^{\alpha}\widehat{u}_{\ell}\quad\mbox{and}\quad\|\widehat{u}_{\ell}\|^{2}_{L^{2}(\widehat{\Omega})}=\sum_{\alpha\in\\{0,1\\}^{d}}\|\widehat{\mathbf{Q}}_{p,\bm{\tau}_{\ell}}^{\alpha}\widehat{u}_{\ell}\|^{2}_{L^{2}(\widehat{\Omega})}\quad\forall\widehat{u}_{\ell}\in\widehat{V}_{\ell}.$
(5.2)

The next theorem shows that the splitting is also stable in $H^{2}$.

######  Theorem 5.1.

Let $p\in\mathbb{N}$ with $p\geq 3$. Then there exists a constant $c>0$ such
that

$\displaystyle
c^{-1}\|\widehat{u}_{\ell}\|^{2}_{\bar{\mathcal{B}}}\leq\sum_{\alpha\in\\{0,1\\}^{d}}\|\widehat{\mathbf{Q}}_{p,\bm{\tau}_{\ell}}^{\alpha}\widehat{u}_{\ell}\|^{2}_{\bar{\mathcal{B}}}\leq
c\|\widehat{u}_{\ell}\|^{2}_{\bar{\mathcal{B}}}\quad\forall\widehat{u}_{\ell}\in\widehat{V}_{\ell}.$

###### Proof.

Theorem 4.6 states the stability of $Q^{0}_{p,\bm{\tau}_{\ell}}$ in the
$H^{2}$-seminorm. The stability of $Q^{1}_{p,\bm{\tau}_{\ell}}$ in the
$H^{2}$-seminorm follows using the triangle inequality. The stability of these
statements in the $L^{2}$-norm is obvious. From these observations, the right
inequality follows by arguments that are completely analogous to those of the
proof of Theorem 4.8.

The left inequality follows from (5.2) and the triangle inequality. ∎

For notational convenience, we restrict the setup of the smoother to the two
dimensional case. For notational convenience, we write the splitting (5.1) as

$\widehat{V}_{\ell}=S^{00}_{p,\bm{\tau}_{\ell}}\oplus
S^{01}_{p,\bm{\tau}_{\ell}}\oplus S^{10}_{p,\bm{\tau}_{\ell}}\oplus
S^{11}_{p,\bm{\tau}_{\ell}},\quad\text{where}\quad
S^{\alpha_{1},\alpha_{2}}_{p,\bm{\tau}_{\ell}}=S^{\alpha_{1}}_{p,\bm{\tau}_{\ell,1}}\otimes
S^{\alpha_{2}}_{p,\bm{\tau}_{\ell,2}}.$

Following the ideas of [14, 25], we construct local smoothers $L_{\alpha}$ for
any of the spaces $V_{\ell,\alpha}:=S^{\alpha}_{p,\bm{\tau}_{\ell}}$. These
local contributions are chosen such that they satisfy the corresponding local
condition

$\bar{\mathcal{B}}_{\ell,\alpha}+\beta\widehat{\mathcal{M}}_{\ell,\alpha}\leq
L_{\ell,\alpha}\leq
c(\bar{\mathcal{B}}_{\ell,\alpha}+(\beta+h^{-4})\widehat{\mathcal{M}}_{\ell,\alpha}),$

where

$\bar{\mathcal{B}}_{\ell,\alpha}:=\mathbf{P}_{\ell,\alpha}^{T}\bar{\mathcal{B}}_{\ell}\mathbf{P}_{\ell,\alpha}\qquad\mbox{and}\qquad\widehat{\mathcal{M}}_{\ell,\alpha}:=\mathbf{P}_{\ell,\alpha}^{T}\widehat{\mathcal{M}}_{\ell}\mathbf{P}_{\ell,\alpha}$

and $\mathbf{P}_{\ell,\alpha}$ is the matrix representation of the canonical
embedding $V_{\ell,\alpha}\rightarrow V_{\ell}$. The canonical embedding has
tensor product structure, i.e., $P_{\ell,\alpha_{1}}\otimes\cdots\otimes
P_{\ell,\alpha_{d}}$, where the $P_{\ell,\alpha_{i}}$ are the matrix
representations of the corresponding univariate embeddings. In the two-
dimensional case, $\bar{\mathcal{B}}_{\ell}$ and
$\widehat{\mathcal{M}}_{\ell}$ have the representation

$\bar{\mathcal{B}}_{\ell}=B\otimes M+M\otimes
B\quad\mbox{and}\quad\widehat{\mathcal{M}}_{\ell}=M\otimes M,$

where $B$ and $M$ are the corresponding univariate stiffness and mass matrices
(not necessarily equal for both spacial directions). For notational
convenience, we do not indicate the spacial direction and the grid level for
these matrices. Restricting $\bar{\mathcal{B}}_{\ell}$ to the subspace
$V_{\ell,(\alpha_{1},\alpha_{2})}$ gives

$\bar{\mathcal{B}}_{\ell,(\alpha_{1},\alpha_{2})}=B_{\alpha_{1}}\otimes
M_{\alpha_{2}}+M_{\alpha_{1}}\otimes B_{\alpha_{2}},$

where $B_{\alpha_{i}}=P^{T}_{\ell,\alpha_{i}}BP_{\ell,\alpha_{i}}$ and
$M_{\alpha_{i}}=P^{T}_{\ell,\alpha_{i}}MP_{\ell,\alpha_{i}}$. We define

$\bar{\mathcal{A}}_{\ell}:=\bar{\mathcal{B}}_{\ell}+\beta\widehat{\mathcal{M}}_{\ell}\quad\text{and}\quad\bar{\mathcal{A}}_{\alpha_{1},\alpha_{2}}:=\bar{\mathcal{B}}_{\alpha_{1},\alpha_{2}}+\beta\widehat{\mathcal{M}}_{\alpha_{1},\alpha_{2}}.$

The inverse inequality for $S^{0}_{p,\bm{\tau}_{\ell,i}}$ (Theorem 4.5),
allows us to estimate

$B_{0}\leq\sigma M_{0},$

where $\sigma=\sigma_{0}h_{\ell,\mathrm{min}}^{-4}$ and $\sigma_{0}=144$.
Using this, we define the smoothers $L_{\alpha_{1},\alpha_{2}}$ as follows and
obtain estimates for them as follows:

$\displaystyle\bar{\mathcal{A}}_{00}$
$\displaystyle\leq(2\sigma+\beta)M_{0}\otimes M_{0}$
$\displaystyle=:L_{00}\leq
c(\bar{\mathcal{A}}_{00}+h^{-4}\widehat{\mathcal{M}}_{00}),$
$\displaystyle\bar{\mathcal{A}}_{01}$ $\displaystyle\leq
M_{0}\otimes\left((\sigma+\beta)M_{1}+B_{1}\right)$ $\displaystyle=:L_{01}\leq
c(\bar{\mathcal{A}}_{01}+h^{-4}\widehat{\mathcal{M}}_{01}),$
$\displaystyle\bar{\mathcal{A}}_{10}$
$\displaystyle\leq\left(B_{1}+(\sigma+\beta)M_{1}\right)\otimes M_{0}$
$\displaystyle=:L_{10}\leq
c(\bar{\mathcal{A}}_{10}+h^{-4}\widehat{\mathcal{M}}_{10}),$
$\displaystyle\bar{\mathcal{A}}_{11}$ $\displaystyle=B_{1}\otimes
M_{1}+M_{1}\otimes B_{1}+\beta M_{1}\otimes M_{1}$ $\displaystyle=:L_{11}\leq
c(\bar{\mathcal{A}}_{11}+h^{-4}\widehat{\mathcal{M}}_{11}).$

The extension to three and more dimensions is completely straight-forward (cf.
[14]). For each of the subspaces $V_{\ell,\alpha}$, we have defined a
symmetric and positive definite smoother $L_{\alpha}$. The overall smoother is
given by

$L_{\ell}:=\sum_{\alpha\in\\{0,1\\}^{d}}(\mathbf{Q}^{D,\alpha})^{T}L_{\alpha}\mathbf{Q}^{D,\alpha},$

where
$\mathbf{Q}^{D,\alpha}=\widehat{\mathcal{M}}_{\alpha}^{-1}\mathbf{P}^{T}_{\ell,\alpha}\widehat{\mathcal{M}}_{\ell}$
is the matrix representation of the $L^{2}$-projection from $V_{\ell}$ to
$V_{\ell,\alpha}$. Completely analogous to [14, Section 5.2], we obtain

$L_{\ell}^{-1}=\sum_{\alpha\in\\{0,1\\}^{d}}\mathbf{P}_{\ell,\alpha}L^{-1}_{\alpha}\mathbf{P}^{T}_{\ell,\alpha}.$

######  Theorem 5.2.

Let $d\in\mathbb{N}$ and $p\in\mathbb{N}$ with $p\geq 3$. The subspace
corrected mass smoother $L_{\ell}$, satisfies (3.6), i.e.,

$(\mathcal{A}_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})\leq\frac{1}{\tau_{\ell}}(L_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})\leq
C_{S}\,\lambda_{\ell}\,((\mathcal{A}_{\ell}+h^{-4}\mathcal{M}_{\ell})\underline{u}_{\ell},\underline{u}_{\ell})\quad\forall\,\underline{u}_{\ell}\in\mathbb{R}^{\dim
V_{\ell}}$

for all $\tau\in(0,\tau_{0})$, where $\tau_{0}>0$ is some constant.

###### Proof.

The inequality

$(\bar{\mathcal{A}}_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})\leq(L_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})\leq
c((\bar{\mathcal{A}}_{\ell}+h^{-4}\widehat{\mathcal{M}}_{\ell})\underline{u}_{\ell},\underline{u}_{\ell})$

was shown in [25, Theorem 17] for $\beta=0$. Note that no part of that proof
requires uniform grids. So, the proof can be used almost verbatim also in the
context of this paper. Using (5.2), the extension to $\beta>0$ is straight
forward. Using this and Lemma 2.1, we get

$\displaystyle(\widehat{\mathcal{A}}_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})$
$\displaystyle\leq
d(\bar{\mathcal{A}}_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})\leq\frac{d}{\tau_{\ell}}(L_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})\leq
c((\widehat{\mathcal{A}}_{\ell}+h^{-4}\widehat{\mathcal{M}}_{\ell})\underline{u}_{\ell},\underline{u}_{\ell})$

for some constant $c>0$. Using (2.4), we obtain

$\displaystyle({\mathcal{A}}_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})$
$\displaystyle\leq\frac{c_{1}}{\tau_{\ell}}(L_{\ell}\underline{u}_{\ell},\underline{u}_{\ell})\leq
c_{2}((\mathcal{A}_{\ell}+h^{-4}\mathcal{M}_{\ell})\underline{u}_{\ell},\underline{u}_{\ell})$

for some constants $c_{1},c_{2}>0$, which finishes the proof since
$\lambda_{\ell}$ is bounded from below by a constant (Lemma 4.1). ∎

######  Corollary 5.1.

Suppose that we solve the linear system (2.8) using a multigrid solver as
outlined in Section 3 and using the subspace corrected mass smoother as
outlined in Section 5, then the convergence of the multigrid solver is
described by the relation

$\left((I-B^{s}_{L}\mathcal{A}_{L})\underline{u}_{L},\underline{u}_{L}\right)_{\mathcal{A}_{L}}\leq\left(1-\frac{1}{CL}\right)\left(\underline{u}_{L},\underline{u}_{L}\right)_{\mathcal{A}_{L}},$
(5.3)

where the constant $C$ is independent of the grid sizes $h_{\ell}$, the number
of levels $L$, the spline degree $p$ and the choice of the scaling parameter
$\beta$. It may depend on $d$, the constants $c_{1}$, $c_{2}$, $c_{q}$, and
$c_{r}$ and the shape of $\Omega$, cf. Notation 2.1.

###### Proof.

We use Theorem 3.1, whose assumptions are shown by Theorem 4.9 and the
combination of Lemma 3.1 and Theorem 5.2. ∎

######  Remark 5.1.

The operator $L^{-1}_{\ell}$ can be applied efficiently because all of the
local contributions $L_{00}$, $L_{01}$ and $L_{10}$ can be inverted
efficiently because they are tensor products. For example, we have
$L_{00}^{-1}=\frac{1}{2\sigma+\beta}(M^{-1}_{0}\otimes I)(I\otimes
M^{-1}_{0})$, where both $M^{-1}_{0}\otimes I$ and $I\otimes M^{-1}_{0}$ can
be realized by applying direct solvers for the univariate mass matrix to
several right-hand sides. The operator $L_{11}$ is the sum of two tensor
products. So, it has to be inverted as a whole. However, the dimension of the
corresponding space is so small that the corresponding computational costs are
negligible. More details on how to realize the smoother computationally
efficient, are given in [14, Section 5]. There, it is outlined where an
efficient realization of the subspace corrected mass smoother is also possible
in case of more than two dimensions.

###  5.2 Symmetric Gauss-Seidel smoother and a hybrid smoother

The second smoother we consider is a symmetric Gauss-Seidel smoother
consisting of one forward sweep and one backward sweep. It can be shown that
this smoother satisfies Condition (3.6), where the constant $C_{S}$ depends on
the spline degree, see [25]. This means that also the overall convergence
result (5.3) holds, where again $C$ depends on the spline degree. The
symmetric Gauss-Seidel smoother works well for domains with a nontrivial
geometry transformations, but degenerated for large spline degrees (cf. [8,
16]).

Since the symmetric Gauss-Seidel smoother works well for nontrivial geometry
transformations and the subspace corrected mass smoother is robust with
respect to the spline degree, we combine these smoothers into a hybrid
smoother, which was first introduced in [25]. This hybrid smoother consists of
one forward Gauss-Seidel sweep, followed by one step of the subspace corrected
mass smoother, finally followed by one backward Gauss-Seidel sweep.

##  6 Numerical experiments

In this section, we present the results of numerical experiments performed
with the proposed algorithm. As computational domains, we first consider the
unit square, then we consider the nontrivial geometries displayed in Figures 2
(two dimensional domain) and 2 (three-dimensional domain). We consider the
problem

$\displaystyle\begin{split}\beta u+\Delta^{2}u&=f\quad\text{in}\quad\Omega,\\\
u&=g_{1}\quad\text{on}\quad\partial\Omega,\\\ \Delta
u&=g_{2}\quad\text{on}\quad\partial\Omega,\end{split}$

where

$\displaystyle f(x)=(\beta+d^{2}\pi^{4})\prod^{d}_{k=1}\sin(\pi x_{k}),\quad
g_{1}(x)=\prod^{d}_{k=1}\sin(\pi x_{k}),\quad
g_{2}(x)=-d\pi^{2}\prod^{d}_{k=1}\sin(\pi x_{k}).$

The discretization space on the parameter domain is the space of tensor-
product B-splines. On the coarsest level ($\ell=0$), we choose

$\bm{\tau}_{0,i}=(0,\,1/3,\,1/2,\,4/5,\,1),$ (6.1)

for all spacial directions $i=1,\ldots,d$. The discretization on level $\ell$
is obtained by preforming $\ell$ uniform $h$-refinement steps. The spline
spaces have maximum continuity and spline degree $p$. We solve the resulting
system using the preconditioned conjugate gradient (PCG) with a V-cycle
multigrid method with 1 pre and 1 post smoothing step, as preconditioner. A
random initial guess is used and the stopping criteria is

$\|\underline{r}^{(k)}_{L}\|\leq 10^{-8}\|\underline{r}^{(0)}_{L}\|,$

where
$\underline{r}^{(k)}_{L}:=\underline{f}_{L}-\mathcal{A}_{L}\underline{x}^{(k)}_{L}$
is the residual at step $k$ and $\|\cdot\|$ denotes the Euclidean norm. All
numerical experiments are implemented using the G+Smo library [17].

###  6.1 Numerical experiments on parameter domain

We start with the unit square as the domain, that is, $\Omega=(0,1)^{2}$. Note
that $g_{1}(x)=g_{2}(x)=0$ for this domain. For now, we consider the symmetric
Gauss-Seidel smoother and the subspace corrected mass smoother. For both
smoothers, we choose $\tau=1$. The iteration counts are displayed in Table 1
for $\beta=1$, and in Table 2 for $\beta=10^{7}$.

$\ell\;\diagdown\;p$ |  3 |  4 |  5 |  6 |  7 |  8 |  9
---|---|---|---|---|---|---|---
Symmetric Gauss-Seidel
5 | 10 | 16 | 28 | 45 | 71 | 120 | 210
6 | 10 | 16 | 27 | 44 | 71 | 119 | 209
7 | 10 | 16 | 27 | 44 | 72 | 117 | 212
8 | 11 | 16 | 27 | 45 | 72 | 120 | 221
Subspace corrected mass smoother, $\sigma^{-1}_{0}=0.02$
e 5 | 126 | 122 | 114 | 105 | 98 | 93 | 85
6 | 131 | 129 | 123 | 116 | 110 | 105 | 100
7 | 132 | 133 | 127 | 121 | 116 | 110 | 106
8 | 133 | 134 | 130 | 124 | 118 | 114 | 110
Table 1: Iteration counts for 2D parametric domain, $\beta=1$ $\ell\;\diagdown\;p$ |  3 |  4 |  5 |  6 |  7 |  8 |  9
---|---|---|---|---|---|---|---
Symmetric Gauss-Seidel
5 | 10 | 16 | 28 | 45 | 71 | 119 | 211
6 | 10 | 16 | 27 | 44 | 71 | 118 | 208
7 | 10 | 16 | 27 | 44 | 72 | 117 | 212
8 | 11 | 16 | 27 | 45 | 72 | 119 | 221
Subspace corrected mass smoother, $\sigma^{-1}_{0}=0.02$
5 | 124 | 121 | 113 | 104 | 96 | 92 | 85
6 | 131 | 129 | 123 | 116 | 110 | 105 | 99
7 | 132 | 133 | 127 | 120 | 116 | 110 | 106
8 | 133 | 134 | 130 | 124 | 116 | 118 | 114
Table 2: Iteration counts for 2D parametric domain, $\beta=10^{7}$

From the tables, we see that the symmetric Gauss-Seidel smoother preforms well
for small spline degrees, but degenerates for larger spline degrees. These
results are not surprising since it is known that standard smoothers do not
work well for large spline degrees (cf. [8, 16]). Due to Corollary 5.1, the
multigrid solver with subspace corrected mass smoother is robust with respect
to the spline degree. The tables do reflex this. However, the iteration
numbers are relatively high. Table 3 shows the iteration numbers when using an
uniform grid with spacing $1/4$ on the coarsest level ($\ell=0$), rather than
the grid (6.1). The numbers in Table 3 are significantly smaller. This implies
that the subspace corrected mass smoother is sensitive to the quasi-uniformity
constant $c_{q}$.

$\ell\;\diagdown\;p$ |  3 |  4 |  5 |  6 |  7 |  8 |  9
---|---|---|---|---|---|---|---
Subspace corrected mass smoother, $\sigma^{-1}_{0}=0.015$
5 | 41 | 40 | 39 | 37 | 35 | 34 | 33
6 | 41 | 41 | 39 | 37 | 36 | 35 | 34
7 | 42 | 42 | 40 | 39 | 37 | 35 | 35
8 | 42 | 42 | 41 | 39 | 37 | 37 | 35
Table 3: Iteration counts for 2D parametric domain with uniform grid,
$\beta=1$

###  6.2 Numerical experiments on physical domain

Now, we consider a domain with a nontrivial geometry transformation as
displayed in Figures 2 and 2. The convergence of subspace corrected mass
smoother degrades significantly due to the nontrivial geometry mapping. To
combat this, we consider the hybrid smoother described in Section 5.2.

Figure 1: The two-dimensional domain

Figure 2: The three-dimensional domain

Table 4 and Table 5 display the iteration numbers for the 2D and 3D physical
domains, respectively. These iteration numbers are relatively small and seam
to be robust with respect to both grid size and spline degree. Although the
hybrid smoother is more expensive, as one smoothing step can be view as two
smoothing steps, the reduction of iteration numbers outweigh this cost for
larger spline degrees $p>4$. For smaller spline degrees, the symmetric Gauss-
Seidel smoother is ideal choice.

$\ell\;\diagdown\;p$ |  3 |  4 |  5 |  6 |  7 |  8 |  9
---|---|---|---|---|---|---|---
Hybrid smoother, $\beta=1$
5 | 28 | 23 | 23 | 24 | 26 | 27 | 27
6 | 28 | 23 | 22 | 25 | 24 | 26 | 26
7 | 29 | 23 | 22 | 23 | 24 | 24 | 24
8 | 28 | 22 | 21 | 21 | 22 | 22 | 22
Hybrid smoother, $\beta=10^{7}$
5 | 27 | 23 | 23 | 24 | 26 | 27 | 28
6 | 28 | 23 | 22 | 25 | 25 | 26 | 26
7 | 29 | 23 | 22 | 23 | 24 | 24 | 24
8 | 28 | 22 | 21 | 21 | 22 | 22 | 22
Table 4: Iteration counts for 2D Physical domain, $\sigma^{-1}_{0}=0.015$, $\tau=0.1$ $\ell\;\diagdown\;p$ |  3 |  4 |  5 |  6 |  7
---|---|---|---|---|---
Hybrid smoother, $\beta=1$
1 | 16 | 18 | 21 | 27 | 30
2 | 31 | 28 | 26 | 29 | 32
3 | 46 | 37 | 33 | 33 | 35
4 | 50 | 41 | 34 | 34 | mem
Hybrid smoother, $\beta=10^{7}$
1 | 10 | 11 | 13 | 17 | 20
2 | 12 | 16 | 20 | 25 | 29
3 | 16 | 19 | 22 | 24 | 28
4 | 29 | 28 | 28 | 28 | mem
Table 5: Iteration counts for 3D Physical domain, $\sigma^{-1}_{0}=0.020$,
$\tau=0.1$

######  Remark 6.1.

All experiments have also been performed for the choice $\beta=0$. In this
case, one obtains iteration numbers that are identical than those obtained for
$\beta=1$. Therefore, we chose to only display the results for $\beta=1$.

Acknowledgements. This research was funded by the Austrian Science Fund (FWF):
P31048.

## References

  * [1] A. Beigl, J. Sogn, and W. Zulehner, Robust preconditioners for multiple saddle point problems and applications to optimal control problems, SIAM Journal on Matrix Analysis and Applications, 41 (2020), pp. 1590–1615.
  * [2] J. H. Bramble, Multigrid methods, Routledge, 2018.
  * [3] J. H. Bramble, J. E. Pasciak, J. P. Wang, and J. Xu, Convergence estimates for multigrid algorithms without regularity assumptions, Mathematics of Computation, 57 (1991), pp. 23–45.
  * [4] S. C. Brenner, An optimal-order nonconforming multigrid method for the biharmonic equation, SIAM Journal on Numerical Analysis, 26 (1989), pp. 1124–1138.
  * [5] L. Chen, J. Hu, and X. Huang, Multigrid methods for Hellan–Herrmann–Johnson Mixed Method of Kirchhoff Plate Bending Problems, Journal of Scientific Computing, (2015), pp. 1–24.
  * [6] P. G. Ciarlet, The finite element method for elliptic problems, SIAM, 2002.
  * [7] A. P. de la Riva, F. J. Gaspar, and C. Rodrigo, On the robust solution of an isogeometric discretization of bilaplacian equation by using multigrid methods, Computers & Mathematics with Applications, 80 (2020), pp. 386–394.
  * [8] K. P. S. Gahalaut, J. K. Kraus, and S. K. Tomar, Multigrid methods for isogeometric discretization, Computer methods in applied mechanics and engineering, 253 (2013), pp. 413–425.
  * [9] V. Girault and P.-A. Raviart, Finite element methods for Navier-Stokes equations: theory and algorithms, vol. 5, Springer Science & Business Media, 2012.
  * [10] P. Grisvard, Singularities in boundary value problems, Springer, Berlin, 1992.
  * [11] P. Grisvard, Elliptic Problems in Nonsmooth Domains. Reprint of the 1985 hardback ed., Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM), 2011.
  * [12] W. Hackbusch, Multi-Grid Methods and Applications, vol. 4, Springer Science & Business Media, 2013.
  * [13] M. R. Hanisch, Multigrid preconditioning for the biharmonic Dirichlet problem, SIAM Journal on Numerical Analysis, 30 (1993), pp. 184–214.
  * [14] C. Hofreither and S. Takacs, Robust multigrid for isogeometric analysis based on stable splittings of spline spaces, SIAM Journal on Numerical Analysis, 4 (2017), pp. 2004–2024.
  * [15] C. Hofreither, S. Takacs, and W. Zulehner, A robust multigrid method for isogeometric analysis in two dimensions using boundary correction, Computer Methods in Applied Mechanics and Engineering, 316 (2017), pp. 22–42.
  * [16] C. Hofreither and W. Zulehner, Spectral analysis of geometric multigrid methods for isogeometric analysis, in International Conference on Numerical Methods and Applications, Springer, 2014, pp. 123–129.
  * [17] A. Mantzaflaris, J. Sogn, S. Takacs, and others (see website), G+Smo.  https://github.com/gismo/gismo, 2021.
  * [18] K.-A. Mardal, B. F. Nielsen, and M. Nordaas, Robust preconditioners for PDE-constrained optimization with limited observations, BIT, 57 (2017), pp. 405–431.
  * [19] K.-A. Mardal, J. Sogn, and S. Takacs, Robust preconditioning and error estimates for optimal control of the convection-diffusion-reaction equation with limited observation in isogeometric analysis, SIAM Journal on Numerical Analysis, (2022).  To appear. (Preprint available: https://arxiv.org/abs/2012.13003).
  * [20] M. A. Olshanskii and A. Reusken, On the convergence of a multigrid method for linear reaction-diffusion problems, Computing, 65 (2000), pp. 193–202.
  * [21] K. Rafetseder and W. Zulehner, A decomposition result for Kirchhoff plate bending problems and a new discretization approach, SIAM Journal on Numerical Analysis, 56 (2018), pp. 1961–1986.
  * [22] E. Sande, C. Manni, and H. Speleers, Sharp error estimates for spline approximation: Explicit constants, $n$-widths, and eigenfunction convergence, Mathematical Models and Methods in Applied Sciences, 29 (2019), pp. 1175–1205.
  * [23] E. Sande, C. Manni, and H. Speleers, Explicit error estimates for spline approximation of arbitrary smoothness in isogeometric analysis, Numerische Mathematik, (2020), pp. 1–41.
  * [24] J. Sogn, Schur complement preconditioners for multiple saddle point problems and applications, PhD thesis, Johannes Kepler University Linz, 2018\.
  * [25] J. Sogn and S. Takacs, Robust multigrid solvers for the biharmonic problem in isogeometric analysis, Computers & Mathematics with Applications, 77 (2019), pp. 105–124.
  * [26] J. Sogn and W. Zulehner, Schur complement preconditioners for multiple saddle point problems of block tridiagonal form with application to optimization problems, IMA Journal of Numerical Analysis, 39 (2018), pp. 1328–1359.
  * [27] S. Takacs, Robust approximation error estimates and multigrid solvers for isogeometric multi-patch discretizations, Mathematical Models and Methods in Applied Sciences, 28 (2018), pp. 1899–1928.
  * [28] S. Takacs and T. Takacs, Approximation error estimates and inverse inequalities for B-splines of maximum smoothness, Mathematical Models and Methods in Applied Sciences, 26 (2016), pp. 1411–1445.
  * [29] S. Zhang and J. Xu, Optimal solvers for fourth-order pdes discretized on unstructured grids, SIAM Journal on Numerical Analysis, 52 (2014), pp. 282–307.