#  Space-time least squares approximation for Schrödinger equation and
efficient solver

Andrea Bressan<EMAIL_ADDRESS>Alen Kushova
<EMAIL_ADDRESS>Giancarlo Sangalli
<EMAIL_ADDRESS>Mattia Tani<EMAIL_ADDRESS>Dipartimento di
Matematica “F. Casorati”, Università di Pavia, Via A. Ferrata, 5, 27100 Pavia,
Italy.  Istituto di Matematica Applicata e Tecnologie Informatiche, “E.
Magenes” del CNR, Via A. Ferrata, 1, 27100 Pavia, Italy.

###### Abstract

In this work we present a space-time least squares isogeometric discretization
of the Schrödinger equation and propose a preconditioner for the arising
linear system in the parametric domain. Exploiting the tensor product
structure of the basis functions, the preconditioner is written as the sum of
Kronecker products of matrices. Thanks to an extension to the classical Fast
Diagonalization method, the application of the preconditioner is efficient and
robust w.r.t. the polynomial degree of the spline space. The time required for
the application is almost proportional to the number of degrees-of-freedom,
for a serial execution.

###### keywords:

Isogeometric Analysis , splines , Schrödinger equation , space-time Least
squares formulation , Fast Diagonalization.

##  1 Introduction

Space-time finite element methods originated in the papers [1, 2, 3, 4], where
standard finite elements are ascribed an extra dimension for the time and,
typically, adopt a discontinuous approximation in time, since this produces a
time marching algorithm with a traditional step-by-step format (see e.g. [5]).
Over the years, the theory of space-time methods has been developed mainly for
evolutionary equations of the parabolic type and hyperbolic type, whereas, for
quantum mechanics, and more precisely for Schrödinger’s equation, there are
few contributions and the methods are still in development.

To our knowledge, one of the first works concerning space-time variational
formulations for the (nonlinear) Schrödinger equation, is [6], in which
Karakashian and Makridakis proposed a space–time method combining a conforming
Galerkin discretization in space and an upwind DG time-stepping. This method
reduces to a Radau IIA Runge-Kutta time discretization in the case of constant
potentials. In [7], for the linear Schrödinger equation the authors propose
two variational formulations that are proved to be well posed: a strong
formulation, with no relaxation of the original equation, and an ultraweak
formulation, that transfers all derivatives onto test functions. The proposed
discretization for the ultraweak form is based on a discontinuous Petrov-
Galerkin (DPG) method, that addresses optimal stability, and quasi-optimal
error rates in $L^{2}$-norm. In [8] a space–time ultraweak Trefftz
discontinuous Galerkin (DG) method for the Schrödinger equation with
piecewise-constant potential is proposed and analyzed, proving well-posedness
and stability of the method, and optimal high-order $h$-convergence error
estimates in a skeleton norm, for the one and two dimensional cases. Recently,
in [9], Hain and Urban proposed a space–time ultraweak variational formulation
with optimal inf-sup constant. The formulation in [9] is related to the
ultraweak DPG method in [7], but differs in the choice of the test and trial
spaces. Hain and Urban first fix a conforming test space, and then construct
an optimal trial space, while Demkowicz et al. first constructs a trial space
and then a suitable test space. The discretization proposed in [9] uses high
order B-splines with maximum regularity and can be extended to the
Isogeometric Analysis (IgA) framework.

Introduced in [10], see also the book [11], IgA, is an evolution of the
classical finite element methods. In IgA, both the approximation of the
solution of the partial differential equation that models the problem, and the
representation of the computational domain, are accomplished using B-spline
functions, or their generalizations (NURBS). This is meant to simplify the
interoperability between computer aided design and numerical simulations. IgA
also benefits from the approximation properties of splines, whose high-
continuity yields higher accuracy when compared to $C^{0}$ piecewise
polynomials, see e.g., [12, 13, 14].

In this paper we focus on the linear time dependent Schrödinger equation
without potential. Starting from the well posed space-time strong formulations
in [7], we derive a well posed space-time isogeometric Petrov-Galerkin
discretization, that is essentially a Galerkin approximation of the space-time
least squares variational formulation of the model problem. The matrix
associated to the discrete linear system can be written as sum of Kronecker
products, and has the same structure of the one arising from [9]. The main
contribution of this paper, is the development of a stable preconditioner that
leads to a fast solver for the problem modeled in the parametric domain. As it
was done in [15, 16] for parabolic problems, our preconditioner exploits the
Kronecker structure of the linear system, and makes use of Fast
Diagonalization method (FD) [17]. In this work, FD is applied among the space
direction only. Although, the computational cost of the setup of the resulting
preconditioner is $O(N_{dof})$ FLOating-Point operations (FLOPs), while its
application is $O(N_{dof}^{1+1/d})$ FLOPs, where $d$ is the number of spatial
dimensions and $N_{dof}$ denotes the total number of degrees-of-freedom
(assuming, for simplicity, to have the same number of degrees-of-freedom in
time and in each spatial direction). We remark that global space-time methods,
in principle, facilitate the full parallelization of the solver, see [18, 19,
20].

The outline of the paper is as follows. The model problem is introduced in
Section 2. In Section 3 we present the basics of B-splines based IgA and the
main properties of the Kronecker product operation. The isogeometric least
squares discretization is introduced in Section 4, while in Section 5 we
define the preconditioner for the parametric domain and we discuss its
application. We present the numerical results assessing the performance of the
proposed preconditioner in Section 6. Finally, in the last section we draw
some conclusions and we highlight some future research directions.

##  2 Model problem

We consider a bounded domain ${\Omega}\subset\mathbb{R}^{d}$, usually
$d=1,2,3$, with Lipschitz boundary, and a time interval $(0,T)$, where $T>0$
is the final time. The space-time domain is denoted by
$\mathcal{Q}:=(0,T)\times{\Omega}$. Assuming Dirichlet boundary conditions,
denote by $\Gamma_{D}:=(0,T)\times\partial{\Omega}$ the Dirichlet boundary of
the space-time cylinder $\mathcal{Q}$, while
$\mathcal{Q}_{0}=\\{0\\}\times{\Omega}$ is the initial side. Our model problem
is the Schrödinger equation with homogeneous boundary and initial conditions:
we look for a solution $u$ such that

$\left\\{\begin{array}[]{rcllc}\mathrm{i}\partial_{t}u-\nu\Delta
u&=&f&\mbox{in }&\mathcal{Q},\\\\[1.0pt] u&=&0&\mbox{on
}&\Gamma_{D},\\\\[1.0pt] u&=&0&\mbox{in }&\mathcal{Q}_{0},\end{array}\right.$
(2.1)

where $\mathrm{i}$ is the imaginary unit and $\nu>0$ is a constant coefficient
usually depending on Planck’s constant $\mathfrak{\hbar}$ and the mass of the
modeled physical particle. We assume that $f\in L^{2}(\mathcal{Q})$ and denote
by $\mathbb{S}:=\mathrm{i}\partial_{t}-\nu\Delta$ the Schrödinger operator,
$\mathbb{S}^{*}$ its formal adjoint, and $(\cdot,\cdot)$ the complex scalar
product in $L^{2}(\mathcal{Q})$.

###  2.1 Space-time variational formulation

Let us introduce the Hilbert spaces

$\mathcal{V}\leavevmode\hbox{\set@color\resizebox{}{}{{\leavevmode\hbox{\set@color{$\displaystyle:=\left\\{v\in
L^{2}(\mathcal{Q}):\mathbb{S}v\in L^{2}(\mathcal{Q})\text{ and
}(\mathbb{S}^{*}w,v)-(w,\mathbb{S}v)=0\,\,\forall
w\in\mathcal{C}^{\infty}_{0}(\mathbb{R}^{d+1}):w|_{\Gamma_{D}\cup(\\{T\\}\times\Omega)}=0\right\\},$}}}}}$
$\mathcal{W}:=L^{2}(\mathcal{Q}),$

endowed with the following norms

$\|v\|_{\mathcal{V}}^{2}:=\|v\|^{2}_{L^{2}(\mathcal{Q})}+\|\mathbb{S}v\|^{2}_{L^{2}(\mathcal{Q})}\
\text{ and }\ \|w\|_{\mathcal{W}}:=\|w\|_{L^{2}(\mathcal{Q})},$

respectively. Then, the space-time variational formulation of (2.1) reads:

$\text{Find }u\in\mathcal{V}\text{ such that
}\mathcal{A}({u},v)=\mathcal{F}(v)\quad\,\forall v\in\mathcal{W},$ (2.2)

where the sesquilinear form $\mathcal{A}(\cdot,\cdot)$ and the linear form
$\mathcal{F}(\cdot)$ are defined $\forall v\in\mathcal{V}\text{ and }\forall
w\in\mathcal{W}$ as

$\mathcal{A}(v,w):=\int_{\Omega}\int_{0}^{T}\left(\mathbb{S}v\right)\overline{w}\,\mathrm{dt}\,\mathrm{d}\Omega\quad\text{and}\quad\mathcal{F}(w):=\int_{\Omega}\int_{0}^{T}f\,\overline{w}\,\mathrm{dt}\,\mathrm{d}\Omega.$

The well-posedness of the variational formulation above, for $d=1$, is in [7],
but the generalization to $d>1$ is straightforward, see the Appendix of this
paper.

###  2.2 Extensions

The previous setting can be generalized to non-homogeneous initial and
boundary conditions. For example, suppose that in (2.1) we have the initial
condition $u=u_{0}$ in $\mathcal{Q}_{0}$ with $u_{0}\in L^{2}(\Omega).$ Then,
we consider a lifting $\underline{u}_{0}$ of $u_{0}$ such that
$\underline{u}_{0}\in L^{2}(\mathcal{Q})$, see e.g. [21]. Finally, we split
the solution $u$ as $u=\underline{u}+\underline{u}_{0}$, where
$\underline{u}\in\mathcal{X}$ is the solution of the following Schrödinger
equation with homogeneous initial and boundary conditions:

$\left\\{\begin{array}[]{rcllc}\mathrm{i}\partial_{t}\underline{u}-\nu\Delta\underline{u}&=&\underline{f}&\mbox{in
}&\mathcal{Q},\\\\[1.0pt] \underline{u}&=&0&\mbox{on }&\Gamma_{D},\\\\[1.0pt]
\underline{u}&=&0&\mbox{in }&\mathcal{Q}_{0},\end{array}\right.$

where $\underline{f}:=f-\mathbb{S}\underline{u}_{0}$.

##  3 Isogeometric framework and preliminaries

###  3.1 B-Splines

Given $m$ and $p$ two positive integers, a knot vector in $[0,1]$ is a
sequence of non-decreasing points
$\Xi:=\left\\{0=\xi_{1}\leq\dots\leq\xi_{m+p+1}=1\right\\}$. We consider open
knot vectors, i.e. we set $\xi_{1}=\dots=\xi_{p+1}=0$ and
$\xi_{m}=\dots=\xi_{m+p+1}=1$, and denote by
$Z=\\{\zeta_{1},\dots,\zeta_{r}\\}$ the vector of breakpoints, that is the
vector of knots without repetition.

Univariate B-splines $\widehat{b}_{i,p}:(0,1)\rightarrow\mathbb{R}$ are
piecewise polynomials defined according to Cox-De Boor recursion formulas (see
[22]). The univariate spline space is defined as

$\widehat{\mathcal{S}}_{h}^{p}:=\mathrm{span}\\{\widehat{b}_{i,p}\\}_{i=1}^{m},$

where $h$ denotes the mesh-size, i.e. $h:=\max\\{|\xi_{i+1}-\xi_{i}|\ |\
i=1,\dots,m+p\\}$. The interior knot multiplicity influences the smoothness of
the B-splines at the knots (see [22]). For more details on B-splines
properties and their use in IgA we refer to [11].

Multivariate B-splines are defined as tensor product of univariate B-splines.
We consider functions that depend on $d$ spatial variables and the time
variable. Therefore, given positive integers $m_{l},p_{l}$ for $l=1,\dots,d$
and $m_{t},p_{t}$, we introduce $d+1$ univariate knot vectors
$\Xi_{l}:=\left\\{\xi_{l,1}\leq\dots\leq\xi_{l,m_{l}+p_{l}+1}\right\\}$ , with
associated breakpoints $Z_{l}=\\{\zeta_{l,1},\dots,\zeta_{l,r_{l}}\\}$, for
$l=1,\ldots,d$ and
$\Xi_{t}:=\left\\{\xi_{t,1}\leq\dots\leq\xi_{t,m_{t}+p_{t}+1}\right\\}$ with
$Z_{t}=\\{\zeta_{t,1},\dots,\zeta_{t,r_{t}}\\}$. Let $h_{l}$ be the mesh-size
associated to the knot vector $\Xi_{l}$ for $l=1,\dots,d$, let
$h_{s}:=\max\\{h_{l}\ |\ l=1,\dots,d\\}$ be the maximal mesh-size in all
spatial knot vectors and let $h_{t}$ be the mesh-size of the time knot vector.
Let also $\boldsymbol{p}$ be the vector that contains the degree indexes, i.e.
$\boldsymbol{p}:=(p_{t},\boldsymbol{p}_{s})$, where
$\boldsymbol{p}_{s}:=(p_{1},\dots,p_{d})$. For simplicity, we assume to have
the same polynomial degree in all spatial directions, i.e., with abuse of
notations, we set $p_{1}=\dots=p_{d}=:p_{s}$, but the general case is similar.

We assume that the following local quasi-uniformity of the knot vectors holds.

######  Assumption 1.

There exists $\theta\geq 1$, independent of $h_{s}$ and $h_{t}$, such that
$\theta^{-1}\leq\zeta_{l,i}/\zeta_{l,i+1}\leq\theta$ for $i=1,\dots,r_{l}$,
$l=1,\dots,d$ and $\theta^{-1}\leq\zeta_{t,i}/\zeta_{t,i+1}\leq\theta$ for
$i=1,\dots,r_{t}$.

The multivariate B-splines are defined as

$\widehat{B}_{\boldsymbol{i},\boldsymbol{p}}(\tau,\boldsymbol{\eta}):=\widehat{b}_{i_{t},p_{t}}(\tau)\widehat{B}_{\boldsymbol{i_{s}},\boldsymbol{p}_{s}}(\boldsymbol{\eta}),$

where

$\widehat{B}_{\boldsymbol{i_{s}},\boldsymbol{p}_{s}}(\boldsymbol{\eta}):=\widehat{b}_{i_{1},p_{s}}(\eta_{1})\ldots\widehat{b}_{i_{d},p_{s}}(\eta_{d}),$
(3.1)

$\boldsymbol{i_{s}}:=(i_{1},\dots,i_{d})$,
$\boldsymbol{i}:=(i_{t},\boldsymbol{i_{s}})$ and
$\boldsymbol{\eta}=(\eta_{1},\ldots,\eta_{d})$. The corresponding spline space
is defined as

$\widehat{\mathcal{S}}^{\boldsymbol{p}}_{{h}}:=\mathrm{span}\left\\{\widehat{B}_{\boldsymbol{i},\boldsymbol{p}}\
\middle|\ i_{t}=1,\dots,m_{t};\ i_{l}=1,\dots,m_{l}\text{ for
}l=1,\dots,d\right\\},$

where $h:=\max\\{h_{t},h_{s}\\}$. We have that
$\widehat{\mathcal{S}}^{\boldsymbol{p}}_{{h}}=\widehat{\mathcal{S}}^{p_{t}}_{h_{t}}\otimes\widehat{\mathcal{S}}^{\boldsymbol{p}_{s}}_{{h}_{s}},$
where

$\widehat{\mathcal{S}}^{\boldsymbol{p}_{s}}_{h_{s}}:=\mathrm{span}\left\\{\widehat{B}_{\boldsymbol{i_{s}},\boldsymbol{p}_{s}}\
\middle|\ i_{l}=1,\dots,m_{l};l=1,\dots,d\right\\}$

is the space of tensor-product splines on $\widehat{\Omega}:=(0,1)^{d}$.
Finally, we make the following regularity assumption.

######  Assumption 2.

We assume that $p_{t}\geq 1$, $p_{s}\geq 2$ and that
$\widehat{\mathcal{S}}^{{p}_{t}}_{h_{t}}\subset C^{0}\left((0,1)\right)$ and
$\widehat{\mathcal{S}}^{\boldsymbol{p}_{s}}_{h_{s}}\subset
C^{1}(\widehat{\Omega})$.

###  3.2 Isogeometric spaces

The space-time computational domain that we consider is $(0,T)\times\Omega$,
where $T>0$ is the final time and $\Omega\subset\mathbb{R}^{d}$ is the space
domain. The choice of considering the time as first variable will be clarified
in Section 5.2. The following assumptions asserts the regularity of the
parametrization.

######  Assumption 3.

We assume that $\Omega$ is parametrized by
$\boldsymbol{F}:\widehat{\Omega}\rightarrow{\Omega}$, with
$\boldsymbol{F}\in{\left[\widehat{\mathcal{S}}^{\boldsymbol{p}_{s}}_{{h}_{s}}\right]^{d}}$
on the closure of $\widehat{\Omega}$. Moreover, we assume that
$\boldsymbol{F}^{-1}$ has piecewise bounded derivatives of any order.

Denote by
$\boldsymbol{x}=(x_{1},\dots,x_{d}):=\boldsymbol{F}(\boldsymbol{\eta})$ and
$t:=T\tau$. Then the space-time domain is given by the parametrization
$\boldsymbol{G}:(0,1)\times\widehat{\Omega}\rightarrow(0,T)\times\Omega$, such
that
$\boldsymbol{G}(\tau,\boldsymbol{\eta}):=(T\tau,\boldsymbol{F}(\boldsymbol{\eta}))=(t,\boldsymbol{x}).$

We introduce the spline space with initial and boundary conditions, in
parametric coordinates, as

$\widehat{\mathcal{X}}_{h,0}:=\left\\{\widehat{v}_{h}\in\widehat{\mathcal{S}}^{\boldsymbol{p}}_{h}\
\middle|\ \widehat{v}_{h}=0\text{ on }\\{0\\}\times\widehat{\Omega}\text{ and
}\widehat{v}_{h}=0\text{ on }(0,1)\times\partial\widehat{\Omega}\right\\}.$

Notice that
$\widehat{\mathcal{X}}_{h,0}=\widehat{\mathcal{X}}_{t,h_{t},0}\otimes\widehat{\mathcal{X}}_{s,h_{s}}$,
where

$\displaystyle\widehat{\mathcal{X}}_{t,h_{t},0}$
$\displaystyle:=\left\\{\widehat{w}_{h}\in\widehat{\mathcal{S}}^{p_{t}}_{h_{t}}\
\middle|\ \widehat{w}_{h}(0)=0\right\\}\ =\
\text{span}\left\\{\widehat{b}_{i_{t},p_{t}}\ \middle|\ i_{t}=2,\dots,m_{t}\
\right\\},$ $\displaystyle\widehat{\mathcal{X}}_{s,h_{s}}$
$\displaystyle:=\left\\{\widehat{w}_{h}\in\widehat{\mathcal{S}}^{\boldsymbol{p}_{s}}_{h_{s}}\
\middle|\ \widehat{w}_{h}=0\text{ on }\partial\widehat{\Omega}\right\\}\ $
$\displaystyle\ =\
\text{span}\left\\{\widehat{b}_{i_{1},p_{s}}\dots\widehat{b}_{i_{d},p_{s}}\
\middle|\ i_{l}=2,\dots,m_{l}-1;\ l=1,\dots,d\ \right\\}.$

Analogously, the spline space with final and boundary conditions, in
parametric coordinates, is

$\displaystyle\widehat{\mathcal{X}}_{h,T}$
$\displaystyle:=\left\\{\widehat{v}_{h}\in\widehat{\mathcal{S}}^{\boldsymbol{p}}_{h}\
\middle|\ \widehat{v}_{h}=0\text{ on }\\{T\\}\times\widehat{\Omega}\text{ and
}\widehat{v}_{h}=0\text{ on }(0,1)\times\partial\widehat{\Omega}\right\\}$
$\displaystyle=\widehat{\mathcal{X}}_{t,h_{t},T}\otimes\widehat{\mathcal{X}}_{s,h_{s}}$

where

$\widehat{\mathcal{X}}_{t,h_{t},T}:=\left\\{\widehat{w}_{h}\in\widehat{\mathcal{S}}^{p_{t}}_{h_{t}}\
\middle|\ \widehat{w}_{h}(T)=0\right\\}\ =\
\text{span}\left\\{\widehat{b}_{i_{t},p_{t}}\ \middle|\ i_{t}=1,\dots,m_{t}-1\
\right\\}.$

More precisely, by reordering the basis functions, i.e. for the space with
initial and boundary homogeneous conditions, we write

$\displaystyle\ \widehat{\mathcal{X}}_{t,h_{t},0}$ $\displaystyle=\
\text{span}\left\\{\widehat{b}_{i,p_{t}}\ \middle|\ i=1,\dots,n_{t}\
\right\\},$ $\displaystyle\widehat{\mathcal{X}}_{s,h_{s}}$ $\displaystyle=\
\text{span}\left\\{\widehat{b}_{i_{1},p_{s}}\dots\widehat{b}_{i_{d},p_{s}}\
\middle|\ i_{l}=1,\dots,n_{s,l};\ l=1,\dots,d\ \right\\}$ $\displaystyle\
=\text{span}\left\\{\widehat{B}_{i,\boldsymbol{p}_{s}}\ \middle|\
i=1,\dots,N_{s}\ \right\\},$

and then

$\widehat{\mathcal{X}}_{h,0}=\text{span}\left\\{\widehat{B}_{{i},\boldsymbol{p}}\
\middle|\ i=1,\dots,N_{dof}\right\\},$ (3.3)

where we defined $n_{t}:=m_{t}-1$, $n_{s,l}:=m_{l}-2$ for $l=1,\dots,d$,
$N_{s}:=\prod_{l=1}^{d}n_{s,l}$ and $N_{dof}:=n_{t}N_{s}$. We can proceed
analogously with the space with final conditions.

Finally, the isogeometric space we consider is the isoparametric push-forward
of (3.3) through the geometric map $\boldsymbol{G}$, i.e.

$\mathcal{X}_{h,0}:=\text{span}\left\\{B_{i,\boldsymbol{p}}:=\widehat{B}_{i,\boldsymbol{p}}\circ\boldsymbol{G}^{-1}\
\middle|\ i=1,\dots,N_{dof}\right\\}.$ (3.4)

Notice that
$\mathcal{X}_{h,0}=\mathcal{X}_{t,h_{t},0}\otimes\mathcal{X}_{s,h_{s}}$, where

$\mathcal{X}_{t,h_{t},0}:=\text{span}\left\\{{b}_{i,p_{t}}:=\widehat{b}_{i,p_{t}}(\cdot/T)\
\middle|\ i=1,\dots,n_{t}\right\\},$

and

$\mathcal{X}_{s,h_{s}}:=\text{span}\left\\{{B}_{i,\boldsymbol{p}_{s}}:=\widehat{B}_{i,\boldsymbol{p}_{s}}\circ\boldsymbol{F}^{-1}\
\middle|\ i=1,\dots,N_{s}\right\\}.$

Analogously, we define $\mathcal{X}_{h,T}$, the isogeometric space with
homogeneous Dirichlet and final conditions.

###  3.3 Kronecker product

The Kronecker product of two matrices $\mathbf{C}\in\mathbb{C}^{n_{1}\times
n_{2}}$ and $\mathbf{D}\in\mathbb{C}^{n_{3}\times n_{4}}$ is defined as

$\mathbf{C}\otimes\mathbf{D}:=\begin{bmatrix}[\mathbf{C}]_{1,1}\mathbf{D}&\dots&[\mathbf{C}]_{1,n_{2}}\mathbf{D}\\\
\vdots&\ddots&\vdots\\\
[\mathbf{C}]_{n_{1},1}\mathbf{D}&\dots&[\mathbf{C}]_{n_{1},n_{2}}\mathbf{D}\end{bmatrix}\in\mathbb{C}^{n_{1}n_{3}\times
n_{2}n_{4}},$

where $[\mathbf{C}]_{i,j}$ denotes the $ij$-th entry of the matrix
$\mathbf{C}$. For extensions and properties of the Kronecker product we refer
to [23]. In particular, when a matrix has a Kronecker product structure, the
matrix-vector product can be efficiently computed. For this purpose, for
$m=1,\dots,d+1$ we introduce the $m$-mode product of a tensor
$\mathfrak{X}\in\mathbb{C}^{n_{1}\times\dots\times n_{d+1}}$ with a matrix
$\mathbf{J}\in\mathbb{C}^{\ell\times n_{m}}$, that we denote by
$\mathfrak{X}\times_{m}\mathbf{J}$. This is a tensor of size
$n_{1}\times\dots\times n_{m-1}\times\ell\times n_{m+1}\times\dots n_{d+1}$,
whose elements are defined as

$\left[\mathfrak{X}\times_{m}\mathbf{J}\right]_{i_{1},\dots,i_{d+1}}:=\sum_{j=1}^{n_{m}}[\mathfrak{X}]_{i_{1},\dots,i_{m-1},j,i_{m+1},\dots,i_{d+1}}[\mathbf{J}]_{i_{m},j}.$

Then, given $\mathbf{J}_{i}\in\mathbb{C}^{\ell_{i}\times n_{i}}$ for
$i=1,\dots,d+1$, it holds

$\left(\mathbf{J}_{d+1}\otimes\dots\otimes\mathbf{J}_{1}\right)\text{vec}\left(\mathfrak{X}\right)=\text{vec}\left(\mathfrak{X}\times_{1}\mathbf{J}_{1}\times_{2}\dots\times_{d+1}\mathbf{J}_{d+1}\right),$
(3.5)

where the vectorization operator “vec” applied to a tensor stacks its entries
into a column vector as

$[\text{vec}(\mathfrak{X})]_{j}=[\mathfrak{X}]_{i_{1},\dots,i_{d+1}}\text{ for
}i_{l}=1,\dots,n_{l}\text{ and for }l=1,\dots,d+1,$

where $j:=i_{1}+\sum_{k=2}^{d+1}\left[(i_{k}-1)\Pi_{l=1}^{k-1}n_{l}\right]$.

##  4 Space-time discretization of the Schrödinger equation

###  4.1 Instability of the space-time Galerkin method

Let $\mathcal{V}_{h}:=\mathcal{X}_{h,0}$ be the isogeometric space defined in
(3.4) endowed with the $\|\cdot\|_{\mathcal{V}}$-norm, and choose
$\mathcal{W}_{h}:=\mathcal{X}_{h,0}$ endowed with the
$\|\cdot\|_{\mathcal{W}}$-norm. Consider the following Galerkin method for
(2.2):

$\text{Find }u_{h}\in\mathcal{V}_{h}\text{ such that
}\mathcal{A}(u_{h},w_{h})=\mathcal{F}(w_{h})\quad\,\forall
w_{h}\in\mathcal{W}_{h}.$ (4.1)

The stability and the well-posedness of formulation (4.1) are not guaranteed,
indeed the inf-sup constant

${\alpha}_{h}:=\inf_{v_{h}\in\mathcal{V}_{h}}\sup_{w_{h}\in\mathcal{W}_{h}}\dfrac{|\mathcal{A}(v_{h},w_{h})|}{\|v_{h}\|_{\mathcal{V}}\|w_{h}\|_{\mathcal{W}}}>0,$

depends on the mesh size $h$ and degenerates under mesh refinement, as shown
in Figure LABEL:fig:galerkin_infsup_test.

(a)  (b)

Figure 1: Inf-sup test for the space-time discretizations.

###  4.2 Least squares space-time method

In order to retrieve a well posed space-time discretization to (2.2), notice
that, given the quadratic functional $\mathcal{J}:\mathcal{V}\to\mathbb{R}$,
defined as

$\mathcal{J}(u):=\dfrac{1}{2}\|\mathbb{S}u-f\|^{2}_{L^{2}(\mathcal{Q})},$

we can write the least squares space-time formulation of (2.1): find
$u\in\mathcal{V}$ such that

$u=\operatorname*{arg\,min}_{v\in\mathcal{V}}\mathcal{J}(v),$

and its Euler-Lagrange equations are

$\left(\mathbb{S}u,\mathbb{S}v\right)=\left(f,\mathbb{S}v\right),\quad\forall
v\in\mathcal{V}.$

This suggests to consider the following least-square discretization method for
(2.2). Let $\mathcal{V}_{h}:=\mathcal{X}_{h,0}$ be the isogeometric space
defined in (3.4) endowed with the $\|\cdot\|_{\mathcal{V}}$-norm, and choose
$\mathcal{W}_{h}:=\mathbb{S}(\mathcal{V}_{h})$ endowed with the
$\|\cdot\|_{\mathcal{W}}$-norm. Consider the following Petrov-Galerkin
approximation method for 2.2:

$\text{Find }u_{h}\in\mathcal{V}_{h}\text{ such that
}\mathcal{A}(u_{h},w_{h})=\mathcal{F}(w_{h})\quad\,\forall
w_{h}\in\mathcal{W}_{h}.$ (4.2)

Notice that, $\mathbb{S}$ is a bijection between the two discrete spaces,
which means, for any $h>0$ it exists the inf-sup constant
$\underline{\alpha}_{h}>0$. Moreover, for this discretization, the inf-sup
$\alpha_{h}$ is uniformly bounded from below by a positive constant
$\alpha>0$, as it is investigated numerically in Figure
LABEL:fig:least_squares_infsup_test.

######  Theorem 1.

There exists a unique solution $u_{h}\in\mathcal{V}_{h}$ to the discrete
problem (4.2). Moreover, it holds

$\|u-u_{h}\|_{\mathcal{V}}\leq\dfrac{1}{\alpha}\inf_{v_{h}\in\mathcal{V}_{h}}\|u-v_{h}\|_{\mathcal{V}},$

where $u\in\mathcal{V}$ is the solution of (2.2).

We have then the following a-priori estimate for $h$-refinement.

######  Theorem 2.

Let $q_{t}$ and $q_{s}$ be two integers, such that $q_{t}\geq 1$ and
$q_{s}\geq 2$. If $u\in\mathcal{V}\cap\left(H^{q_{t}}(0,T)\otimes
H^{2}(\Omega)\right)\cap\left(H^{1}(0,T)\otimes H^{q_{s}}(\Omega)\right)$ is
the solution of (2.2) and $u_{h}\in\mathcal{V}_{h}$ is the solution of (4.2),
then it holds

$\|u-u_{h}\|_{\mathcal{V}}\leq
C\left(h_{t}^{k_{t}-1}\|u\|_{H^{k_{t}}(0,T)\otimes
H^{2}(\Omega)}+h_{s}^{k_{s}-2}\|u\|_{H^{1}(0,T)\otimes
H^{k_{s}}(\Omega)}\right),$ (4.3)

where $k_{t}:=\min\\{q_{t},p_{t}+1\\}$, $k_{s}:=\min\\{q_{s},p_{s}+1\\}$, and
$C$ is a constant that depends only on $p_{t},p_{s},\theta$ and the
parametrization $\boldsymbol{G}$.

###### Proof.

The result follows from the anisotropic error estimates developed in [24]. We
report here only the main steps, since the proof is similar to the one of [15,
Proposition 4]. The generalization of [24, Theorem 5.1] to the $d+1$
dimensional case, gives the existence of a projection
$\Pi_{h}:\mathcal{V}\cap\left(H^{q_{t}}(0,T)\otimes
H^{2}(\Omega)\right)\cap\left(H^{1}(0,T)\otimes
H^{q_{s}}(\Omega)\right)\to\mathcal{V}_{h}$, such that

$\begin{split}\|u-\Pi_{h}u\|_{L^{2}(0,T)\otimes L^{2}(\Omega)}&\leq
C_{1}\left(h_{t}^{k_{t}-1}\|u\|_{H^{k_{t}-1}(0,T)\otimes
L^{2}(\Omega)}+h_{s}^{k_{s}-2}\|u\|_{L^{2}(0,T)\otimes
H^{k_{s}-2}(\Omega)}\right),\\\ \|u-\Pi_{h}u\|_{H^{1}(0,T)\otimes
L^{2}(\Omega)}&\leq C_{2}\left(h_{t}^{k_{t}-1}\|u\|_{H^{k_{t}}(0,T)\otimes
L^{2}(\Omega)}+h_{s}^{k_{s}-2}\|u\|_{H^{1}(0,T)\otimes
H^{k_{s}-2}(\Omega)}\right),\\\ \|u-\Pi_{h}u\|_{L^{2}(0,T)\otimes
H^{2}(\Omega)}&\leq C_{3}\left(h_{t}^{k_{t}-1}\|u\|_{H^{k_{t}-1}(0,T)\otimes
H^{2}(\Omega)}+h_{s}^{k_{s}-2}\|u\|_{L^{2}(0,T)\otimes
H^{k_{s}}(\Omega)}\right).\end{split}$ (4.4)

From the following inequality

$\displaystyle\begin{split}\|u-v_{h}\|_{\mathcal{V}}^{2}&=\|u-v_{h}\|^{2}_{L^{2}(\mathcal{Q})}+\|\mathbb{S}(u-v_{h})\|^{2}_{L^{2}(\mathcal{Q})}\\\
&\leq\|u-v_{h}\|^{2}_{L^{2}(\mathcal{Q})}+2\|\partial_{t}(u-v_{h})\|^{2}_{L^{2}(\mathcal{Q})}+2\nu\|\Delta(u-v_{h})\|^{2}_{L^{2}(\mathcal{Q})}\\\
&\leq\|u-v_{h}\|_{L^{2}(0,T)\otimes
L^{2}(\Omega)}^{2}+2\|u-v_{h}\|_{H^{1}(0,T)\otimes
L^{2}(\Omega)}^{2}+2\nu\|u-v_{h}\|_{L^{2}(0,T)\otimes
H^{2}(\Omega)}^{2},\end{split}$

with the choice $v_{h}=\Pi_{h}u$, and by (4.4) with obvious upper bounds on
the right hand side, it holds

$\|u-\Pi_{h}u\|_{\mathcal{V}}\leq
C\left(h_{t}^{k_{t}-1}\|u\|_{H^{k_{t}}(0,T)\otimes
H^{2}(\Omega)}+h_{s}^{k_{s}-2}\|u\|_{H^{1}(0,T)\otimes
H^{k_{s}}(\Omega)}\right),$

therefore, (4.3) follows from Theorem 1. ∎

######  Remark 1.

In Theorem 1, the degrees $p_{t}$, $p_{s}$ and the mesh-sizes $h_{t}$, $h_{s}$
play a similar role. This motivates our choice $p_{t}=p_{s}=:p$ and
$h_{t}=h_{s}=:h$ for the numerical tests in Section 6. In this case, and if
the solution $u$ is smooth, (4.3) yields $h$-convergence of order $p-1$.

###  4.3 Discrete system

The least-squares space-time discretization (4.2) can be written as:

$\text{Find }u_{h}\in\mathcal{V}_{h}\text{ such that
}\mathcal{A}(u_{h},\mathbb{S}v_{h})=\mathcal{F}(\mathbb{S}v_{h})\quad\,\forall
v_{h}\in\mathcal{V}_{h},$

and in particular, for all $v_{h}\in\mathcal{V}_{h}$, we point out that

$\leavevmode\hbox{\set@color\resizebox{}{}{{\leavevmode\hbox{\set@color{$\displaystyle$}}}}}\begin{split}\mathcal{A}(u_{h},\mathbb{S}v_{h})&=\int_{\Omega}\int_{0}^{T}\left(\mathbb{S}u_{h}\right)\overline{\left(\mathbb{S}v_{h}\right)}\,\mathrm{dt}\,\mathrm{d}\Omega\\\
&=\int_{\Omega}\int_{0}^{T}\partial_{t}u_{h}\overline{\partial_{t}v_{h}}+\nu^{2}\Delta
u_{h}\overline{\Delta v_{h}}+i\nu\partial_{t}\nabla u_{h}\cdot\overline{\nabla
v_{h}}-i\nu\nabla u_{h}\cdot\overline{\partial_{t}\nabla
v_{h}}\,\mathrm{dt}\,\mathrm{d}\Omega,\end{split}$

and

$\mathcal{F}(\mathbb{S}v_{h})=\int_{\Omega}\int_{0}^{T}f\,\overline{\left(\mathbb{S}v_{h}\right)}\,\mathrm{dt}\,\mathrm{d}\Omega=\int_{\Omega}\int_{0}^{T}f\,\overline{\left(i\partial_{t}v_{h}-\Delta
v_{h}\right)}\,\mathrm{dt}\,\mathrm{d}\Omega.$

Therefore, the linear system associated to (4.2) is

$\mathbf{A}\mathbf{u}=\mathbf{f},$ (4.5)

where
$[\mathbf{A}]_{i,j}=\mathcal{A}\left(B_{j,\boldsymbol{p}},\mathbb{S}(B_{i,\boldsymbol{p}})\right)$
and
$[\mathbf{f}]_{i}=\mathcal{F}\left(\mathbb{S}(B_{i,\boldsymbol{p}})\right)$.
The tensor-product structure of the isogeometric space (3.4) allows to write
the system matrix $\mathbf{A}$ as sum of Kronecker products of matrices as

$\mathbf{A}\
=\mathbf{M}_{s}\otimes\mathbf{L}_{t}+\nu^{2}\mathbf{B}_{s}\otimes\mathbf{M}_{t}+\nu\mathbf{L}_{s}\otimes\left(\mathbf{W}_{t}+\mathbf{W}_{t}^{*}\right),$
(4.6)

where for $i,j=1,\dots,N_{s}$

$\begin{split}[\mathbf{L}_{s}]_{i,j}=\int_{\Omega}\nabla
B_{i,p_{s}}(\boldsymbol{x})\cdot\nabla&B_{j,p_{s}}(\boldsymbol{x})\
\mathrm{d}\Omega,\quad[\mathbf{M}_{s}]_{i,j}=\int_{\Omega}B_{i,p_{s}}(\boldsymbol{x})\
B_{j,p_{s}}(\boldsymbol{x})\ \mathrm{d}\Omega,\\\
[\mathbf{B}_{s}]_{i,j}&=\int_{\Omega}\Delta B_{i,p_{s}}(\boldsymbol{x})\Delta
B_{j,p_{s}}(\boldsymbol{x})\ \mathrm{d}\Omega.\end{split}$ (4.7a)

while for $i,j=1,\dots,n_{t}$

$\begin{split}[\mathbf{L}_{t}]_{i,j}=\int_{0}^{T}b^{\prime}_{j,{p}_{t}}(t)\,&b^{\prime}_{i,{p}_{t}}(t)\,\mathrm{dt},\quad[\mathbf{M}_{t}]_{i,j}=\int_{0}^{T}\,b_{i,p_{t}}(t)\,b_{j,p_{t}}(t)\,\mathrm{dt},\\\
[\mathbf{W}_{t}]_{i,j}&=\mathrm{i}\int_{0}^{T}b^{\prime}_{j,{p}_{t}}(t)\,b_{i,{p}_{t}}(t)\,\mathrm{dt},\end{split}$
(4.7b)

###  4.4 Ultraweak space-time method

Here we recall also the following ultraweak discretization that has been
proposed in [9]. Let $\mathcal{W}_{h}:=\mathcal{X}_{h,T}$ be the isogeometric
space with final conditions endowed with the $\|\cdot\|_{\mathcal{W}}$-norm,
and fix $\mathcal{V}_{h}:=\mathbb{S}(\mathcal{W}_{h})$ endowed with the
$\|\cdot\|_{L^{2}(\mathcal{Q})}$-norm. Notice that, $\forall
v_{h}\in\mathcal{V}_{h},\,w_{h}\in\mathcal{W}_{h}$, it holds

$\mathcal{A}(v_{h},w_{h})=\left(\mathbb{S}(v_{h}),w_{h}\right)=\left(v_{h},\mathbb{S}(w_{h})\right)-\mathrm{i}\left(v_{h}(\cdot,0),{w}_{h}(\cdot,0)\right)_{L^{2}(\Omega)},$

with $\left(\cdot,\cdot\right)_{L^{2}(\Omega)}$ denoting the complex scalar
product in $L^{2}(\Omega)$. Therefore, introducing the sesquilinear form

$\mathcal{A}_{\mathtt{uw}}(v_{h},w_{h}):=\left(v_{h},\mathbb{S}(w_{h})\right),\quad\forall
v_{h}\in\mathcal{V}_{h},\,w_{h}\in\mathcal{W}_{h},$ (4.8)

we have the following ultraweak method for (2.2):

$\text{Find }u_{h}\in\mathcal{V}_{h}\text{ such that
}\mathcal{A}_{\mathtt{uw}}(u_{h},w_{h})=\mathcal{F}(w_{h})+\mathrm{i}\left(u_{0},{w}_{h}(\cdot,0)\right)_{L^{2}(\Omega)}\quad\,\forall
w_{h}\in\mathcal{W}_{h},$ (4.9)

where now the right hand side contains eventually the initial data $u_{0}$. As
regards the well posedness and stability of (4.9) we refer to [9].

######  Remark 2.

The discrete system associated to (4.9) has the same Kronecker structure as
(4.6), with appropriate final conditions at $T$.

##  5 Fast solver for the parametric domain

In this section we focus on the case $\boldsymbol{F}={Id}$, that is
$\mathcal{Q}=(0,T)\times(0,1)^{d}$ is the parametric domain in space and a
finite interval in time direction. In this context we are able to introduce a
stable and fast solver for problem (4.2). We introduce, for the system (4.5),
the preconditioner

$\widehat{\mathbf{P}}:=\widehat{\mathbf{M}}_{s}\otimes\mathbf{L}_{t}+\nu^{2}\widehat{\mathbf{L}}_{s}^{T}\widehat{\mathbf{M}}_{s}^{-1}\widehat{\mathbf{L}}_{s}\otimes\mathbf{M}_{t}+\nu\widehat{\mathbf{L}}_{s}\otimes(\mathbf{W}_{t}+\mathbf{W}_{t}^{*}),$
(5.1)

where the matrices $\mathbf{L}_{t},\mathbf{M}_{t}$ and $\mathbf{W}_{t}$ are
defined in (4.7b), while $\widehat{\mathbf{L}}_{s}$ and
$\widehat{\mathbf{M}}_{s}$ are

$\widehat{\mathbf{L}}_{s}=\sum_{l=1}^{d}\widehat{\mathbf{M}}_{d}\otimes\dots\otimes\widehat{\mathbf{M}}_{l+1}\otimes\widehat{\mathbf{L}}_{l}\otimes\widehat{\mathbf{M}}_{l-1}\otimes\dots\otimes\widehat{\mathbf{M}}_{1},\quad\text{and}\quad\widehat{\mathbf{M}}_{s}=\widehat{\mathbf{M}}_{d}\otimes\dots\otimes\widehat{\mathbf{M}}_{1},$

and for $l=1,\dots,d$, with indexes $i,j=1,\dots,n_{l}$, it holds

$[\widehat{\mathbf{L}}_{l}]_{i,j}:=\int\limits_{0}^{1}\widehat{b}^{\prime}_{j,p}(x_{l})\widehat{b}^{\prime}_{i,p}(x_{l})\mathrm{d}x_{l},\quad\text{
and
}\quad[\widehat{\mathbf{M}}_{l}]_{i,j}:=\int\limits_{0}^{1}\widehat{b}_{j,p}(x_{l})\widehat{b}_{i,p}(x_{l})\mathrm{d}x_{l}.$

The efficient application of the proposed preconditioner, that is, the
solution of a linear system with matrix $\widehat{\mathbf{P}}$, should exploit
the structure highlighted above. When the pencils
$({\widehat{\mathbf{L}}_{1}},\widehat{\mathbf{M}}_{1}),\ldots,({\widehat{\mathbf{L}}_{d}},\widehat{\mathbf{M}}_{d})$
admit a stable generalized eigendecomposition, a possible approach is the Fast
Diagonalization (FD) method, see [25] and [17] for details.

###  5.1 Stable factorization of
${(\widehat{\mathbf{L}}_{l},\widehat{\mathbf{M}}_{l})}$ for ${l=1,\dots,d}$

The spatial stiffness and mass matrices ${\widehat{\mathbf{L}}_{l}}$ and
$\widehat{\mathbf{M}}_{l}$ are symmetric and positive definite for
$l=1,\dots,d$. Thus, the pencils
$({\widehat{\mathbf{L}}_{l}},\widehat{\mathbf{M}}_{l})$ for $l=1,\dots,d$
admit the generalized eigendecomposition

${\widehat{\mathbf{L}}_{l}}\mathbf{U}_{l}=\widehat{\mathbf{M}}_{l}\mathbf{U}_{l}\mathbf{\Lambda}_{l},$

where the matrices $\mathbf{U}_{l}$ contain in each column the
$\widehat{\mathbf{M}}_{l}$-orthonormal generalized eigenvectors and
$\mathbf{\Lambda}_{l}$ are diagonal matrices whose entries contain the
generalized eigenvalues. Therefore we have for $l=1,\dots,d$ the
factorizations

$\mathbf{U}^{T}_{l}{\widehat{\mathbf{L}}_{l}}\mathbf{U}_{l}=\mathbf{\Lambda}_{l}\quad\text{
and
}\quad\mathbf{U}^{T}_{l}\widehat{\mathbf{M}}_{l}\mathbf{U}_{l}=\mathbb{I}_{n_{s,l}},$
(5.2)

where $\mathbb{I}_{n_{s,l}}$ denotes the identity matrix of dimension
$n_{s,l}\times n_{s,l}$. Figure 2 shows the shape of the generalized
eigenvectors in $\mathbf{U}_{l}$, with associated eigenvalue in $\Lambda_{l}$,
for a fixed univariate direction $l=1,\dots,d$ discretized with degree
$p_{s}=3$ B-Splines and uniform partition. The stability of the decomposition
is expressed by the condition number of the eigenvector matrix. In particular
$\mathbf{U}^{T}_{l}\widehat{\mathbf{M}}_{l}\mathbf{U}_{l}=\mathbb{I}_{n_{s,l}}$
implies that

$\kappa_{2}(\mathbf{U}_{l}):=\|\mathbf{U}_{l}\|_{2}\|\mathbf{U}_{l}^{-1}\|_{2}=\sqrt{\kappa_{2}(\widehat{\mathbf{M}}_{l})},$

where $\|\cdot\|_{2}$ is the norm induced by the Euclidean vector norm. The
condition number $\kappa_{2}(\widehat{\mathbf{M}}_{l})$ has been studied
theoretically in [26] and numerically in [15] and it does not depend on the
mesh-size, but it depends on the polynomial degree. Indeed, we report in Table
1 the behavior of $\kappa_{2}(\mathbf{U}_{l})$ for different values of spline
degree $p_{s}$ and for different uniform discretizations with number of
elements denoted by $n_{el}$. We observe that $\kappa_{2}(\mathbf{U}_{l})$
exhibits a dependence only on $p_{s}$, but stays moderately low for all low
polynomial degrees that are in the range of interest.

Figure 2: Generalized eigenvectors in space with associated eigenvalues for $p_{s}=3$ and $n_{el}=32$ elements. The real part is expressed in blue, while the imaginary part (null) is in red. $n_{el}$ | $p_{s}=2$ | $p_{s}=3$ | $p_{s}=4$ | $p_{s}=5$ | $p_{s}=6$ | $p_{s}=7$ | $p_{s}=8$
---|---|---|---|---|---|---|---
0032 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
0064 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
0128 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
0256 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
0512 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
1024 | $2.7\cdot 10^{0}$ | $4.5\cdot 10^{0}$ | $7.6\cdot 10^{0}$ | $1.3\cdot 10^{1}$ | $2.1\cdot 10^{1}$ | $3.5\cdot 10^{1}$ | $5.7\cdot 10^{1}$
Table 1: $\kappa_{2}(\mathbf{U}_{l})$ for different polynomial degrees $p_{s}$
and number of elements $n_{el}$.

Moreover, in [27] it is shown that there is spectral equivalence between
$\widehat{\mathbf{B}}_{s}$ and
$\widehat{\mathbf{L}}_{s}^{T}\widehat{\mathbf{M}}_{s}^{-1}\widehat{\mathbf{L}}_{s}$.
We investigate numerically this spectral equivalence, and Figure
LABEL:fig:equivalenza_in_spazio shows the eigenvalues of
$(\widehat{\mathbf{L}}_{s}^{T}\widehat{\mathbf{M}}_{s}^{-1}\widehat{\mathbf{L}}_{s})^{-1}\widehat{\mathbf{B}}_{s}$
are clustered and close to 1, for splines of degree $p=2,3,4$ and different
uniform partitionas with $n_{el}=8,16,32,64,128$. In conclusion the spectral
equivalence is stable under mesh refinement.

As regards the time pencils, the spectral equivalence between
$\widehat{\mathbf{L}}_{t}$ and
$\widehat{\mathbf{W}}_{t}^{*}\widehat{\mathbf{M}}_{t}^{-1}\widehat{\mathbf{W}}_{t}$
is unstable under mesh refinement, see Figure LABEL:fig:equivalenza_in_tempo
where we performed the analogous test, therefore we kept the full structure of
the time pencils in the preconditioner.

(a)

(b)

Figure 3: Numerical investigation of spectral equivalence.

###  5.2 Application of the preconditioner

The application of the preconditioner involves the solution of the linear
system

$\widehat{\mathbf{P}}\mathbf{s}=\mathbf{r},$ (5.3)

where $\widehat{\mathbf{P}}$ has the structure (5.1). We are able to
efficiently solve system (5.3) by the Fast Diagonalization method. The
starting point, is the setup of the preconditioner, that is the factorizations
(5.2) of the pencils $({\widehat{\mathbf{L}}_{l}},\widehat{\mathbf{M}}_{l})$
for $l=1,\dots,d$.

Then, define $\mathbf{U}_{s}:=\mathbf{U}_{d}\otimes\dots\otimes\mathbf{U}_{1}$
and
$\mathbf{\Lambda}_{s}:=\sum_{l=1}^{d}\mathbb{I}_{n_{s,d}}\otimes\dots\otimes\mathbb{I}_{n_{s,l+1}}\otimes\mathbf{\Lambda}_{l}\otimes\mathbb{I}_{n_{s,l-1}}\otimes\dots\otimes\mathbb{I}_{n_{s,1}}$.
Notice that
$\widehat{\mathbf{M}}_{s}^{-1}={\mathbf{U}}_{s}{\mathbf{U}}_{s}^{T}$,
therefore the matrix
$\widehat{\mathbf{L}}_{s}\widehat{\mathbf{M}}_{s}^{-1}\widehat{\mathbf{L}}_{s}$
admits the stable factorization

$\mathbf{U}^{T}_{s}\widehat{\mathbf{L}}_{s}\widehat{\mathbf{M}}_{s}^{-1}\widehat{\mathbf{L}}_{s}\mathbf{U}_{s}=\mathbf{\Lambda}_{s}^{2}.$

The preconditioner $\widehat{\mathbf{P}}$ admits the following factorization

$\widehat{\mathbf{P}}=\left(\mathbf{U}_{s}^{T}\otimes\mathbb{I}_{n_{t}}\right)^{-1}\left(\mathbb{I}_{N_{s}}\otimes\mathbf{L}_{t}+\nu^{2}\mathbf{\Lambda}_{s}^{2}\otimes\mathbf{M}_{t}+\nu\mathbf{\Lambda}_{s}\otimes\left(\mathbf{W}_{t}+\mathbf{W}^{*}_{t}\right)\right)\left(\mathbf{U}_{s}\otimes\mathbb{I}_{n_{t}}\right)^{-1}.$
(5.4)

Note that the second factor in (5.4) that is

$\mathbf{H}:=\left(\mathbb{I}_{N_{s}}\otimes\mathbf{L}_{t}+\nu^{2}\mathbf{\Lambda}_{s}^{2}\otimes\mathbf{M}_{t}+\nu\mathbf{\Lambda}_{s}\otimes\left(\mathbf{W}_{t}+\mathbf{W}^{*}_{t}\right)\right)$

is sum of three Kronecker matrices, whose space factors are diagonal matrices.
We have the following block diagonal structure

$\mathbf{H}=\begin{bmatrix}\mathbf{H}_{1}&&\\\ &\ddots&\\\
&&\mathbf{H}_{N_{s}}\end{bmatrix},$

where $\mathbf{H}_{i}$, for $i=1,\dots,N_{s}$, are banded matrices with
bandwidth $2p_{t}+1$ defined as

$\mathbf{H}_{i}:=\mathbf{L}_{t}+{\nu^{2}}[\mathbf{\Lambda}_{s}^{2}]_{i,i}\mathbf{M}_{n_{t}}+\nu[\mathbf{\Lambda}_{s}]_{i,i}\left(\mathbf{W}_{t}+\mathbf{W}^{*}_{t}\right).$

In order to invert $\mathbf{H}$, it is now sufficient to invert the following
independent $N_{s}$ problems of size $n_{t}\times n_{t}$:

$\mathbf{H}_{i}\mathbf{x}_{i}=\mathbf{y}_{i}\quad\text{ for }i=1,\dots,N_{s}.$
(5.5)

Summarizing, the solution of (5.3) can be computed by the following algorithm.

Algorithm 1 Fast Diagonalization

1:Compute the factorizations (5.2).

2:Compute $\mathbf{y}=(\mathbf{U}_{s}^{T}\otimes\mathbb{I}_{t})\mathbf{r}$.

3:Compute $\mathbf{x}_{i}=\mathbf{H}_{i}^{-1}\mathbf{y}_{i}\quad\text{ for
}i=1,\dots,N_{s}.$

4:Compute $\mathbf{s}=(\mathbf{U}_{s}\otimes\mathbb{I}_{t})\
\widetilde{\mathbf{s}}.$

We conclude with the following remark for a possible parallel implementation
of Algorithm 1.

######  Remark 3.

The decision to consider time as the first variable allows us to write the
matrix $\mathbf{H}$ in a block diagonal form. In view of an efficient parallel
implementation, this natural diagonal block structure does not require data
shuffling, decreasing the communication cost between nodes.

###  5.3 Computational cost and memory requirements

In this section we discuss the computational costs and memory requirements in
the implementation of Algorithm 1. First, notice that the matrix $\mathbf{A}$
in (4.6) is symmetric positive definite therefore we choose Conjugate
Gradients (CG) as linear solver for solving the system (4.5). Clearly, the
computational cost of each iteration of the CG solver depends on both the
preconditioner setup and application cost.

We assume for simplicity that, for each univariate direction $l=1,\dots,d$,
the space matrices have dimension $n_{s}\times n_{s}$, while the time matrices
involved in the preconditioners, have dimension $n_{t}\times n_{t}$. Thus the
total number of degrees-of-freedom is $N_{dof}=N_{s}n_{t}=n_{s}^{d}n_{t}$.

The setup of $\widehat{\mathbf{P}}$ includes the operations performed in Step
1 of Algorithm 1, i.e. $d$ spatial eigendecompositions, that have a total cost
of $O(dn_{s}^{3})$ FLOPs, and the construction of the block diagonal matrix
$\mathbf{H}$, which costs $O(p_{t}n_{t}N_{s})=O(p_{t}N_{dof})$. We remark that
the setup of the preconditioners has to be performed only once, since the
matrices involved do not change during the iterative procedure.

The application of the preconditioner is performed by Steps 2-4 of Algorithm
1. Exploiting (3.5), Step 2 and Step 4 costs
$O(dn_{s}^{d+1}n_{t})=O(dn_{s}N_{dof})$ FLOPs. The cost of solving each sparse
problem (5.5) makes the cost for Step 3 equal to
$O(p_{t}^{2}n_{t}N_{s})=O(p_{t}^{2}N_{dof})$ FLOPs.

In conclusion, the total cost of Algorithm 1 is
$O(dn_{s}^{3})+O(dn_{s}N_{dof})+O(p_{t}^{2}N_{dof})$ FLOPs. The non-optimal
dominant cost of Step 2 and Step 4 is determined by the dense matrix-matrix
products. However, these operations are usually implemented on modern
computers in a very efficient way and the overall serial computational time
grows almost as $O(N_{dof})$, see i.e. [15, 16]

The other dominant computational cost in a CG iteration is the cost of the
residual computation. In Algorithm 1, this involves the multiplication of the
matrix $\mathbf{A}$ with a vector. This multiplication is done by exploiting
the special structure (4.6), that allows a matrix-free approach and the use of
formula (3.5). With the matrix-free approach, noting that the time matrices
$\mathbf{L}_{t},\mathbf{M}_{t},\mathbf{W}_{t}$ are banded matrices with
bandwidth $2p_{t}+1$, and the spatial matrices
$\mathbf{J}_{s},\mathbf{L}_{s},\mathbf{M}_{s}$ have a number of non-zeros per
row equal to $(2p_{s}+1)^{d}$, the computational cost of a single matrix-
vector product is $O(N_{dof}p^{d})$ FLOPs, if we assume $p=p_{s}\approx
p_{t}$. The dominant cost in the iterative solver is represented by the
residual computation. This is a typical behaviour of the FD-based
preconditioning strategies, see [15, 28, 29].

We now investigate the memory consumption of the preconditioning strategy
proposed. For the preconditioner, we have to store the eigenvector spatial
matrices, $\mathbf{U}_{1},\dots,\mathbf{U}_{d}$, the diagonal matrices
$\mathbf{\Lambda}_{1},\dots,\mathbf{\Lambda}_{d}$ and the banded time pencils
$\mathbf{L}_{t},\mathbf{M}_{t}$ and $\mathbf{W}_{t}$ of size $n_{t}\times
n_{t}$. The memory required is roughly

$O(d(n_{s}^{2}+n_{s}))+O(p_{t}N_{dof})+O(p_{t}n_{t}).$

For the system matrix $\mathbf{A}$, in addition to the time factors
$\mathbf{L}_{t},\mathbf{M}_{t}$ and $\mathbf{W}_{t}$, we need to store the
spatial factors $\mathbf{M}_{s},\mathbf{B}_{s}$ and $\mathbf{L}_{s}$. Thus the
memory further required is roughly

$O(p_{s}^{d}N_{s}).$

These numbers show that memory-wise our space-time strategy is very appealing
when compared to other approaches, even when space and time variables are
discretized separately, e.g., with finite differences in time or other time-
stepping schemes. For example if we assume $d=3$, $p_{t}\approx p_{s}=p$ and
$n_{t}^{2}\leq Cp^{3}N_{s}$, then the total memory consumption is
$O(p^{3}N_{s}+N_{dof})$, that is equal to the sum of the memory needed to
store the Galerkin matrices associated to spatial variables and the memory
needed to store the solution of the problem.

We remark that we could avoid storing the factors of $\mathbf{A}$ by using the
matrix-free approach of [14]. The memory and the computational cost of the
iterative solver would significantly improve, both for the setup and the
matrix-vector multiplications. However, we do not pursue this strategy, as it
is beyond the scope of this paper.

##  6 Numerical Results

This section is devoted to the computation of the solution of Schrödinger
problem (2.1), and to its extension to non-homogeneous conditions, with the
discretization proposed in (4.2). We first present the numerical experiments
that assess the convergence behavior of the least squares Petrov-Galerkin
approximation and then we analyze the performance of the preconditioners.

We consider only sequential executions and we force the use of a single
computational thread in a Intel Core i5-1035G1 processor, running at 1 GHz and
with 16 GB of RAM.

The tests are performed with Matlab R2023a and GeoPDEs toolbox [30]. We use
the eig Matlab function to compute the generalized eigendecompositions present
in Step 1 of Algorithm 1,while Tensorlab toolbox [31] is employed to perform
the multiplications with Kronecker matrices occurring in Step 2 and Step 4.
The solution of the linear systems (5.5) in Step 3 is performed pagewise by
Matlab direct solver (pagewise backslash operator pagemldivide). The linear
system is solved by CG, with tolerance equal to $10^{-8}$ and with the null
vector as initial guess in all tests.

(a)  (b)  (c)

Figure 4: Smooth solution and error convergence of the space-time least
squares discretization.

According to Remark 1, we use the same mesh-size in space and in time
$h_{s}=h_{t}=:h$, and use splines of maximal continuity and same degree in
space and in time $p_{t}=p_{s}=:p$. For the sake of simplicity, we also
consider uniform knot vectors, and denote the number of elements in each
parametric direction by $n_{el}:=\frac{1}{h}$.

In out tables, the symbol $``\ast\ast"$ denotes that the invertion of the
matrix $\mathbf{A}$ in (4.6), by Matlab direct solver backslash operator
$``\backslash"$, requires more than 2 hours of computational time, while the
symbol $``\ast"$ indicates that the number of iterations in the CG solver
exceeds the upper bound set to $200$ iterates. We remark that in all the
tables the total solving time of the iterative strategies includes also the
setup time of the considered preconditioner.

###  6.1 Orders of convergence

Consider the Schrödinger equation as modeled in (2.1), for the space time
domain $\mathcal{Q}=(0,T)\times(0,1)$, with $T=2$. The reference solution is
the complex gaussian

$u(t,x)=\dfrac{\alpha\beta}{\sqrt{\beta^{2}-\mathrm{i}\gamma
t}}\exp\left\\{-\dfrac{x^{2}}{\beta^{2}-\mathrm{i}\gamma t}\right\\},$ (6.1)

where $\alpha=\beta=1.5$ and $\gamma=2.5$. Here the Dirichlet boundary
condition is $u|_{\Gamma_{D}}$, the initial condition is $u(0,x)$ and the
right hand side is $f=Au$. The problem is discretized with a uniform mesh in
both space and time directions. The solution for $p=3$ is shown in Figure
LABEL:fig:soluzione_gaussiana. In Figure LABEL:fig:convergenza_gaussiana it is
shown the convergence analysis of the error under $h$-refinement and for
different polynomial degrees, for instance $p=2,\dots,4$. The errors are
computed both with $\|\cdot\|_{\mathcal{V}}$-norm, for which the convergence
Theorem 1 holds, and with $\|\cdot\|_{L^{2}(\mathcal{Q})}$-norm, even if this
case is not covered by theoretical results. For this smooth solution, the
error study reveals optimal convergence under $h$-refinement in
$\mathcal{V}$-norm, and the $L^{2}(\mathcal{Q})$ error is lower than the
residual component in the $\mathcal{V}$-norm.

(a)  (b)

Figure 5: Non-regular solution and error convergence of the space-time
discretization.

The second test considers the following example from [7]. Consider the space
domain $\Omega=(0,1)$, the final time $T=2$, and the space-time domain
$\mathcal{Q}=(0,T)\times(0,1)$. Homogeneous Dirichlet boundary conditions are
considered on $\Gamma_{D}$. Let us denote by $e_{k}$ and $\omega^{2}_{k}$,
which is, for $k=1,2,\dots$, an eigenpair of

$-\Delta e_{k}=\omega^{2}_{k}e_{k},\quad\text{a.e. in }\Omega.$

By normalizing $e_{k}$ such that $\|e_{k}\|_{L^{2}(\Omega)}=1$, we consider
$f(t,x)=\sum_{k=1}^{+\infty}f_{k}(t)e_{k}(x),$ where $f_{k}(t)$ are the
Fourier coefficients of $f$ decomposed in the orthonormal basis $e_{k}$ at a
given time $t$. By the following specific choice of coefficients

$f_{k}(t)=\dfrac{1}{k}\exp\left\\{i\omega_{k}^{2}t\right\\}\quad\text{for
}k=1,2,\dots,$

we considered as right hand side in (2.1) the following high mode truncated
expansion

$f_{M}=\sum_{k=1}^{M}\dfrac{1}{k}\exp\left\\{i\omega_{k}^{2}t\right\\}e_{k}(x),$

with $M\gg 0$. Notice that, the solution to (2.1) with this specific right-
hand side, is

$u(t,x)=\sum_{k=1}^{M}\dfrac{-it}{k}\exp\left\\{i\omega_{k}^{2}t\right\\}e_{k}(x).$

We computed the solution for different polynomial degrees, on a uniform mesh,
for an high mode right hand side $f_{M}$, with $M=625$. In
LABEL:fig:soluzione_demkowicz it is plotted the real part of the numerical
solution for splines with degree $p=4$, together with the real part of the
explicit solution. The solution of such a problem is non-regular and it can be
shown that $u(t,\cdot)\in H^{1/2}(\Omega)$, while $u(\cdot,x)\in
H^{1/4}(0,T)$. Figure LABEL:fig:convergenza_soluzione_demkowicz shows the
error convergence for the high mode right hand side $f_{M}$, with $M=652$
modes, that is optimal for each polynomial degree $p=2,3,4$.

###  6.2 Performance of the preconditioner in the parametric domain

The computational space domain is $\Omega=(0,1)^{2}$ and the space time domain
is $\mathcal{Q}=(0,T)\times\Omega$ with $T=1$. The reference solution is a
traveling wave, that is

$u(t,\boldsymbol{x})=a\exp\left\\{-i\dfrac{|\boldsymbol{x}|^{2}+t^{2}}{\omega^{2}}\right\\},$
(6.2)

with wave number $\omega=0.2$ and amplitude $a=\sqrt[4]{2/\omega^{2}}$. Here
the Dirichlet boundary conditions are $u|_{\Gamma_{D}}$. The right hand side
is the following

$f(t,\boldsymbol{x})=a\left(\dfrac{4i}{\omega^{2}}+\dfrac{4|\boldsymbol{x}|^{2}}{\omega^{4}}+\dfrac{2t}{\omega^{2}}\right)\exp\left\\{-i\dfrac{|\boldsymbol{x}|^{2}+t^{2}}{\omega^{2}}\right\\}.$

The numerical solution on a mesh of $64$ elements per univariate direction is
shown in Figure 6 for different time frames.

(a)

(b)

(c)

(d)

Figure 6: Real part of the numerical approximation of (6.2) at four different
time frames. Table 2: Parametric domain. Performance of
$\widehat{\mathbf{P}}$. Performance of preconditioner
---
$N_{e}$ | $N_{dof}$ |  backslash (time) |  $\widehat{\mathbf{P}}$ (iter / time) | ICHOL (iter / time)
Degree $p=2$
8 | 1000 | 0.0306 | 7 / 0.0426 | 12 / 0.0256
16 | 5832 | 0.7974 | 7 / 0.0659 | 38 / 0.1338
32 | 39304 | 27.46 | 7 / 0.2992 | 174 / 3.9153
64 | 287496 | 1148 | 7 / 2.7304 | ∗
Degree $p=3$
8 | 1331 | 0.0618 | 8 / 0.0302 | 10 / 0.0172
16 | 6859 | 1.6415 | 9 / 0.0982 | 32 / 0.2884
32 | 42875 | 118 | 8 / 0.4703 | 148 / 8.1482
64 | 300763 | ∗∗ | 8 / 4.7142 | ∗
Degree $p=4$
8 | 1728 | 0.1587 | 10 / 0.0427 | 10 / 0.0548
16 | 8000 | 2.6433 | 10 / 0.2289 | 28 / 0.4363
32 | 46656 | 419 | 10 / 1.2329 | 127 / 14.1742
64 | 314432 | ∗∗ | 10 / 9.0219 | 191 / 176

We analyze the performance of the proposed preconditioner
$\widehat{\mathbf{P}}$ for a variety of uniform partitions up to $n_{el}=64$
per univariate directions, and polynomial degree $p=2,3,4$. In Table 2 it is
reported the computational clock time cost of solving the linear system both
directly using Matlab backslash and iteratively by CG solver, reporting also
the number of iterations for this latter case. When solving with CG we
investigate the performance of the solver with preconditioner
$\widehat{\mathbf{P}}$, and compare it with a classical algebraic
preconditioner as incomplete Cholesky factorization (ICHOL). Matlab backslash
is clearly inefficient, since performing gaussian elimination requires
$N_{dof}^{3}$ FLOPs. Using classical preconditioners in CG iterative solvers
is a reasonable approach for small size problem, but the number of iterations
grows with the size of the problem. The performance of the preconditioner
$\widehat{\mathbf{P}}$ with CG, is identical among $h$-refined meshes, and
seams reasonably $p$-robust. The number of iterations never exceeds $10$, and
the total amount of time required to solve the discrete problem, is always
cheaper than the other approaches we tested.

##  7 Conclusions

In this work we proposed and studied a space-time least square method for the
Schrödinger equation in the framework of isogeometric analysis. Our scheme is
based on smooth spline in space and time, that allows, in the particular case
of the parametric domain, to introduce a suitable preconditioner for the
arising linear system. Our preconditioner $\widehat{\mathbf{P}}$ is
represented by a sum of Kronecker products of matrices, that makes the
computational cost of its construction (setup) and application, as well as the
storage cost, very appealing. In particular the construction of the
preconditioner exploits a spectral equivalence between the space matrices
$\mathbb{B}_{s}$ and $\mathbf{L}_{s}^{T}\mathbf{M}_{s}^{-1}\mathbf{L}_{s}$
that , thanks to the FD technique, admits a stable block-diagonal
factorization.

The application cost for a serial execution is almost equal to $O(N_{dof})$,
and the block-diagonal structure is suitable for parallel implementation on
distributed memory machines, and this will be an interesting future direction
of study.

At the same time, the storage cost is roughly the same that we would have by
discretizing separately in space and in time, if we assume $n_{t}\leq
Cp^{d}N_{s}$. Indeed, in this case the memory used for the whole iterative
solver is $O(p^{d}N_{s}+N_{dof})$. Although, our approach could be coupled
with a matrix-free idea, and this is expected to further improve the
efficiency of the overall method.

As a final comment, it would be interesting to further exploits the structure
of time pencils, in order to achive a full factorization of the proposed
preconditioner. This may also give a hint in proposing an ad-hoc
preconditioner for the isogeometric framework, which we are still working on.

##  Appendix A Well-posedness of the space-time variational formulation

Here we extend the results presented in [7] on the well posedness of (2.2).
First we introduce a suitable notation, such that this appendix can be read
independently from the paper. Let us recall $\mathcal{Q}=(0,T)\times\Omega$,
with $\Omega\subset\mathbb{R}^{d}$, $d=1,2,3$, while
$\Gamma_{D}=(0,T)\times\partial\Omega$. Consider
$\Gamma_{0}=\Gamma_{D}\cup(\\{0\\}\times\Omega)$ and
$\Gamma_{T}=\Gamma_{D}\cup(\\{T\\}\times\Omega)$ and let us define
$\mathcal{D}_{0}:=\left\\{\phi\in
C^{\infty}_{0}(\mathbb{R}^{d+1}):\phi|_{\Gamma_{0}}=0\right\\}$, which is the
space of smooth functions of $\mathbb{R}^{d+1}$ with compact support such that
restricted to $\mathcal{Q}$ satisfy both homogeneous Dirichlet and initial
conditions. Analogously define $\mathcal{D}_{T}:=\left\\{\psi\in
C^{\infty}_{0}(\mathbb{R}^{d+1}):\psi|_{\Gamma_{T}}=0\right\\}$, that instead
satisfies homogeneous Dirichlet and final conditions. Recall
$\mathbb{S}:=\mathrm{i}\partial_{t}-\nu\Delta$, and notice that integration by
parts gives:

$\int_{\Omega}\int_{0}^{T}\left(\mathbb{S}\phi\right)\overline{\psi}\,\mathrm{dt}\,\mathrm{d}\Omega=\int_{\Omega}\int_{0}^{T}\phi\left(\overline{\mathbb{S}\psi}\right)\,\mathrm{dt}\,\mathrm{d}\Omega,\quad\forall\phi\in\mathcal{D}_{0},\text{
and }\forall\psi\in\mathcal{D}_{T}.$

The space $\mathcal{V}$ in (2.2) is the domain of
$\mathbb{S}:\mathcal{V}\subset L^{2}(\mathcal{Q})\to L^{2}(\mathcal{Q})$, that
can be written as:

$\mathcal{V}:=\left\\{v\in L^{2}(\mathcal{Q}):\mathbb{S}v\in
L^{2}(\mathcal{Q})\text{ and
}(\mathbb{S}\psi,v)-(\psi,\mathbb{S}v)=0\,\,\forall\psi\in\mathcal{D}_{T}\right\\},$
(A.1)

and we have $\mathcal{C}^{\infty}_{0}(\mathcal{Q})\subset\mathcal{V}\subset
L^{2}(\mathcal{Q})$, that is $\mathbb{S}$ is densely defined. Denoting by
$\mathbb{S}^{*}:\mathcal{V}^{*}\subset L^{2}(\mathcal{Q})\to
L^{2}(\mathcal{Q})$ the adjoint operator, whose domain is given by

$\mathcal{V}^{*}:=\left\\{w\in L^{2}(\mathcal{Q}):\exists g\in
L^{2}(\mathcal{Q})\text{ such that }(\mathbb{S}v,w)=(v,g)\,\,\forall
v\in\mathcal{V}\right\\},$ (A.2)

we have $\mathbb{S}^{*}w:=g$, with $\mathbb{S}^{*}=\mathbb{S}$, and
$\mathcal{C}^{\infty}_{0}(\mathcal{Q})\subset\mathcal{V}^{*}\subset
L^{2}(\mathcal{Q})$. Notice that we are identifying
$L^{2}(\mathcal{Q})^{\prime}\equiv L^{2}(\mathcal{Q})$ through Riesz
isomorphism. We endow both $\mathcal{V}$ and $\mathcal{V}^{*}$ with the norms
$\|\cdot\|_{\mathcal{V}}$ and $\|\cdot\|_{\mathcal{V}^{*}}$ respectively, such
that

$\|v\|_{\mathcal{V}}^{2}:=\|v\|_{L^{2}(\mathcal{Q})}^{2}+\|\mathbb{S}v\|_{L^{2}(\mathcal{Q})}^{2},\quad\text{and}\quad\|v\|_{\mathcal{V}^{*}}^{2}:=\|v\|_{L^{2}(\mathcal{Q})}^{2}+\|\mathbb{S}^{*}v\|_{L^{2}(\mathcal{Q})}^{2}.$

Define the boundary operators $B:\mathcal{V}\to(\mathcal{V}^{*})^{\prime}$ and
$B^{*}:\mathcal{V}^{*}\to(\mathcal{V})^{\prime}$, such that

$\langle
Bv,w\rangle:=(\mathbb{S}v,w)_{L^{2}(\mathcal{Q})}-(v,\mathbb{S}^{*}w)_{L^{2}(\mathcal{Q})},$
(A.3a) $\langle
B^{*}w,v\rangle:=(\mathbb{S}^{*}w,v)_{L^{2}(\mathcal{Q})}-(w,\mathbb{S}v)_{L^{2}(\mathcal{Q})},$
(A.3b)

hold true for all $v,w\in L^{2}({\mathcal{Q}})$ such that
$\mathbb{S}v,\mathbb{S}^{*}w\in L^{2}(\mathcal{Q})$. From [7, Lemma A.2], we
have

$\mathcal{V}^{*}=(B(\mathcal{V}))^{\bot},$ (A.4)

and, from (A.1), it holds

$\mathcal{V}=(B^{*}(\mathcal{D}_{T}))^{\bot}.$ (A.5)

In particular, from [7, Lemma 2.1], it holds
$\mathcal{D}_{0}\subset\mathcal{V}$ and
$\mathcal{D}_{T}\subset\mathcal{V}^{*}$, and in addition we make the following
density assumption.

######  Assumption 4.

We assume that
$\overline{\mathcal{D}_{0}}^{\|\cdot\|_{\mathcal{V}}}=\mathcal{V}$ and
$\overline{\mathcal{D}_{T}}^{\|\cdot\|_{\mathcal{V}^{*}}}=\mathcal{V}^{*}$.

Under Assumption 4, from [7, Lemma 2.2], it holds

$\mathcal{V}^{*}=(B(\mathcal{V}))^{\bot}=(B(\mathcal{D}_{0}))^{\bot},$ (A.6)

and

$\mathcal{V}=(B^{*}(\mathcal{V}^{*}))^{\bot}=(B^{*}(\mathcal{D}_{T}))^{\bot}.$
(A.7)

The next result states the well posedness of the space time variational
formulation.

######  Theorem 3.

Under Assumption 4, the linear Schrödinger operator $\mathbb{S}:\mathcal{V}\to
L^{2}(\mathcal{Q})$ is a continuous bijections, that is problem (2.2) is well
posed.

The proof of well posedness is given in [7, Theorem 2.4], and Assumption 4,
but more precisely (A.6), is needed only to prove injectivity of $\mathbb{S}$.
We now prove that, Assumption (4) is verified for every integer $d\geq 1$.

######  Lemma 1.

Given $\mathcal{Q}=(0,T)\times\Omega$, with $\Omega=[0,1]^{d}$ and integer
$d\geq 1$, then Assumption (4) holds true.

###### Proof.

We prove that $\mathcal{D}_{0}$ is dense in $\mathcal{V}$, the other stated
density result is analogous. The case $d=1$ is in [7, Theorem 3.1], thus we
fix an integer $d>1$. Consider $v\in\mathcal{V}$, first we extend $v$ to the
whole $\mathbb{R}^{d+1}$ domain.

  1. 1.

Extending along space directions: Let us extend the space-time domain among
the space directions as follows. Denote by $\mathcal{Q}_{1,c}=\mathcal{Q}$,
$\mathcal{Q}_{1,l}:=[0,T]\times[-1,0]\times[0,1]^{d-1}$ and
$\mathcal{Q}_{1,r}:=[0,T]\times[1,2]\times[0,1]^{d-1}$. Analogously, for
$i=2,\dots,d$, we introduce
$\mathcal{Q}_{i,c}:=\bigcup_{j\in\\{l,c,r\\}}\mathcal{Q}_{i-1,j}$, then
$\mathcal{Q}_{i,l}:=[0,T]\times[-1,2]^{i-1}\times[-1,0]\times[0,1]^{d-i}$ and
$\mathcal{Q}_{i,r}:=[0,T]\times[-1,2]^{i-1}\times[1,2]\times[0,1]^{d-i}$,
considering $[0,1]^{0}=\emptyset$. Finally let us call
$\Omega_{E}:=\bigcup_{j\in\\{l,c,r\\}}\mathcal{Q}_{d,j}$ the enlarged space-
time cylinder. Then, we introduce the intermediate extension operators,
$E_{i}:\mathcal{Q}_{i,c}\to\mathcal{Q}_{i+1,c}$ for $i=1,\dots,d-1$, and
$E_{d}:\mathcal{Q}_{d,c}\to\mathcal{Q}_{E}$, such that

$E_{i}f(t,\boldsymbol{x}):=\begin{cases}-f(t,x_{1},\dots,-x_{i},\dots,x_{d})&(t,\boldsymbol{x})\in\mathcal{Q}_{i,l},\\\
f(t,x_{1},\dots,x_{i},\dots,x_{d})&(t,\boldsymbol{x})\in\mathcal{Q}_{i,c},\\\
-f(t,x_{1},\dots,2-x_{i},\dots,x_{d})&(t,\boldsymbol{x})\in\mathcal{Q}_{i,r}.\end{cases}$

We denote the reverse operators by $E_{i}^{\prime}$, defined as
$E_{i}^{\prime}g(t,\boldsymbol{x})=g(t,\boldsymbol{x})-g(t,x_{1},\dots,-x_{i},\dots,x_{d})-g(t,x_{1},\dots,2-x_{i},\dots,x_{d})$,
for $(t,\boldsymbol{x})\in\mathcal{Q}_{i,c}$, and for $i=1,\dots,d$. The
definitions are to be interpreted almost everywere, and finally our extension
operator from $\mathcal{Q}$ to $\mathcal{Q}_{E}$ is $E:=E_{d}\circ\dots,\circ
E_{1}$, while its reverse operator from $\mathcal{Q}_{E}$ to $\mathcal{Q}$ is
$E^{\prime}=E_{1}^{\prime}\circ\dots\circ E_{d}^{\prime}$. It is easy to see
by a change of variable that

$(Ef,g)_{L^{2}(\mathcal{Q}_{E})}=(f,E^{\prime}g)_{L^{2}(\mathcal{Q})},\quad\forall
f\in L^{2}(\mathcal{Q}),\forall g\in L^{2}(\mathcal{Q}_{E}).$

Next, we claim that

$\mathbb{S}Ev=E\mathbb{S}v,\quad\forall v\in\mathcal{V}.$

Clearly, $Ev$ is in $L^{2}(\mathcal{Q}_{E})$ and notice that
$E^{\prime}\mathbb{S}\varphi=\mathbb{S}E^{\prime}\varphi$ for all
$\varphi\in\mathcal{C}^{\infty}_{0}(\mathcal{Q}_{E})$. Therefore, it holds

$\begin{split}\langle\mathbb{S}Ev,\varphi\rangle_{\mathcal{C}^{\infty}_{0}(\mathcal{Q}_{E})}&=(Ev,\mathbb{S}\varphi)_{L^{2}(\mathcal{Q}_{E})}=(v,E^{\prime}\mathbb{S}\varphi)_{L^{2}(\mathcal{Q})}=\\\
&=(v,\mathbb{S}E^{\prime}\varphi)_{L^{2}(\mathcal{Q})}=(\mathbb{S}v,E^{\prime}\varphi)_{L^{2}(\mathcal{Q})}-\langle
Bv,E^{\prime}\varphi\rangle.\end{split}$

Now, since $E^{\prime}\varphi|_{\Gamma_{T}}=0$, we have
$E^{\prime}\varphi\in\mathcal{V}^{*}$, and thus $\langle
Bv,E^{\prime}\varphi\rangle=0$ by (A.4). It follows that

$\langle\mathbb{S}Ev,\varphi\rangle_{\mathcal{C}^{\infty}_{0}(\mathcal{Q}_{E})}=(\mathbb{S}v,E^{\prime}\varphi)_{L^{2}(\mathcal{Q})}=(E\mathbb{S}v,E\varphi)_{L^{2}(\mathcal{Q}_{E})},$

completing the proof of the claim. We also conclude that $\mathbb{S}Ev$ is in
$L^{2}(\mathcal{Q}_{E})$ whenever $v\in\mathcal{V}$.

  2. 2.

Extending along time direction: Let $\widetilde{E}$ denote the extension of
$E$ by zero to $\mathbb{R}^{d+1}$, and $\tau_{\delta}$ be the translation
operator in $t$ direction, i.e.,
$\tau_{\delta}w(t,\boldsymbol{x})=w(t-\delta,\boldsymbol{x})$. From [32] it
holds

$\lim_{\delta\to
0}\|\tau_{\delta}w-w\|_{L^{2}(\mathbb{R}^{d+1})}=0,\quad\forall w\in
L^{2}({\mathbb{R}^{d+1}}).$ (A.8)

Introducing $\mathcal{Q}_{E,\delta}=(-\delta,T+\delta)\times(-1,2)^{d}$, by a
change of variables, it holds

$(\tau_{\delta}\widetilde{E}f,g)_{L^{2}(\mathcal{Q}_{E,\delta})}=(Ef,\tau_{-\delta}g)_{L^{2}(\mathcal{Q}_{E})},\quad\forall
f\in L^{2}(\mathcal{Q}),\forall g\in L^{2}(\mathcal{Q}_{E,\delta}).$

Denoting by $R_{\delta}$ the restriction operator of function on
$\mathbb{R}^{d+1}$ to $\mathcal{Q}_{E,\delta}$, we now claim that

$\mathbb{S}R_{\delta}\tau_{\delta}\widetilde{E}v=R_{\delta}\tau_{\delta}\widetilde{E}\mathbb{S}v,\quad\forall
v\in\mathcal{V}.$

The proof is analogous to the one in Step 1. Given
$\varphi\in\mathcal{C}^{\infty}_{0}(\mathcal{Q}_{E,\delta})$ it holds

$\begin{split}\langle\mathbb{S}R_{\delta}\tau_{\delta}\widetilde{E}v,\varphi\rangle_{\mathcal{C}^{\infty}_{0}(\mathcal{Q}_{E,\delta})}&=(\tau_{\delta}\widetilde{E}v,\mathbb{S}\varphi)_{L^{2}(\mathcal{Q}_{E,\delta})}=(Ev,\mathbb{S}\tau_{-\delta}\varphi)_{L^{2}(\mathcal{Q}_{E})}\\\
&=(v,E^{\prime}\mathbb{S}\tau_{-\delta}\varphi)_{L^{2}(\mathcal{Q})}=(v,\mathbb{S}E^{\prime}\tau_{-\delta}\varphi)_{L^{2}(\mathcal{Q})}\\\
&=(\mathbb{S}v,E^{\prime}\tau_{-\delta}\varphi)_{L^{2}(\mathcal{Q})}-\langle
Bv,E^{\prime}\tau_{-\delta}\varphi\rangle.\end{split}$

Now, since $E^{\prime}\tau_{-\delta}\varphi|_{\Gamma_{T}}=0$, we have
$E^{\prime}\tau_{-\delta}\varphi\in\mathcal{V}^{*}$, and thus $\langle
Bv,E^{\prime}\tau_{-\delta}\varphi\rangle=0$ by (A.4). It follows that

$\langle\mathbb{S}R_{\delta}\tau_{\delta}\widetilde{E}v,\varphi\rangle_{\mathcal{C}^{\infty}_{0}(\mathcal{Q}_{E,\delta})}=(\widetilde{E}\mathbb{S}v,\tau_{-\delta}\varphi)_{L^{2}(\mathcal{Q}_{E})}=(\tau_{-\delta}\widetilde{E}\mathbb{S}v,\varphi)_{L^{2}(\mathcal{Q}_{E,\delta})},$

which proves the claim.

  3. 3.

Mollify: Consider the mollifier
$\rho_{\epsilon}\in\mathcal{C}^{\infty}_{0}(\mathbb{R}^{d+1})$, defined by

$\rho_{\epsilon}(t,\boldsymbol{x}):=\epsilon^{-d-1}\rho_{1}(\epsilon^{-1}t,\epsilon^{-1}x_{1},\dots,\epsilon^{-1}x_{d}),\quad\text{for}\epsilon>0,$

where

$\rho_{1}(t,\boldsymbol{x}):=\begin{cases}ke^{-1/{(1-|(t,\boldsymbol{x})|^{2})}},&\text{if
}|(t,\boldsymbol{x})|^{2}<1,\\\ 0&\text{if }|(t,\boldsymbol{x})|^{2}\geq
1,\end{cases}$

with $|\cdot|$ denoting the Euclidean norm in $\mathbb{R}^{d+1}$, and $k$ is a
constant chosen such that $\int_{\mathbb{R}^{d+1}}\rho_{1}=1$. Notice that,
given $\delta>0$ small enough, i.e., $\delta<\min\\{T/2,1/2\\}$, the
convolutions $v_{\epsilon}:=\rho_{\epsilon}\ast\tau_{\delta}\widetilde{E}v$
and $s_{\epsilon}:=\rho_{\epsilon}\ast\tau_{\delta}\widetilde{E}\mathbb{S}v$
are smooth functions that satisfy

$\lim_{\epsilon\to
0}\|v_{\epsilon}-\tau_{\delta}\widetilde{E}v\|_{L^{2}(\mathbb{R}^{d+1})}=0,\quad\text{and}\quad\lim_{\epsilon\to
0}\|s_{\epsilon}-\tau_{\delta}\widetilde{E}\mathbb{S}v\|_{L^{2}(\mathbb{R}^{d+1})}=0.$
(A.9)

Moreover, the smooth function $\mathbb{S}v_{\epsilon}$ need not coincide to
$s_{\epsilon}$ everywere, but they coincide on $\mathcal{Q}$ whenever
$\epsilon<\delta/2$. Thus, consider $\delta=3\epsilon$, and let
$\epsilon<\min\\{T/6,1/6\\}$ go to zero. We have

$\|\mathbb{S}v_{\epsilon}-\mathbb{S}v\|_{L^{2}(\mathcal{Q})}=\|s_{\epsilon}-\mathbb{S}v\|_{L^{2}(\mathcal{Q})}\leq\|s_{\epsilon}-\tau_{\delta}\widetilde{E}\mathbb{S}v\|_{L^{2}(\mathbb{R}^{d+1})}+\|\tau_{\delta}\widetilde{E}\mathbb{S}v-\widetilde{E}\mathbb{S}v\|_{L^{2}(\mathbb{R}^{d+1})}$
(A.10) and
$\|v_{\epsilon}-v\|_{L^{2}(\mathcal{Q})}\leq\|v_{\epsilon}-\tau_{\delta}\widetilde{E}v\|_{L^{2}(\mathcal{Q})}+\|\tau_{\delta}\widetilde{E}v-\widetilde{E}v\|_{L^{2}(\mathcal{Q})}$

Using (A.8) and (A.9), it follows that

$\lim_{\epsilon\to 0}\|v_{\epsilon}-v\|_{\mathcal{V}}=0.$

To conclude, we examine the value of $v_{\epsilon}$ at the edges of the space-
time cylinder. We have

$\displaystyle
v_{\epsilon}(t,0,x_{2},\dots,x_{d})=\int_{\mathbb{R}^{d+1}}\rho_{\epsilon}(t-\sigma,-r_{1},x_{2}-r_{2},\dots,x_{d}-r_{d})\tau_{\delta}\widetilde{E}v(\sigma,r_{1},\dots,r_{d})\mathrm{d}r_{1}\dots\mathrm{d}r_{d}\mathrm{d}\sigma,$

with the integrand the inner integral being the product of and even function
$\rho_{\epsilon}$, with respect to $r_{1}$, and an odd function
$\tau_{\delta}\widetilde{E}v$ of $r_{1}$. Thus
$v_{\epsilon}(t,0,x_{2},\dots,x_{d})=0$ and the same holds for
$v_{\epsilon}(t,1,x_{2},\dots,x_{d})$ and the other univariate space
dierctions. Moreover since $\tau_{\delta}\widetilde{E}v$ is identically zero
in a neighborhood of $(0,\boldsymbol{x})$, we conclude that
$v_{\epsilon}|_{\Gamma_{0}}=0$.

∎

Next we extend this result to the isogeometric case domain.

######  Theorem 4.

Given $\mathcal{Q}=(0,T)\times\Omega$, with $\Omega=\mathbf{F}([0,1]^{d})$ and
integer $d\geq 1$, then Assumption (4) holds true.

###### Proof.

We prove that $\mathcal{D}_{0}$ is dense in $\mathcal{V}$, the other stated
density result is analogous. Given $v\in\mathcal{V}$, recall
$\mathbf{G}:[0,1]^{d+1}\to\mathcal{Q}$ is the parameterization of the space-
time cilinder, such that
$\mathbf{G}(\tau,\boldsymbol{\eta}):=(T\tau,\mathbf{F}(\boldsymbol{\eta}))=(t,\boldsymbol{x})$.
Define $\widehat{v}:=v\circ\mathbf{G}$. Clearly $\widehat{v}\in
L^{2}([0,1]^{d+1})$, and $\mathbb{S}\widehat{v}\in L^{2}([0,1]^{d+1})$.
Moreover
$(\mathbb{S}v,\phi)_{L^{2}([0,1]^{d+1})}-(v,\mathbb{S}\phi)_{L^{2}([0,1]^{d+1})}=0$,
for all $\phi\in\mathcal{C}^{\infty}_{0}(\mathbb{R}^{d+1})$ such that
$\phi|_{\mathbf{G}^{-1}(\Gamma_{T})}=0$. By applying Lemma 1, it exists
$\\{\widehat{v}_{\epsilon}\\}_{\epsilon>0}\subset\mathcal{C}^{\infty}_{0}(\mathbb{R}^{d+1})$
such that $\widehat{v}_{\epsilon}|_{\mathbf{G}^{-1}(\Gamma_{0})}=0$, that
satisfies

$\lim_{\epsilon\to
0}\|\widehat{v}_{\epsilon}-\widehat{v}\|^{2}_{L^{2}([0,1]^{d+1})}+\|\mathbb{S}\widehat{v}_{\epsilon}-\mathbb{S}\widehat{v}\|^{2}_{L^{2}([0,1]^{d+1})}=0.$

Therefore define $v_{\epsilon}:=\widehat{v}_{\epsilon}\circ\mathbf{G}^{-1}$
and notice that $\\{v_{\epsilon}\\}_{\epsilon>0}\subset\mathcal{D}_{0}$.
Moreover, by a change of variable,

$\lim_{\epsilon\to 0}\|{v}_{\epsilon}-{v}\|^{2}_{\mathcal{V}}\leq
C\lim_{\epsilon\to
0}\|\widehat{v}_{\epsilon}-\widehat{v}\|^{2}_{L^{2}([0,1]^{d+1})}+\|\mathbb{S}\widehat{v}_{\epsilon}-\mathbb{S}\widehat{v}\|^{2}_{L^{2}([0,1]^{d+1})}=0.$

This completes the proof. ∎

## Acknowledgments

## References

  * [1] I. Fried, Finite-element analysis of time-dependent phenomena., AIAA Journal 7 (6) (1969) 1170–1173.
  * [2] J. C. Bruch Jr., G. Zyvoloski, Transient two-dimensional heat conduction problems solved by the finite element method, International Journal for Numerical Methods in Engineering 8 (3) (1974) 481–494.
  * [3] J. T. Oden, A general theory of finite elements. I. Topological considerations, International Journal for Numerical Methods in Engineering 1 (2) (1969) 205–221.
  * [4] J. T. Oden, A general theory of finite elements. II. Applications, International Journal for Numerical Methods in Engineering 1 (3) (1969) 247–259.
  * [5] F. Shakib, T. J. R. Hughes, A new finite element formulation for computational fluid dynamics: IX. Fourier analysis of space-time Galerkin/least-squares algorithms, Computer Methods in Applied Mechanics and Engineering 87 (1) (1991) 35–58.
  * [6] O. Karakashian, C. Makridakis, A space-time finite element method for the nonlinear Schrödinger equation: the discontinuous Galerkin method, Mathematics of computation 67 (222) (1998) 479–499.
  * [7] L. Demkowicz, J. Gopalakrishnan, S. Nagaraj, P. Sepulveda, A spacetime DPG method for the Schrödinger equation, SIAM Journal on Numerical Analysis 55 (4) (2017) 1740–1759.
  * [8] S. Gómez, A. Moiola, A space-time Trefftz discontinuous Galerkin method for the linear Schrödinger equation, SIAM Journal on Numerical Analysis 60 (2) (2022) 688–714.
  * [9] S. Hain, K. Urban, An ultra-weak space-time variational formulation for the Schrödinger equation, arXiv preprint arXiv:2212.14398 (2022).
  * [10] T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering 194 (39) (2005) 4135–4195.
  * [11] J. A. Cottrell, T. J. R. Hughes, Y. Bazilevs, Isogeometric analysis: toward integration of CAD and FEA, John Wiley & Sons, 2009.
  * [12] J. A. Evans, Y. Bazilevs, I. Babuška, T. J. R. Hughes, $n$-widths, sup-infs, and optimality ratios for the $k$-version of the isogeometic finite element method, Computer Methods in Applied Mechanics and Engineering 198 (2009) 1726–1741.
  * [13] A. Bressan, E. Sande, Approximation in FEM, DG and IGA: a theoretical comparison, Numerische Mathematik (2019).
  * [14] G. Sangalli, M. Tani, Matrix-free weighted quadrature for a computationally efficient isogeometric $k$-method, Computer Methods in Applied Mechanics and Engineering 338 (2018) 117 – 133.
  * [15] M. Montardini, M. Negri, G. Sangalli, M. Tani, Space-time least-squares isogeometric method and efficient solver for parabolic problems, Mathematics of Computation (accepted for publication) (2019).
  * [16] G. Loli, M. Montardini, G. Sangalli, M. Tani, An efficient solver for space–time isogeometric Galerkin methods for parabolic problems, Computers & Mathematics with Applications 80 (11) (2020) 2586–2603.
  * [17] R. E. Lynch, J. R. Rice, D. H. Thomas, Direct solution of partial difference equations by tensor product methods, Numerische Mathematik 6 (1) (1964) 185–199.
  * [18] C. A. Dorao, H. A. Jakobsen, A parallel time–space least-squares spectral element solver for incompressible flow problems, Applied Mathematics and Computation 185 (1) (2007) 45–58.
  * [19] M. J. Gander, 50 years of time parallel time integration, in: Multiple Shooting and Time Domain Decomposition Methods, Springer, 2015, pp. 69–113.
  * [20] A. M. Kvarving, E. M. Rønquist, A fast tensor-product solver for incompressible fluid flow in partially deformed three-dimensional domains: Parallel implementation, Computers & Fluids 52 (2011) 22–32.
  * [21] L. C. Evans, Partial Differential equations, American Mathematical Society, Berlin, 2010.
  * [22] C. De Boor, A practical guide to splines (revised edition), Applied Mathematical Sciences, Springer, Berlin, 2001.
  * [23] T. G. Kolda, B. W. Bader, Tensor decompositions and applications, SIAM review 51 (3) (2009) 455–500.
  * [24] L. Beirão da Veiga, D. Cho, G. Sangalli, Anisotropic NURBS approximation in isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 209 (2012) 1–11.
  * [25] M. O. Deville, P. F. Fischer, E. H. Mund, High-order methods for incompressible fluid flow, Cambridge University Press, 2002.
  * [26] K. P. S. Gahalaut, S. K. Tomar, C. Douglas, Condition number estimates for matrices arising in NURBS based isogeometric discretizations of elliptic partial differential equations, arXiv preprint arXiv:1406.6808 (2014).
  * [27] J. Henning, D. Palitta, V. Simoncini, K. Urban, An ultraweak space-time variational formulation for the wave equation: Analysis and efficient numerical solution, ESAIM: Mathematical Modelling and Numerical Analysis 56 (4) (2022) 1173–1198.
  * [28] G. Sangalli, M. Tani, Isogeometric preconditioners based on fast solvers for the Sylvester equation, SIAM Journal on Scientific Computing 38 (6) (2016) A3644–A3671.
  * [29] M. Montardini, G. Sangalli, M. Tani, Robust isogeometric preconditioners for the Stokes system based on the Fast Diagonalization method, Computer Methods in Applied Mechanics and Engineering 338 (2018) 162 – 185.
  * [30] R. Vázquez, A new design for the implementation of isogeometric analysis in Octave and Matlab: GeoPDEs 3.0, Computers & Mathematics with Applications 72 (3) (2016) 523–554.
  * [31] L. Sorber, M. Van Barel, L. De Lathauwer, Tensorlab v2. 0, Available online, URL: www.tensorlab.net (2014).
  * [32] H. Brezis, Functional analysis, Sobolev spaces and partial differential equations, Vol. 2, Springer, 2011.