# Efficient infusion of self-supervised representations in Automatic Speech
Recognition

Darshan Prabhu
Sony Research India
<EMAIL_ADDRESS>& Sai Ganesh Mirishkar
Sony Research India
<EMAIL_ADDRESS>& Pankaj Wasnik
Sony Research India
<EMAIL_ADDRESS>Work done during an internship at Sony Research India.
Corresponding to: Sai Ganesh Mirishkar

###### Abstract

Self-supervised learned (SSL) models such as Wav2vec and HuBERT yield state-
of-the-art results on speech-related tasks. Given the effectiveness of such
models, it is advantageous to use them in conventional ASR systems. While some
approaches suggest incorporating these models as a trainable encoder or a
learnable frontend, training such systems is extremely slow and requires a lot
of computation cycles. In this work, we propose two simple approaches that use
(1) framewise addition and  (2) cross-attention  mechanisms to efficiently
incorporate the representations from the SSL model(s) into the ASR
architecture, resulting in models that are comparable in size with standard
encoder-decoder conformer systems while also avoiding the usage of SSL models
during training. Our approach results in faster training and yields
significant performance gains on the Librispeech and Tedlium datasets compared
to baselines. We further provide detailed analysis and ablation studies that
demonstrate the effectiveness of our approach.

##  1 Introduction

Since speech is a complex signal, audio-related tasks like Automatic Speech
Recognition(ASR) rely heavily on having robust speech representations as
input. While standard ASR systems use conventional signal processing
techniques [13] for generating these representations, another set of models,
called representation models, are tasked with learning to generate such
representations using large unlabeled data. These models are trained using the
masked language modeling [3] objective, where the model attempts to
reconstruct parts of the text/audio that are masked. As a result, they can
understand raw text/speech well and generate high-quality text/speech
representations. BERT [3], Wav2Vec [1], and HuBERT [6] are some popular models
that fall under this category. They are also called Self-Supervised
Learned(SSL) models, as they do not need labeled data during training. These
models are popular for downstream tasks as they perform exceptionally well
even with small amounts of data. However, since these models are huge,
training an ASR system retrofitted with such models as a learnable frontend
[16] or an Encoder with limited computational resources becomes challenging,
if not impossible. To address this issue, several works propose freezing the
representation model and using its output as auxiliary information to a custom
encoder. This enables us to entirely forgo the usage of SSL models during
training, as the extraction of these representations can be offloaded to
preprocessing 111This is a one-off step performed once at the beginning of
training.. Early work in this regard is found in Natural Language
Processing(NLP), where BERT representations are integrated into the Neural
Machine Translation system [18]. Recent works perform this integration for
multimodal inputs as well [17, 7]. In speech, fine-tuned Wav2vec2 embeddings
have been used as auxiliary data alongside the domain adversarial setup [9,
10] for the accented ASR setting.

In this work, we propose two approaches that incorporate representations from
pre-trained SSL models into an end-to-end ASR architecture. We also explore
the possibility of introducing multiple such representations into the
architecture and perform a comprehensive analysis of our approach in terms of
convergence, efficiency, and understanding.

##  2 Proposed Methodology

(a) Overall Architecture

(b) Subsampled Framewise Addition (SFA)

(c) Cross Attention (CA)

Figure 1: Overview of our proposed Architecture that integrates the
representations from the self-supervised model into the ASR encoder using two
approaches: Subsampled Framewise Addition(SFA) and Cross Attention(CA).
$\hat{v}$ and $\hat{u}$ are the self-supervised and fbank representations
respectively.

Figure 1 shows our proposed modifications to the existing joint CTC-Attention
[8] framework that allows auxiliary audio representations(extracted from pre-
trained SSL models) to be easily introduced within the end-to-end ASR
architecture without substantially increasing the model size 222As the usage
of SSL model during inference is unavoidable, no computational benefit would
be observed using our approach for inference. . Figure 1(a) shows the overall
architecture that comprises three main modules: Encoder(Enc), Decoder(Dec-Att)
and Connectionist Temporal Classification(Dec-Ctc) [4]. Although all our
changes are made to the encoder, we first briefly explain the overall
architecture, followed by our proposed modifications to the encoder.

Let $\textbf{x}=\\{x_{1},x_{2},\ldots,x_{N}\\}$ be the raw representation of
the input audio. This representation $x$ is passed through  (1) Fbank
generator that generates $T$ length fbank representation
$\textbf{u}=\\{u_{1},u_{2},\ldots,u_{T}\\}$ where $u_{i}\in\mathbb{R}^{d}$ and
(2) Self supervised pre-trained model that generates $T^{\prime}$ length SSL
representation $\textbf{v}=\\{v_{1},v_{2},\ldots,v_{T^{\prime}}\\}$ where
$v_{i}\in\mathbb{R}^{d^{\prime}}$ . The encoder(Enc) jointly reasons on both u
and v to generate contextualized audio representation
$\textbf{h}=\textsc{Enc}(\textbf{u},\textbf{v})=\\{h_{1},h_{2},\ldots
h_{T}\\}$. This contextualized representation h is then consumed by Dec-Att
and Dec-Ctc module which aim at predicting the output token sequence
$\textbf{y}=\\{y_{1},y_{2},\ldots,y_{i},\ldots y_{M}\\}$ using CTC and
Attention loss.

The encoder module begins with a $\mathrm{Convolution\\_Subsampling}$ block
that applies convolution and two-factor subsampling 333Subsampling is a
crucial step that reduces the computation complexity of the encoder while
performing at par with the original encoder [2]. on the fbank representation u
resulting in a $T/2$ length sequence $\hat{\textbf{u}}$. Parallely, we employ
a linear layer and layer normalization on the SSL representation v to obtain a
$T^{\prime}$ length sequence $\hat{\textbf{v}}$. While layer normalization
helps with better generalization, the linear layer reduces the dimension from
$d^{\prime}$ to $d$. These sequences are then fed to the
$\mathrm{Fusion\\_Layer}$ that acts on both these inputs, performing a
deterministic fusion of $\hat{\textbf{u}}$ and $\hat{\textbf{v}}$ to generate
$\textbf{h}^{0}$(In section 2.1, we discuss in detail the different choices
available for the fusion layer along with their merits and demerits.). The
output of the fusion layer is then passed through a stack of $L$ identical
encoder layers. Each encoder layer feeds on the output from the previous layer
$\textbf{h}^{i-1}$ and generates contextualized representation
$\textbf{h}^{i}$ by performing the standard operation of a conformer [5]
encoder layer.

###  2.1 Fusion Layer

In this section, we discuss two simple approaches to generate the
representation $\textbf{h}^{0}$ using sequences $\hat{\textbf{u}}$ and
$\hat{\textbf{v}}$. To reiterate, u is a $T$ length sequence that is passed
through two-factor subsampling to generate the representation
$\hat{\textbf{u}}$ which is a sequence of length $T/2$. $\hat{\textbf{v}}$
which is obtained by passing raw audio through the SSL model followed by layer
normalization is a $T^{\prime}$ length sequence which in the case of Wav2Vec
and HuBERT happens to be equal to $T$. Finally, the output of this fusion
layer is $\textbf{h}^{0}$ which is also a $T/2$ length sequence.

#### Subsampled Framewise Addition (SFA):

Inspired by the work of Jialu et al. [9], we first propose a simple
parameterless approach as shown in Figure 1(b) that relies on the observation
that performing subsampling on $\hat{\textbf{v}}$ by a factor of 2
conveniently leads to both $\hat{\textbf{u}}$ and $\hat{\textbf{v}}$ to be of
equal length. We can then perform framewise addition of both these sequences
to generate $\textbf{h}^{0}$. Mathematically, this generation can be written
as:
$\textbf{h}^{0}_{i}=\hat{\textbf{u}}_{i}\oplus\hat{\textbf{v}}_{min(T,2\times
i)}(\forall\enspace 1\leq i\leq T/2)$ where $\textbf{h}^{0}_{i}$ is the
$i^{th}$ entry in $\textbf{h}^{0}$ and $\oplus$ is the elementwise addition of
two $d$-dimensional vectors. As this is a parameterless operation, it does not
add to the model size but there are two main drawbacks with this approach.
First, we have a predetermined decision on which frame of the corresponding
sequences are to be added. It would be beneficial to let the model determine
this mapping. Second, the approach heavily relies on the lengths of the two
sequences to have some linear relation, and for arbitrary lengths, it becomes
harder to determine this mapping. It becomes exponentially harder if the
lengths of both sequences differ by a large margin.

#### Cross Attention (CA):

To address both the concerns, we introduce a cross-attention layer as shown in
Figure 1(c) that uses the concept of attention to determine how each frame of
$\hat{\textbf{u}}$ wants to attend to the frames of $\hat{\textbf{v}}$.
Mathematically, this generation of $\textbf{h}^{0}$ can be written as:

$\textbf{h}^{0}_{i}=\hat{\textbf{u}}_{i}\oplus\textbf{I}_{i}\quad\forall\quad
1\leq i\leq
T/2\enspace\textrm{where}\enspace\textbf{I}=\mathrm{MultiHeadAttention}(\hat{\textbf{u}},\hat{\textbf{v}},\hat{\textbf{v}})$

where $\mathrm{MultiHeadAttention}(Q,K,V)$ refers to the standard multi-headed
attention proposed by Vaswani et al. [14] with $Q,K$ and $V$ denoting query,
key and value respectively. This approach is free from any predetermined
mapping of the frames of the corresponding sequences, as attention helps the
model learn such mappings.

##  3 Experimental Setup

We run all our experiments on NVIDIA A100 GPUs using the ESPnet toolkit [15].
As is common practice, we add 3-way speed perturbation to both datasets before
training. Unless specified otherwise, in all of our experiments, we train a
conformer model with 12 encoder and 6 decoder layers. We use 256-dimensional
tensors and four heads for attention computation. All the models are trained
for 50 epochs with a batch size of 32, a learning rate of 1.0, and 4 gradient
accumulation steps. The representations from SSL models are dumped beforehand
to speed up the training. Throughout all of our experiments, we use three
self-supervised models: (1) Wav2vec 2.0 Base, (2) HuBERT Base, and (3) HuBERT
Large. All these models are only pre-trained with the SSL objective without
any ASR fine-tuning. We report numbers on Librispeech-100 [11] and Tedlium2
[12] datasets that are publicly available.

##  4 Results

Table 1: Comparison of Performance(WER) of our system with baselines on Librispeech-100 dataset [11]. Our experiments are labeled in the format $X$+$Y$+$Z$, indicating that the architecture $X$ was modified to include representation from the $Y$ SSL model by employing the $Z$ approach. Method | Dev Clean | Dev Other | Test Clean | Test Other
---|---|---|---|---
B1: Transformer (Trans.) | 10.1 | 25.8 | 10.4 | 26.4
B2: Conformer (Conf.) | 7.9 | 21.4 | 8.4 | 22.0
T1: Conf. + w2v-Base \+ SFA | 5.9 | 16.2 | 6.4 | 16.4
T2: Conf. + w2v-Base \+ CA | 5.9 | 16.0 | 6.3 | 16.3
T3: Conf. + HuBERT-Base \+ SFA | 5.2 | 13.5 | 5.4 | 13.5
T4: Conf. + HuBERT-Base \+ CA | 5.1 | 13.0 | 5.4 | 13.3
T5: Conf. + (w2v-Base, HuBERT-Base) + CA | 5.7 | 14.1 | 6.1 | 14.5

Table 1 compares word error rates (WER) between our proposed architectures (T1
to T5) and baselines (B1 and B2) on the Librispeech 100h [11] dataset. We
report two baselines, out of which the Conformer baseline outperforms the
Transformer baseline by a large margin. We show five variants of our
architecture, labeled T1 through T5, each of which uses the conformer
architecture but differs in either the choice of SSL model or the fusion
method. We see that among Wav2vec and HuBERT, HuBERT representations
consistently perform better. Furthermore, we find no added advantage in using
both Wav2Vec and HuBERT representations. Irrespective of the SSL model chosen,
cross-attention(CA) shows significant improvement over the sub-sampled
framewise addition(SFA) approach.

##  5 Analysis

#### Faster Model Convergence:

As shown in Figure 3, our approach outperforms the baseline best results
within only ten epochs of training and continues to improve over the next few
iterations. Moreover, the convergence is much faster in the early iterations,
resulting in a reasonably good model within the first 5-10 epochs of training.

#### Efficacy of our approach:

While directly incorporating HuBERT or Wav2vec as a frontend/Encoder is bound
to perform better, our approach results in a model that is both faster in
training and smaller in size at the cost of a slight degradation in
performance. Furthermore, performing SSL feature extraction in the
preprocessing stage adds a very minimal overhead 444This process can be easily
parallelized. In our experiments, we use sixty parallel jobs to perform
feature extraction. ( approx 4 hrs for Librispeech-100 and 10 hrs for Tedlium2
) compared to using SSL models throughout training.

#### Number of encoder layers:

Table 3 compares the conformer baseline with four variations of our best
system(Conf.+HuBERT-Base+CA) which differ only in the number of encoder layers
used. Even after reducing the number of encoder layers by 80%, we find that
the model still performs much better despite having only half as many
parameters and training time as the baseline model.

Figure 2: Comparison of epochwise model performance(WER) between baseline and
our best setting(Conf.+HuBERT-Base+CA) on the validation split of
Librispeech-100 dataset.

Figure 3: Visualization of attention scores for two samples from the
Librispeech-100 dataset, one from Dev-Other split and Test-Other split
respectively. We use the model from our best setting to obtain these scores.

Table 2: Comparison of parameter count, training time and Performance(WER) of our architecture with baseline on Librispeech 100h dataset [11]. Conf.(E=x) represents conformer architecture with x Encoder layers. Method | # of params | Train time | Dev | Test
---|---|---|---|---
Clean | Other | Clean | Other
Conf.(E=12) | 30.6M | 24h | 7.9 | 21.4 | 8.4 | 22.0
Conf.(E=12) + HuBERT-Base \+ CA | 31.0M | 26h | 5.1 | 13.0 | 5.4 | 13.3
Conf.(E=8) + HuBERT-Base \+ CA | 24.7M | 22h | 5.0 | 12.9 | 5.3 | 13.0
Conf.(E=4) + HuBERT-Base \+ CA | 18.4M | 17h | 5.4 | 12.8 | 5.7 | 12.7
Conf.(E=2) + HuBERT-Base \+ CA | 15.2M | 15h | 5.5 | 12.7 | 5.5 | 12.8

Table 3: Comparison of Performance(WER) of our approach with baseline on the Dev and Test splits of Tedlium dataset [12]. We follow the same heuristics from Table 1 to label the experiments. Method | Dev | Test
---|---|---
Conformer (Conf.) | 10.5 | 8.6
Conf. + w2v-Base \+ CA | 9.6 | 9.3
Conf. + HuBERT-Base \+ CA | 9.4 | 8.9
Conf. + HuBERT-Large \+ CA | 7.6 | 6.8

#### Visualization of Attention:

Figure 3 depicts the attention distribution generated by our Cross-
Attention(CA) fusion layer for few examples. As can be seen, the layer has
learned to predict the alignment between the two representations, indicating
that approaches like Convolution or windowed attention that focuses on local
attention can also be employed to achieve reasonable improvement.

#### Performance on unseen dataset:

Table 3 compares the performance between baselines and different variations of
our proposed architecture on the Tedlium dataset. We deliberately chose the
SSL model that is not exposed to the Tedlium dataset to test the effectiveness
of the generated representations on the unseen datasets. Because the Large
model is exposed to a much more diverse dataset, we find it to be far more
effective than the Base model. Furthermore, we find normalization to be a
critical step for out-of-domain datasets.

##  6 Conclusion

In this work, we propose an end-to-end ASR architecture that integrates self-
supervised model representations into the speech encoder. We accomplish this
with a fusion layer, which can be as simple as framewise addition to a more
complex cross-attention mechanism. To evaluate the efficacy of our approach,
we conduct experiments on the Librispeech and Tedlium datasets, demonstrating
significant reductions in word error rate (WER) compared to standard baseline
models. Furthermore, we perform a thorough analysis of our approach, with a
particular focus on the convergence rate, the impact of the number of
parameters, and the information captured by the attention mechanism. Through
these detailed analyses, we highlight the speed, efficiency, and scalability
that our approach achieves in comparison to other baseline methods.

## References

  * Baevski et al. [2020] Alexei Baevski, Henry Zhou, Abdelrahman Mohamed, and Michael Auli.  wav2vec 2.0: A framework for self-supervised learning of speech representations, 2020.
  * Burchi and Vielzeuf [2021] Maxime Burchi and Valentin Vielzeuf.  Efficient conformer: Progressive downsampling and grouped attention for automatic speech recognition, 2021.
  * Devlin et al. [2019] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding, 2019.
  * Graves et al. [2006] Alex Graves, Santiago Fernández, Faustino Gomez, and Jürgen Schmidhuber.  Connectionist temporal classification: Labelling unsegmented sequence data with recurrent neural networks.  In _Proceedings of the 23rd International Conference on Machine Learning_ , ICML ’06, page 369–376, New York, NY, USA, 2006. Association for Computing Machinery.  ISBN 1595933832.  doi: 10.1145/1143844.1143891.  URL https://doi.org/10.1145/1143844.1143891.
  * Gulati et al. [2020] Anmol Gulati, James Qin, Chung-Cheng Chiu, Niki Parmar, Yu Zhang, Jiahui Yu, Wei Han, Shibo Wang, Zhengdong Zhang, Yonghui Wu, and Ruoming Pang.  Conformer: Convolution-augmented transformer for speech recognition, 2020\.
  * Hsu et al. [2021] Wei-Ning Hsu, Benjamin Bolte, Yao-Hung Hubert Tsai, Kushal Lakhotia, Ruslan Salakhutdinov, and Abdelrahman Mohamed.  Hubert: Self-supervised speech representation learning by masked prediction of hidden units, 2021.
  * Kim and Kang [2022] Donghwa Kim and Pilsung Kang.  Cross-modal distillation with audio–text fusion for fine-grained emotion classification using bert and wav2vec 2.0.  _Neurocomputing_ , 506:168–183, 2022.  ISSN 0925-2312.  doi: https://doi.org/10.1016/j.neucom.2022.07.035.  URL https://www.sciencedirect.com/science/article/pii/S0925231222008931.
  * Kim et al. [2017] Suyoun Kim, Takaaki Hori, and Shinji Watanabe.  Joint ctc-attention based end-to-end speech recognition using multi-task learning, 2017.
  * Li et al. [2021] Jialu Li, Vimal Manohar, Pooja Chitkara, Andros Tjandra, Michael Picheny, Frank Zhang, Xiaohui Zhang, and Yatharth Saraf.  Accent-robust automatic speech recognition using supervised and unsupervised wav2vec embeddings, 2021.
  * Mirishkar et al. [2023] Ganesh S Mirishkar, Vishnu Vidyadhara Raju V, Meher Dinesh Naroju, Sudhamay Maity, Prakash Yalla, and Anil Kumar Vuppala.  Iiith-cstd corpus: Crowdsourced strategies for the collection of a large-scale telugu speech corpus.  _ACM Transactions on Asian and Low-Resource Language Information Processing_ , 22(7):1–26, 2023.
  * Panayotov et al. [2015] Vassil Panayotov, Guoguo Chen, Daniel Povey, and Sanjeev Khudanpur.  Librispeech: an asr corpus based on public domain audio books.  In _Acoustics, Speech and Signal Processing (ICASSP), 2015 IEEE International Conference on_ , pages 5206–5210. IEEE, 2015.
  * Rousseau et al. [2012] Anthony Rousseau, Paul Deléglise, and Yannick Estève.  Ted-lium: an automatic speech recognition dedicated corpus.  In _Conference on Language Resources and Evaluation (LREC)_ , pages 125–129, 2012.
  * Stafford et al. [1998] Kathleen M. Stafford, Christopher G. Fox, and David S. Clark.  Long-range acoustic detection and localization of blue whale calls in the northeast Pacific Ocean.  _The Journal of the Acoustical Society of America_ , 104(6):3616–3625, 12 1998.  ISSN 0001-4966.  doi: 10.1121/1.423944.  URL https://doi.org/10.1121/1.423944.
  * Vaswani et al. [2023] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need, 2023.
  * Watanabe et al. [2018] Shinji Watanabe, Takaaki Hori, Shigeki Karita, Tomoki Hayashi, Jiro Nishitoba, Yuya Unno, Nelson Enrique Yalta Soplin, Jahn Heymann, Matthew Wiesner, Nanxin Chen, Adithya Renduchintala, and Tsubasa Ochiai.  ESPnet: End-to-end speech processing toolkit.  In _Proceedings of Interspeech_ , pages 2207–2211, 2018.  doi: 10.21437/Interspeech.2018-1456.  URL http://dx.doi.org/10.21437/Interspeech.2018-1456.
  * Zeghidour et al. [2021] Neil Zeghidour, Olivier Teboul, Félix de Chaumont Quitry, and Marco Tagliasacchi.  Leaf: A learnable frontend for audio classification, 2021.
  * Zhao et al. [2022] Zihan Zhao, Yanfeng Wang, and Yu Wang.  Multi-level fusion of wav2vec 2.0 and bert for multimodal emotion recognition, 2022.
  * Zhu et al. [2020] Jinhua Zhu, Yingce Xia, Lijun Wu, Di He, Tao Qin, Wengang Zhou, Houqiang Li, and Tie-Yan Liu.  Incorporating bert into neural machine translation, 2020.


# On accelerated Universe expansion

L. V. Verozub<EMAIL_ADDRESS>Kharkov National University

###### Abstract

It is shown that observed peculiarities of the Universe expansion are an
inevitable consequence of the gravitational force properties following from
gauge-invariant gravitation equations considered in detail in an author’s
paper in Annalen der Physik, v.17, 28 (2008).

###### pacs:

04.50+h; 98.80.-k

††preprint: APS/123-QED

Numerous data testifying that the most distant galaxies move away from us with
acceleration where obtained for the last 8 years riess . This fact poses
serious problems both for fundamental physics and astrophysics weinberg . In
the present paper it is shown that the available observational data are an
inevitable consequence of properties of the gravitational force implying from
gauge-invariant gravitation equations Verozub08a . These equations were tested
successfully by binary pulsar $PSR~{}1913+16$ verkoch00 .

In Minkowski’s space-time the radial component of the gravitational force of a
point mass $M$ affecting the free- falling particle of mass $m$ is $m\ddot{r}$
where the acceleration $\ddot{r}=d^{2}r/dt^{2}$ must be found from the
gravitation equations in use. According to Verozub08a  the force is given by

$F=-m\left[c^{2}C^{\prime}/2A+(A^{\prime}/2A-2C^{\prime}/2C)\dot{r}^{2}\right],$
(1)

where

$A=f^{\prime 2}/C,\ C=1-r_{\mathrm{g}}/f,\
f=(r_{\mathrm{g}}^{3}+r^{3})^{1/3},\,f^{\prime}=df/dt.$ (2)

In this equation the dot denotes the derivative with respect to $t$,
$r_{\mathrm{g}}=2GM/c^{2}$, $G$ is the gravitational constant, $c$ is the
speed of light at infinity, the prime denotes the derivative with respect to
$r$.

For particles at rest ($\dot{r}=0$ )

$F=-\frac{GmM}{r^{2}}\left[1-\frac{r_{\mathrm{g}}}{(r^{3}+r_{\mathrm{g}}^{3})^{1/3}}\right]$
(3)

Fig. 1 shows the force $F$ affecting a particle at rest and a particle free
falling from infinity as the function of the distance
$\overline{r}=r/r_{\mathrm{g}}$ from the centre.

Figure 1:  The gravitational force (arbitrary units) affecting a free-falling
particle (the curve 1) and a particle at rest (the curve 2) near an attractive
point mass $M$.

It follows from Fig. 1 that the gravitational force affecting free-falling
particles changes its sign at $r\approx 2r_{\mathrm{g}}$ . Although we have
never yet observed particles motion at distances of the order of $r_{g}$ we
can verify this result for very distant objects in the Universe, at large
cosmological redshifts, because it is well-known that the radius of the
observed region of the Universe is of the order of its Schwarzschild radius.

A magnitude which is related with observations in the expanding Universe is
the relative velocity of distant star objects with respect to an observer. The
radial velocity $v=\dot{R}=dR/dt$ of particles on the surface of a
selfgravitating expanding homogeneous sphere of a radius $R$ can be obtained
from equations of the motion of a test particle Verozub08a :

$v=c\frac{Cf^{2}}{R^{2}}\sqrt{1-\frac{C}{\overline{E}^{2}}},$ (4)

where $C$ are the functions of the distance $R$, $r_{g}=(8/3)\pi c^{-2}G\rho
R^{3}$ is Schwarzshild’s radius of homogeneous matter inside of the sphere and
$\rho$ is the matter density. The parameter $\overline{E}$ is the constant
total energy of a particle divided by $mc^{2}$.

Fig. 2 shows the radial acceleration $\ddot{R}=v^{\prime}\dot{v}$ of a
particle on the surface of the sphere of the radius $R$ in flat space-time

Figure 2:  The acceleration of particles on the surface of an sphere of the
radius $R$ for three value of the parameter $\overline{E}$. The matter density
is equal to $10^{-29}{g}\,{cm}^{-3}$.

Two conclusions can be made from this figure.
1\. At some distance from the observer the relative acceleration changes its
sign. If the $R<2\cdot 10^{27}{cm}$, the radial acceleration of particles is
negative. If $R>2\cdot 10^{27}{cm}$, the acceleration is positive. Hence, for
sufficiently large distances the gravitational force affecting particles is
repulsive and gives rise to a relative radial acceleration of particles with
respect to any observer.

2\. The gravitational force, affecting the particles, tends to zero when $R$
tends to infinity. (The same fact takes place as regards the force acting on
particles in the case of static matter). The reason of the fact is that at a
sufficiently large distance $R$ from the observer the Schwarzschild radius of
the matter inside of the sphere of the radius $R$ become of the same order as
its radius. Approximately at $R\sim 2\,r_{g}$ the gravitational force begin to
decrease. The ratio $R/r_{g}\rightarrow\infty$ tends to zero when $R$ tends to
infinity, and under the circumstances the gravitational force in the theory
under consideration tends to zero.

The above ball can be consedered as a part of the flat accelerating Universe
because like the case of Einstein’s equations a spherically-symmetric matter
layer does not create gravitational field inside itself Verozub08a .
Furthermore, to calculate the velocity of a particle at the distance $R$ from
the observer there is no necessity to demand a global spherical symmetry of
matter outside of the sphere due to the above second properties of the
gravitational force, because the gravitational influence of very distant
matter is neglected small. Therefore a relative velocity of particles at the
distance $R$ from the observer is determined by gravitational field of the
matter inside of the sphere of the radius $R$.

Proceeding from this equation we will find Hubble’s diagram, following mainly
the method being used in zeldovich . Let $\nu_{0}$. be a local frequency in
the proper reference frame of a moving light source at the distance $R$ from
an observer, $\nu_{l}$ be this frequency in a local inertial frame, and $\nu$
be the frequency as measured by the observer in the sphere centre. The
redshift $z=(\nu_{0}-\nu)/\nu$ is caused by both Doppler effect and
gravitational field. The Doppler effect is a consequence of a difference
between the local frequency of the source in inertial and comoving reference
frame, and it is given by landau

$\nu_{l}=\nu_{0}\,[1-\sqrt{(1-v/c)(1+v/c)}]$ (5)

The gravitational redshift is caused by the matter inside of the sphere of the
radius $R$. It is a consequence of the energy conservation for photon.
According to equations of the motion of a test particle Verozub08a  the rest
energy of a particle in gravitation field is given by

$E=mc^{2}\sqrt{C}.$ (6)

Therefore, the difference in two local level $E_{1}$ and $E_{2}$ of an atom
energy in the field is $\Delta E=(E_{2}-E_{1})\sqrt{C}$, so that the local
frequency $\nu_{0}$ at the distance $R$ from an observer are related with the
observed frequency $\nu$ by equality

$\nu=\nu_{l}\sqrt{C},$ (7)

where we take into account that for the observer location $\sqrt{C}=1$. It
follows from (7) and (5) that the relationship between frequency $\nu$ as
measured by the observer and the proper frequency $\nu_{0}$ of the moving
source in the gravitational field takes the form

$\frac{\nu}{\nu_{0}}=\sqrt{C\,\frac{1-v/c}{1+v/c}}$ (8)

Equation (8) yields the quantity $z$ as a function of $R$.

By solving this equation numerically we obtain the dependence $R=R(z)$ of the
measured distance $R$ as a function of the redshift. Therefore the distance
modulus weinberg2  to a star object is given by

$\mu=5\,log_{10}[R(z)\,(z+1)]-5$ (9)

where $R(1+z)$ is a bolometric distance (in $pc$) to the object.

If (4) is a correct equation for the radial relative velocity of distant star
objects in the expansive Universe, it must to lead to the Hubble law at small
distances R. Under this condition the Schwarzschild radius $r_{g}=(8/3)\pi
G\rho R^{3}$ of the matter inside of the sphere is very small compared with
$R$. For this reason $f\approx r$, and $C=1-r_{g}/r$. Therefore, at
$\overline{E}=1$, we obtain from (4) that

$v=HR,$ (10)

where

$H=\sqrt{(8/3)\pi G\rho}.$ (11)

If $\overline{E}\neq 1$ equation (4) does not lead to the Hubble law, since
$v$ does not tend to zero when $R\rightarrow 0$. For this reason we set
$\overline{E}=1$ and look for the value of the density at which a good
accordance with observation data can be obtained.

The fig. (3) show the Hubble diagram based on eq. (9) compared with
observations data riess .

Figure 3:  The distance modulus $\mu$ vs. the redshift $z$ for the density
$\rho=4.5\cdot 10^{-30}g\,cm^{-3}$. Small squares denote the observation data
according to Riess et al.

It follows from this figure that the model under consideration are in a good
accordance with observation data.

For the value of the density $\rho=4.5\cdot 10^{-30}g\,cm^{-3}$ we obtain from
(11) that

$H=1.59\cdot 10^{-18}c^{-1}=49\,km\,c^{-1}\,Mpc.$ (12)

Fig. 4 shows the dependence of the radial velocity $v$ on the redshift. It
follows from this figure that at $z>1$ the Universe expands with an
acceleration. At $R\rightarrow\infty$ the velocity and acceleration tend to
zero.

Figure 4:  The radial velocity vs. redshift $z$ for the density $\rho=4.5\cdot
10^{-30}g\,cm^{-3}$

## References

  * (1) A. Riess et al., ApJ 607, 665 (2004)
  * (2) S. Weinberg, E-print astro-ph/0005265 (2000)
  * (3) L.V̇erozub, Ann. Phys. (Berlin), 17, 28 (2008)
  * (4) L. Verozub & A. Kochetov, Grav. and Cosm., 6, 246 (2000)
  * (5) Ia. Zel’dovich, & I. Novikov,  Relativistic Astrophysics, v. 2: The Structure and Evolution of the Universe (University of Chicago Press, 1971)
  * (6) L. Landau & E. Lifshitz, The Classical Theory of Field (Addison - Wesley, Massachusetts, 1971)
  * (7) S. Weinberg, S. 1972, Gravitation and Cosmology, (J. Wiley & Son Inc., 1972))


# An efficient Chorin-Temam projection proper orthogonal decomposition based
reduced-order model for nonstationary Stokes equations

Xi Li111School of Mathematics, Sichuan University, Chengdu, Sichuan 610064,
China (li_xi@stu.scu.edu.cn). The work of this author was supported by the
National Natural Science Foundation of China(Grant No. 11971337)., Yan
Luo222School of Mathematical Sciences, University of Electronic Science and
Technology of China, Chengdu, Sichuan 611731, China (luoyan_16@126.com). The
work of this author was supported by the Young Scientists Fund of the National
Natural Science Foundation of China(Grant No. 11901078) and by the Fundamental
Research Funds for the Central Universities(Grant No. ZYGX2020J021). and Minfu
Feng333Corresponding author. School of Mathematics, Sichuan University,
Chengdu, Sichuan 610064, China (fmf@scu.edu.cn). The work of this author was
supported by the National Natural Science Foundation of China(Grant No.
11971337).

###### Abstract

In this paper, we propose an efficient proper orthogonal decomposition based
reduced-order model(POD-ROM) for nonstationary Stokes equations, which
combines the classical projection method with POD technique. This new scheme
mainly owns two advantages: the first one is low computational costs since the
classical projection method decouples the reduced-order velocity variable and
reduced-order pressure variable, and POD technique further improves the
computational efficiency; the second advantage consists of circumventing the
verification of classical LBB/inf-sup condition for mixed POD spaces with the
help of pressure stabilized Petrov-Galerkin(PSPG)-type projection method,
where the pressure stabilization term is inherent which allows the use of non
inf-sup stable elements without adding extra stabilization terms. We first
obtain the convergence of PSPG-type finite element projection scheme, and then
analyze the proposed projection POD-ROM’s stability and convergence. Numerical
experiments validate out theoretical results.
Keywords: Projection method; PSPG stabilization; Proper orthogonal
decomposition.

##  1 Introduction

Reduced Order Models(ROMs) have become widespread in the scientific community
and owns constant attention in academic community because of its strength and
efficiency in alleviating the huge computational costs needed in many modern
engineering applications and academic purposes. As one of the most popular ROM
approaches, the proper orthogonal decomposition(POD) strategy aims to use
singular value decomposition(SVD) to extract the dominant modes(in the sense
of energy) from a high-resolution numerical scheme(commonly referred to full
order models(FOMs)) to form POD bases/modes. Onto these reduced-order bases, a
Galerkin projection can be used to project into the space spanned by these
bases/modes to obtain the reduced-order numerical solution, this pioneering
work is done by Kunisch and Volkwein in [19, 20] which has laid the foundation
of numerical analysis for Galerkin-POD. Other references about Galerkin-POD
can be found in [30, 22, 21, 32] and so on. Applications of Galerkin-POD on
incompressible flows, especially with the spatial variables discretized by
finite element method(FE-POD), has a long history, from Luo et al.[22]
obtaining FE-POD-ROM for nonstationary Navier-Stokes equations based on finite
difference method and mixed finite element method to discretize temporal and
spatial variables, Ravindran[26] analyzing FE-POD-ROM for nonstationary
Boussinesq equations, to other finite element methods, like discontinuous
Galerkin and hybrid discontinuous Galerkin, for spatial discrete are applied
to get FE-POD-ROM[12].
Traditional mixed finite element POD-ROM had encountered a problem, that is,
when considering the classical mixed finite element scheme as the snapshot
source of the POD, the pressure variable is eliminated because of the
classical LBB/inf-sup condition, so only the velocity finite element solution
is used as the snapshot of the POD technique to construct the POD bases/modes,
which makes the final obtained FE-POD-ROM does not contain the pressure
variable. This makes it improper to use POD technique directly to reduce the
model order in some engineering and academic models that need the pressure
variables. Among that problem, the pressure field must be reconstructed a
posterior using some pressure recovery techniques, like pressure Poisson
equation[23], Supremizer stabilized technique[2] and so on. However, just as
stated in [11] and [5], there are some deficiencies in these a posterior
pressure recovery techniques, such as the usage of those techniques still need
the fulfillment of LBB/inf-sup condition which restricts some flexible pairs
of mixed finite element spaces, and the un-physical and unclear boundary
condition of pressure Poisson equation poses a challenge for the utilization
of this technique. With those in mind, recently, a new research direction
related to FE-POD has aroused some researchers’ interest, we refer as
stabilized FE-POD, which draws inspiration from the successful applications of
stabilized techniques in finite element method. Caiazzo et al.[5] firstly
proposed to combine stabilization technique, like residual-based streamline-
upwind Petrov-Galerkin(SUPG) stabilization, with FE-POD to preserve the finite
element pressure solution, so that the POD reduced-order pressure can be
obtained from the snapshot of finite element pressure solution. From then on,
other stabilized techniques were introduced into FE-POD, like local projection
stabilization(LPS) FE-POD[27, 24], artificial compressible(AC) FE-POD[11],
streamline diffusion(SD) FE-POD[1] and so on.
In this paper, our aim is to combine the classical projection method and POD
technique to propose an efficient projection POD-ROM. As one of the fast
decoupled numerical algorithms for solving incompressible fluid, projection
method has become a classical one after it was proposed by Chorin and Temam
[8, 31] and developed by Shen and Guermond [29, 15, 28, 14] etc. Apart from
that, as mentioned by Rannacher in [25], another remarkable feature of this
method is that the original Chorin-Temam projection method owns an inherent
mechanism of pressure stability. This mechanism was proposed by Rannacher [25]
and was realized by Frutos et al. [10]. These two advantages of the classical
projection method are merged into the POD technique to get our projection POD-
ROM. The main difference in our approach with respect to those stabilized FE-
POD-ROM is that the pressure stabilized Petrov-Galerkin(PSPG)-like FE-POD-ROM
has the inherent pressure stabilization, which means no additional pressure
stabilization term is needed to circumvent the LBB/inf-sup condition of mixed
POD spaces which is still an open problem to verify.
The main contribution of this paper is that a more efficient ROM numerical
scheme for solving incompressible fluid is obtained by combining the
advantageous computational efficiency of both projection method and POD
technique. This POD-ROM scheme stems from replacing the finite element
functions in the fully discrete scheme of the projection method with the POD
reduced-order functions. To this end, we analyze the fully discrete scheme of
classical projection method in the form of with inherent pressure stabilized
term, and numerically verify some theoretical results concerning with this
finite element scheme. Furthermore, we take the finite element solution as
snapshots in order to use method of snapshots to construct POD modes and POD
spaces, and obtain our projection POD scheme. The stability and convergence of
this scheme are analyzed, and which are also validated by some numerical
experiments.
The outline of this paper is as follows: In section 2, we introduce some
essential notations. In section 3, the classical projection method will be
presented, followed by its stability and convergence analysis. In section 4,
we propose the projection POD numerical scheme and prove its stability and
convergence. In section 5, some numerical experiments will be investigated to
confirm the theoritical results analyzed before. Finally, in section 6, we
draw conclusions to complete this paper.

##  2 Preliminaries and notations

we consider the nonstationary Stokes equations

$\left\\{\begin{aligned} \partial_{t}\bm{u}-\nu\Delta\bm{u}+\nabla
p&=\bm{f},\quad\text{in}\;\Omega,\\\
\nabla\cdot\bm{u}&=0,\quad\text{in}\;\Omega,\\\
\bm{u}&=\bm{0},\quad\text{on}\;\partial\Omega,\\\
\bm{u}(0,x)&=\bm{u}_{0}(x),\quad\text{in}\;\partial\Omega,\end{aligned}\right.$
(2.1)

where $\Omega$ is an open bounded domain in $\mathbb{R}^{2}$ with a
sufficiently smooth boundary $\partial\Omega$. The unknowns are the vector
function $\bm{u}$(velocity) and the scalar function $p$(pressure).
In order to state the above problem properly, we need some settings. Let
$L^{p}(\Omega),\;1\leq p\leq\infty$ denote the space of standard p-th power
absolutely integrable functions with respect to the Lebesgue measure. In
particular, $L^{2}(\Omega)$ is a Hilbert space endowed with the scalar product
$(\cdot,\cdot)$ and its induced norm $\|\cdot\|_{0}$. Moreover, we use
$\|\cdot\|^{2}_{\ell^{2}(L^{2})}:=\Delta t\sum^{N}_{n=1}\|\cdot\|^{2}_{0}$ to
denote the discrete integral in time with respect to $\|\cdot\|^{2}_{0}$. The
Sobolev space $W^{m,p}(\Omega)$ and its norm $\|\cdot\|_{m,p}$ are also
standard in the sense of [9]. Furthermore, we use the abbreviation
$H^{m}(\Omega):=W^{m,2}(\Omega)$ and $\|\cdot\|_{m}:=\|\cdot\|_{m,2}$, it’s a
Hilbert space with a scalar product. Moreover, we denote by
$H^{1}_{0}(\Omega)$ the space of functions of $H^{1}(\Omega)$ with vanishing
trace on $\partial\Omega$ and by $H^{-1}(\Omega)$ its dual space. The symbol
$\langle\cdot,\cdot\rangle$ in general denotes the duality pairing between the
space and it’s dual one. Vector analogues of the Sobolev spaces along with
vector valued functions are denoted by boldface letters, for instance
$\bm{H}^{m}(\Omega):=\left(H^{m}(\Omega)\right)^{2}$. We all know that by
Poincaré inequality[3, p.135], the seminorm $|u|_{1}=\|\nabla u\|_{0}$ in
$H^{1}(\Omega)$ is a norm in $H^{1}_{0}(\Omega)$.
Throughout this paper, we use $C$ to denote a positive constant independent of
$\Delta t$, $h$, not necessarily the same at each occurrence.
Let us introduce some convenient spaces. The first one is

$\bm{H}({\rm
div};\Omega):=\\{\bm{u}\in\bm{L}^{2}(\Omega)|\;\nabla\cdot\bm{u}\in\bm{L}^{2}({\Omega})\\},$

which is a Hilbert space with norm $\|\bm{u}\|_{{\rm
div}}:=\|\bm{u}\|_{0}+\|\nabla\cdot\bm{u}\|_{0}$, and

$\bm{H}_{0}({\rm div};\Omega):=\\{\bm{u}\in\bm{H}({\rm
div};\Omega)|\;\bm{u}\cdot\bm{n}|_{\partial\Omega}=0\\}.$

For future use, we also need the following spaces, which play a key role in
Helmholtz decomposition,

$\bm{J}_{0}:=\left\\{\bm{v}\in\bm{H}^{1}_{0}(\Omega):\;\nabla\cdot\bm{v}=0\right\\},$

and

$\bm{J}_{1}:=\left\\{\bm{v}\in\bm{L}^{2}(\Omega):\;\nabla\cdot\bm{v}=0,\;\text{and}\;\bm{v}\cdot\bm{n}|_{\partial\Omega}=0\right\\}.$

By means of spaces above, we give the classical Helmholtz decomposition[13,
p.29] as

$\bm{L}^{2}(\Omega)=\bm{J}_{1}(\Omega)\oplus[\bm{J}_{1}(\Omega)]^{\perp}=\bm{J}_{1}(\Omega)\oplus\\{\nabla
q:q\in H^{1}(\Omega)\\}.$ (2.2)

In order to give a variational formulation of problem (2.1), we consider the
velocity space $\bm{V}:=\bm{H}^{1}_{0}$ and the pressure space
$Q:=L^{2}_{0}=\\{q\in L^{2}(\Omega)|\int_{\Omega}q\,{\rm d}x=0\\}$. For a
given time constant $T$, a weak formulation of the nonstationary Stokes
equations are expressed as follows: for almost $\forall t\in(0,T]$, find
$(\bm{u},p):(0,T]\rightarrow\bm{V}\times Q$, such that

$\left\\{\begin{aligned}
(\partial_{t}\bm{u},\bm{v})+\nu(\nabla\bm{u},\nabla\bm{v})+(\nabla
p,\bm{v})&=\langle\bm{f},\bm{v}\rangle,\quad\forall\bm{v}\in\bm{V},\\\
(\nabla\cdot\bm{u},q)&=0,\quad\forall q\in Q.\\\
\bm{u}(x,0)&=\bm{u}^{0}(x).\end{aligned}\right.$ (2.3)

Let $\\{\mathcal{T}_{h}\\}$ be a uniformly regular family of triangulation of
$\overline{\Omega}$(see [9, p.111]) and
$h:=\max_{K\in\mathcal{T}_{h}}\\{h_{K}|\;h_{K}\\\ :={\rm diam}(K)\\}$. Since
we want to construct the POD space where snapshots are deriving from both the
finite element approximated velocity and also pressure, contrary to the
general choice which is to choose the inf-sup stable mixed finite element
spaces, we choose equal-order mixed finite element spaces and then it’s
obviously not stable in the sense of the classical inf-sup condition, as

$Y^{l}_{h}:=\\{v_{h}\in
C^{0}(\overline{\Omega})|\;v_{h}|_{K}\in\mathbb{P}_{l}(K),\;\forall
K\in\mathcal{T}_{h}\\},\;l\geq 1,$

where $\mathbb{P}_{l}(K)$ is the space of polynomials up to degree $l$ on $K$.
Since there are two types of velocities in the temporal semi-discrete
projection scheme, i.e., intermediate velocity $\widetilde{\bm{u}}^{n+1}$ and
end-of-step velocity $\bm{u}^{n+1}$, we need the three finite element spaces
to locate three full discrete variables of three unknowns. To this end, we
define the equal-order finite element spaces for end-of-step velocity,
intermediate velocity and pressure

$\bm{V}_{h}:=\bm{Y}^{l}_{h}\cap\bm{V},\quad\bm{W}_{h}:=\bm{Y}^{l}_{h}\cap\bm{J}_{1},\quad
Q_{h}:=Y^{l}_{h}\cap Q.$

Since the triangulations are assumed to be shape regular, the following
inverse inequality holds for each $\bm{v}_{h}\in\bm{V}_{h}$ on each mesh cell
$K\in\\{\mathcal{T}_{h}\\}$(see [3, Theorem 4.5.11]),

$\|\bm{v}_{h}\|_{m,p,K}\leq C_{{\rm
inv}}h^{l-m+2(\frac{1}{p}-\frac{1}{q})}_{K}\|\bm{v}_{h}\|_{l,q,K},$ (2.4)

where $0\leq l\leq m\leq 1$, $0\leq q\leq p\leq\infty$ and $h_{k}$ is the size
of the mesh cell $K\in\\{\mathcal{T}_{h}\\}$. Then, $I_{h}$(resp. $J_{h}$)
denotes a bounded linear interpolation operator
$I_{h}:\;\bm{W}^{t,q}(K)\rightarrow\bm{V}_{h}$(resp.
$J_{h}:\;W^{t,q}(K)\rightarrow Q_{h}$), where $1\leq q\leq\infty$, $0\leq
s\leq t\leq l+1$, that satisfies the following optimal bounds(see [3, Theorem
4.4.4])

$\displaystyle\|\bm{v}-I_{h}\bm{v}\|_{s,q,K}$ $\displaystyle\leq
Ch_{K}^{t-s}\|\bm{v}\|_{t,q,K},$ (2.5) $\displaystyle\|q-J_{h}q\|_{s,q,K}$
$\displaystyle\leq Ch_{K}^{t-s}\|q\|_{t,q,K},$

and also interpolation operator $J_{h}$ has following stability

$\|\nabla J_{h}q\|_{0}\leq C\|\nabla q\|_{0}.$ (2.6)

The Sobolev inclusion $L^{2}(\Omega)\hookrightarrow H^{-1}(\Omega)$ relation
implies the following inequality holds:

$\|v\|_{-1}\leq C\|v\|_{0},\quad\forall v\in L^{2}(\Omega).$ (2.7)

##  3 Classical projection method

###  3.1 classical projection scheme

The original projection method consists of firstly finding an intermediate
velocity $\widetilde{\bm{u}}^{n+1}$ from the momentum equation without the
pressure term, and then utilizing the classical Helmholtz decomposition (2.2)
to get an end-of-step velocity $\bm{u}^{n+1}$ which is solenoidal and its
normal component $\bm{u}^{n+1}\cdot\bm{n}$ vanishes. Meanwhile, as a “by-
product” of the above decomposition, a discrete pressure $p^{n+1}$ can be also
get for which some investigations have been made to confirm that this “by-
product” pressure is indeed the proper approximation of the exact pressure
$p(t_{n+1})$, see [25].
For a given timestep $\Delta t$, we consider a uniform discretization of the
time interval $[0,T]$ as $0=t_{0}<t_{1}<\cdots<t_{N}$ where $N=\lfloor
T/\Delta t\rfloor$ is the rounding $T/\Delta t$ down. In mathematics, the
above procedure can be interpreted as the following classical Chorin-Teman
projection method[7, 31]:
(Scheme 1 - Semi-discrete of classical projection method) For
$n=0,1,\ldots,N-1$,
Step 1\- Given $\bm{u}^{n}\in\bm{J}_{0}$, find
$\widetilde{\bm{u}}^{n+1}\in\bm{V}$ that satisfies

$\frac{\widetilde{\bm{u}}^{n+1}-\bm{u}^{n}}{\Delta
t}-\nu\Delta\widetilde{\bm{u}}^{n+1}=\bm{f}^{n+1},\quad\text{in}\;\Omega.$

Step 2\- Given $\widetilde{\bm{u}}^{n+1}\in\bm{V}$, compute
$(\bm{u}^{n+1},p^{n+1})\in\bm{H}_{0}({\rm div},\Omega)\times\left[Q\cap
H^{1}(\Omega)\right]$ from

$\left\\{\begin{aligned} \bm{u}^{n+1}+\Delta t\nabla
p^{n+1}&=\widetilde{\bm{u}}^{n+1},\quad\text{in}\;\Omega,\\\
\nabla\cdot\bm{u}^{n+1}&=0,\quad\text{in}\;\Omega,\\\
\bm{u}^{n+1}\cdot\bm{n}&=0,\quad\text{on}\;\partial\Omega.\end{aligned}\right.$

The usage of classical Helmholtz decomposition is reflected in decomposing the
intermediate velocity $\widetilde{\bm{u}}^{n+1}$ from step 1 to get a
solenoidal velocity $\bm{u}^{n+1}$ and a gradient of a $H^{1}$-function, which
can be stated as:
$\mathbf{Step\ 2^{\prime}}$\- Projecting $\widetilde{\bm{u}}^{n+1}$ onto the
space $\bm{J}_{0}$

$\bm{u}^{n+1}=P_{\bm{J}_{0}}(\widetilde{\bm{u}}^{n+1}).$

######  Remark 3.1.

Let us explain the relation between $\bm{Step2}$ and $\bm{Step2^{\prime}}$.
Followed by (2.2), decomposing $\widetilde{\bm{u}}^{n+1}$ would directly gets:

$\widetilde{\bm{u}}^{n+1}=\bm{u}^{n+1}+\nabla\phi^{n+1},$ (3.8)

where $\bm{u}^{n+1}\in\bm{J}_{1}$ implies $\nabla\cdot\bm{u}^{n+1}=0$ and
$\bm{u}^{n+1}\cdot\bm{n}|_{\partial\Omega}=0$, which is exactly the second and
third equations in $\bm{Step\ 2}$. Multiplying the previous decomposition by
$1/\Delta t$ and adding it to the first equation in $\bm{Step\ 1}$, we get
$(\bm{u}^{n+1}-\bm{u}^{n})/\Delta
t-\nu\Delta\widetilde{\bm{u}}^{n+1}+\nabla(\phi^{n+1}/\Delta t)=\bm{f}^{n+1}$,
which means $\phi^{n+1}/\Delta t$ can be regarded as a proper approximation of
$p(t_{n+1})$, i.e., we can define $p^{n+1}:=\phi^{n+1}/\Delta t$. Thus, on the
one hand, we can compute $\phi^{n+1}$ in the former decomposition (3.8) and
using $p^{n+1}:=\phi^{n+1}/\Delta t$ to get $p^{n+1}$; on the other hand, we
can also substitute $\phi^{n+1}$ with $\Delta tp^{n+1}$ in (3.8) to directly
compute $p^{n+1}$, which is exactly what the first equation in step 2 does.

Or, Step 2 is equivalent to the following pressure Poisson equation using
$\nabla\cdot\bm{u}^{n+1}=0$:
$\mathbf{Step\ 2^{\prime\prime}}$\- Given $\widetilde{\bm{u}}^{n+1}\in\bm{V}$,
compute $p^{n+1}\in Q\cap H^{2}(\Omega)$ from

$\left\\{\begin{aligned} \Delta p^{n+1}&=\frac{1}{\Delta
t}\nabla\cdot\widetilde{\bm{u}}^{n+1},\quad\text{in}\;\Omega,\\\
\frac{\partial
p^{n+1}}{\partial\bm{n}}&=0,\quad\text{on}\;\partial\Omega.\end{aligned}\right.$

Based on the previous finite element spaces, we derive the fully discrete
scheme of projection method.
(Scheme 2 - Fully-discrete of classical projection scheme) Let
$\bm{u}^{0}_{h}$ be the Lagrangian interpolation or Ritz projection onto
$\bm{V}_{h}$ of $\bm{u}_{0}$, then for $n=0,1,\ldots,N-1$, we can compute
$(\widetilde{\bm{u}}^{n+1}_{h},\bm{u}^{n+1}_{h},p^{n+1}_{h})\in$
$\bm{V}_{h}\times\bm{W}_{h}\times Q_{h}$ iteratively by

$\left\\{\begin{aligned}
\left(\frac{\widetilde{\bm{u}}^{n+1}_{h}-\bm{u}^{n}_{h}}{\Delta
t},\bm{v}_{h}\right)+\nu(\nabla\widetilde{\bm{u}}^{n+1}_{h},\nabla\bm{v}_{h})&=\langle\bm{f}^{n+1},\bm{v}_{h}\rangle,\quad\forall\bm{v}_{h}\in\bm{V}_{h},\\\
(\nabla\cdot\widetilde{\bm{u}}^{n+1}_{h},q_{h})+\Delta t(\nabla
p^{n+1}_{h},\nabla q_{h})&=0,\quad\forall q_{h}\in Q_{h},\\\
\bm{u}^{n+1}_{h}&=\widetilde{\bm{u}}^{n+1}_{h}-\Delta t\cdot\nabla
p^{n+1}_{h}.\end{aligned}\right.$ (3.9)

As pointed in [25], it is not necessary to compute the projected velocities
$\\{\bm{u}^{n}_{h}\\}^{N}_{n=0}$ since these quantities can be eliminated. To
see this, replacing $\bm{u}^{n}_{h}$ in the first equation with the third
equation at $t_{n}$ in Scheme 3, we reach the following Proj-PSPG-FE scheme.
(Scheme 3 - Proj-PSPG-FE) Let $\bm{u}^{0}_{h}$ be the Lagrangian interpolation
or Ritz projection onto $\bm{V}_{h}$ of $\bm{u}_{0}$ and $p^{0}_{h}$ denote
some proper approximation of initial pressure $p(t_{0})$, we set
$\widetilde{\bm{u}}^{0}_{h}=\bm{u}^{0}_{h}+\Delta t\nabla p^{0}_{h}$ and
replace $\bm{u}^{n}_{h}$ with the third equality in Scheme 2, we obtain the
following equivalent scheme: for $n=0,1,\ldots,N-1$, find
$(\widetilde{\bm{u}}^{n+1}_{h},p^{n+1}_{h})\in\bm{V}_{h}\times Q_{h}$

$\left\\{\begin{aligned}
\left(\frac{\widetilde{\bm{u}}^{n+1}_{h}-\widetilde{\bm{u}}^{n}_{h}}{\Delta
t},\bm{v}_{h}\right)+\nu(\nabla\widetilde{\bm{u}}^{n+1}_{h},\nabla\bm{v}_{h})+(\nabla
p^{n}_{h},\bm{v}_{h})&=\langle\bm{f}^{n+1},\bm{v}_{h}\rangle,\quad\forall\bm{v}_{h}\in\bm{V}_{h},\\\
(\nabla\cdot\widetilde{\bm{u}}^{n+1}_{h},q_{h})+\Delta t(\nabla
p^{n+1}_{h},\nabla q_{h})&=0,\quad\forall q_{h}\in Q_{h}.\end{aligned}\right.$
(3.10)

######  Remark 3.2.

For the initial discrete pressure $p^{0}_{h}$, we simply take $p^{0}_{h}=0$ in
later numerical experiments, the reasons behind this are twofold: firstly, it
can simplify the solving of initial discrete pressure $p^{0}_{h}=0$, compared
to other techniques, like pressure Poisson equation used in [16, 17];
secondly, it can also directly get the initial intermediate discrete velocity
$\widetilde{\bm{u}}^{0}_{h}$ since at this point
$\widetilde{\bm{u}}^{0}_{h}=\bm{u}^{0}_{h}$. We may also observe that this
rude treatment would make the first few time steps discrete pressure appear
oscillation. Indeed, numerical experiments have showed that the first few time
steps discrete pressure values have to be abandoned because of big discrete
errors; nevertheless, we will see in later numerical experiments that, only
after a few time steps(for example, from $n=6$ on), the discrete pressure
values can be adopted as the appropriate discrete approximation of the exact
pressure.

######  Remark 3.3.

We note that although the obtained finite element velocity solution in scheme
3 is $\widetilde{\bm{u}}^{n+1}_{h}$, rather than $\bm{u}^{n+1}_{h}$, we will
see in latter numerical experiment that when choosing $\bm{V}_{h}\times
Q_{h}=P^{1}$-$P^{1}$, the solved $\widetilde{\bm{u}}^{n+1}_{h}$ can be used as
a proper finite element approximation for the exact velocity $\bm{u}^{n+1}$.
If one wants to obtain $\bm{u}^{n+1}_{h}$, we can multiply by the last
equation in scheme 2 some test function
$\bm{w}_{h}\in\bm{W}_{h}=\bm{Y}_{h}^{1}\cap\bm{J}_{1}$ to get

$(\bm{u}^{n+1}_{h},\bm{w}_{h})=(\widetilde{\bm{u}}^{n+1}_{h},\bm{w}_{h})-\Delta
t(\nabla p^{n+1}_{h},\bm{w}_{h}).$ (3.11)

then, solving above equation can finally obtain
$\bm{u}^{n+1}_{h}\in\bm{W}_{h}=\bm{Y}_{h}^{1}\cap\bm{J}_{1}$.

We observe that, in the classical projection scheme (3.9) or (3.10), if we set
$\Delta t=\mathcal{O}(h^{2})$, the existence of PSPG-stabilized term $\Delta
t(\nabla p^{n+1}_{h},\nabla q_{h})$ makes the spatial convergence order of the
error can not achieve more than second order in $L^{2}$ norm of the velocity
and first order of the error in the $L^{2}$ norm of the pressure, so
$P^{1}$-$P^{1}$ for finite element spaces and $(\bm{u},p)\in\bm{H}^{2}\times
H^{1}$ seem to be the best choice concerning about the efficiency of
computation and the regularity for the exact solutions. Thus, We set $l=1$ in
the definition of $Y^{l}_{h}$, thus the finite element spaces for velocity and
pressure would be $\bm{V}_{h}=Y^{1}_{h}\cap\bm{V}$ and $Q_{h}=Y^{1}_{h}\cap
Q$, respectively.
In the sequel we assume $\Delta t=\mathcal{O}(h^{2})$. Specifically, we assume
there exist two positive constants $C_{1},C_{2}$, such that

$C_{1}h^{2}\leq\Delta t\leq C_{2}h^{2}.$ (3.12)

The following modified inf-sup condition relaxes the classical inf-sup
condition and makes many mixed finite element pairs “stable” in the sense of
this modified one, whose detailed proof can be found in [4, Lemma 3] or [17,
Lemma 2.1].

######  Lemma 3.1.

(modified inf-sup condition) Assuming (3.12) holds, then for
$n=0,1,\ldots,N-1$, we have the following pressure stability

$\beta\|q_{h}\|_{0}\leq Ch\|\nabla
q_{h}\|_{0}+\sup_{\bm{v}_{h}\in\bm{V}_{h}}\frac{(q_{h},\nabla\cdot\bm{v}_{h})}{\|\bm{v}_{h}\|_{1}},\quad\forall
q_{h}\in Q_{h}.$ (3.13)

###  3.2 convergence of classical projection scheme

In this subsection, we will cite some error estimation results of the finite
element solution in the classical projection scheme and omit those proofs,
since those estimates have been analyzed in [10] in detail. The necessity of
citing is based on the fact that the following POD snapshots are the finite
element solutions, which means that when we analyze the discretization error
of the POD-based reduced-order solutions, apart from the ROM truncation error,
the discretization error also involve spatial and temporal discretization
errors.

######  Lemma 3.2.

(error estimates for Proj-PSPG-FE) For $n=1,\cdots,N$, Let
$(\bm{u}^{n},p^{n})\in\bm{X}\times Q$ be the solution of continuous
variational form (2.3) at discrete time $t=t_{n}$,
$(\widetilde{\bm{u}}^{n}_{h},p^{n}_{h})\in\bm{X}_{h}\times Q_{h}$ is the
solution obtained with the PSPG-like projection scheme (3.10), and assuming
(3.12) holds, i.e. $\Delta t=\mathcal{O}(h^{2})$, then we have

$\displaystyle\|\bm{u}^{n}-\widetilde{\bm{u}}^{n}_{h}\|_{0}+h\|p^{n}-p^{n}_{h}\|_{0}+h\sqrt{\Delta
t}\|\nabla(p^{n}-p^{n}_{h})\|_{0}$ $\displaystyle\leq C(h^{2}+\Delta t).$
(3.14)
$\displaystyle\|\nabla(\bm{u}-\bm{\widetilde{u}}_{h})\|_{l^{2}(L^{2})}+\|p-p_{h}\|_{l^{2}(L^{2})}$
$\displaystyle\leq C(h+\sqrt{\Delta t}).$

In addition to the error estimates about the intermediate velocity
$\widetilde{\bm{u}}^{n}_{h}$ analyzed in [10], we can also prove the $L^{2}$
error estimates about the end-of-step velocity $\bm{u}^{n}_{h}$ which is not
made in [10]. In other words, we can obtain

######  Lemma 3.3.

For $n=1,\cdots,N$, $\bm{u}^{n}_{h}$ is the velocity solution obtained in
Scheme 2, and $\bm{u}^{n}$ is the velocity solution of continuous variational
form (2.3) at discrete time $t=t_{n}$, then

$\|\bm{u}^{n}-\bm{u}^{n}_{h}\|_{0}\leq C(h^{2}+\Delta t).$ (3.15)

###### Proof.

According to the convergent results given by Lemma 3.2, we can prove the
stability of $\|\nabla p^{n}_{h}\|_{0}$. That is, the first inequality in
(3.14) implies

$\|\nabla(p^{n}-p^{n}_{h})\|_{0}\leq C(h/\sqrt{\Delta t}+\sqrt{\Delta
t}/h)\leq C,$

where we have used (3.12). Then

$\|\nabla p^{n}_{h}\|_{0}\leq\|\nabla(p^{n}-p^{n}_{h})\|_{0}+\|\nabla
p^{n}\|_{0}\leq C,$

where the constant $C$ can be chose to stay nearly fixed with the decrease of
$h$ and $\Delta t$. So, based on the result, testing
$\bm{w}_{h}=\widetilde{\bm{u}}^{n+1}_{h}-\bm{u}^{n+1}_{h}$ in (3.11) and
rearranging to get

$\|\widetilde{\bm{u}}^{n+1}_{h}-\bm{u}^{n+1}_{h}\|_{0}\leq\Delta t\|\nabla
p^{n+1}_{h}\|_{0}\leq C\Delta t.$

Finally, utilizing triangle inequality and Lemma 3.2,

$\|\bm{u}^{n+1}-\bm{u}^{n+1}_{h}\|_{0}\leq\|\bm{u}^{n+1}-\widetilde{\bm{u}}^{n+1}_{h}\|_{0}+\|\widetilde{\bm{u}}^{n+1}_{h}-\bm{u}^{n+1}_{h}\|_{0}\leq
C(h^{2}+\Delta t).$

∎

##  4 Projection POD

###  4.1 POD

In this subsection, we will briefly give some essential ingredients for
constructing POD method, more details for constructing POD space and POD-based
reduced-order models can be found in [21, 32]. As we remarked in Remark 3.2,
the numerical pressure oscillation in the first few time steps caused by
setting $p^{0}_{h}=0$ makes us have to discard the finite element solutions in
the previous time steps, so when using finite element solution obtained by
solving Scheme 4 to construct POD spaces, we also need to consider this
situations. In view of the fact that the different choice of initial pressure
values will affect the discrete error of finite element pressure in the
previous time steps, we uniformly denote $n=n_{0},n_{0}\geq 1$ as the starting
time steps of taking snapshots in Scheme 4; therefore, for a given positive
integer $M$, we choose the finite element solution
$(\widetilde{\bm{u}}^{i}_{h},p^{i}_{h}),\;n_{0}\leq i\leq n_{0}+M-1$ in Scheme
4, and its difference quotients(DQs)
$(\partial\widetilde{\bm{u}}^{i}_{h},\partial p^{i}_{h}),\;n_{0}+1\leq i\leq
n_{0}+M-1$, where the DQs $\partial f^{i}$ are defined by $\partial
f^{i}:=(f^{i}-f^{i-1})/\Delta t$ for some function $f$ at discrete time
$t=t_{i}$, to formulate the snapshot spaces

$\displaystyle\widetilde{\mathcal{U}}$
$\displaystyle:=\langle\widetilde{\bm{u}}^{n_{0}}_{h},\widetilde{\bm{u}}^{n_{0}+1}_{h},\cdots,\widetilde{\bm{u}}^{M+n_{0}-1}_{h},\partial\widetilde{\bm{u}}^{n_{0}+1}_{h},\partial\widetilde{\bm{u}}^{n_{0}+2}_{h},\cdots,\partial\widetilde{\bm{u}}^{M+n_{0}-1}_{h}\rangle,$
$\displaystyle\mathcal{P}$ $\displaystyle:=\langle
p^{n_{0}}_{h},p^{n_{0}+1}_{h},\cdots,p^{M+n_{0}-1}_{h},\partial
p^{n_{0}+1}_{h},\partial p^{n_{0}+2}_{h},\cdots,\partial
p^{M+n_{0}+1}_{h}\rangle.$

where $\langle S\rangle$ denotes the space spanned by the set $S$ and we
denote by $N_{s}=2M-1$ the number of snapshots. Let $d_{\widetilde{u}},d_{p}$
be the dimensions of the spaces $\widetilde{\mathcal{U}}$ and $\mathcal{P}$,
respectively. The symbols $K_{\widetilde{u}},K_{p}$ are correlation matrices
corresponding to the snapshots
$K_{\widetilde{u}}=\left((K^{\widetilde{u}}_{i,j})\right)\in\mathbb{R}^{N_{s}\times
N_{s}}$ and $K_{p}=\left((K^{p}_{i,j})\right)\in\mathbb{R}^{N_{s}\times
N_{s}}$ where

$K^{\widetilde{u}}_{i,j}:=\left(\widetilde{\bm{u}}^{i}_{h},\widetilde{\bm{u}}^{j}_{h}\right),\quad
K^{p}_{i,j}:=\left(p^{i}_{h},p^{j}_{h}\right).$

and $(\cdot,\cdot)$ is the $L^{2}$ inner product. Just as in [19], a singular
value decomposition(SVD) is carried out and the leading generalized
eigenfunctions are chosen as bases, referred to as the POD bases. We denote by
$\lambda_{1}\geq\lambda_{2}\geq\cdots\lambda_{d_{\widetilde{u}}}>0$ the
positive eigenvalues of $K_{\widetilde{u}}$ and by
$\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{d_{\widetilde{u}}}\in\mathbb{R}^{N_{s}}$
the associated eigenvectors. Analogously,
$\\{\gamma_{i},\bm{y}_{i}\\}^{d_{p}}_{i=1}$ are the eigen-pairs of $K_{p}$.
Then, the two POD bases can be written explicitly as

$\displaystyle\bm{\varphi}_{i}$
$\displaystyle=\frac{1}{\sqrt{\lambda_{i}}}\left(\sum^{n_{0}+M-1}_{j=n_{0}}x^{j}_{i}\widetilde{\bm{u}}^{j}_{h}+\sum^{n_{0}+M-1}_{j=n_{0}+1}x^{j+M-1}_{i}\partial\widetilde{\bm{u}}^{j}_{h}\right),\;$
$\displaystyle\text{for}\;i=1,2,\cdots,d_{\widetilde{u}},$ (4.16)
$\displaystyle\psi_{i}$
$\displaystyle=\frac{1}{\sqrt{\gamma_{i}}}\left(\sum^{n_{0}+M-1}_{j=n_{0}}y^{j}_{i}p^{j}_{h}+\sum^{n_{0}+M-1}_{j=n_{0}}y^{j+M-1}_{i}\partial
p^{j}_{h}\right),\;$ $\displaystyle\text{for}\;i=1,2,\cdots,d_{p}.$

In what follows we will denote by

$\displaystyle\bm{V}_{r}=\langle\bm{\varphi}_{1},\bm{\varphi}_{2},\cdots,\bm{\varphi}_{r}\rangle,\quad
Q_{r}=\langle\psi_{1},\psi_{2},\cdots,\psi_{r}\rangle.$

i.e., we choose the first $r$ POD bases to span the POD spaces, and define the
traditional $L^{2}$ projection into the POD velocity space $\bm{V}_{r}$ and
POD pressure space $Q_{r}$ as

######  Definition 4.1.

Let $\Pi^{v}_{r}:L^{2}(\Omega)\rightarrow\bm{V}_{r}$ and
$\Pi^{q}_{r}:L^{2}(\Omega)\rightarrow Q_{r}$ such that

$\displaystyle(\bm{u}-\Pi^{v}_{r}\bm{u},\bm{v}_{r})$
$\displaystyle=0\quad\forall\bm{v}_{r}\in\bm{V}_{r}\quad and$
$\displaystyle(p-\Pi^{q}_{r}p,q_{r})$ $\displaystyle=0\quad\forall q_{r}\in
Q_{r}.$

Then, we can get the following Proj-POD-ROM scheme:
(Scheme 4 - Proj-POD-ROM) For some $n_{0}\geq 1$, let
$(\widetilde{\bm{u}}^{n_{0}}_{r},p^{n_{0}}_{r})=(\Pi^{v}_{r}\widetilde{\bm{u}}^{n_{0}}_{h},\Pi^{q}_{r}p^{n_{0}}_{h})=(\sum^{r}_{i=1}(\widetilde{\bm{u}}^{n_{0}}_{h},\bm{\varphi}_{i})\bm{\varphi}_{i}$,
$\sum^{r}_{i=1}(p^{n_{0}}_{h},\psi_{i})\psi_{i})$ be the Galerkin projection
of $(\widetilde{\bm{u}}^{n_{0}}_{h},p^{n_{0}}_{h})$ onto POD spaces
$\bm{V}_{r}\times Q_{r}$, then for $n=n_{0},n_{0}+1,n_{0}+2,\cdots,N-1$, we
can compute $(\widetilde{\bm{u}}^{n+1}_{r},p^{n+1}_{r})\in\bm{V}_{r}\times
Q_{r}$ from

$\left\\{\begin{aligned}
\left(\frac{\widetilde{\bm{u}}^{n+1}_{r}-\widetilde{\bm{u}}^{n}_{r}}{\Delta
t},\bm{v}_{r}\right)+\nu(\nabla\widetilde{\bm{u}}^{n+1}_{r},\nabla\bm{v}_{r})+(\nabla
p^{n}_{r},\bm{v}_{r})&=\langle\bm{f}^{n+1},\bm{v}_{r}\rangle,\quad\forall\bm{v}_{r}\in\bm{V}_{r},\\\
(\nabla\cdot\widetilde{\bm{u}}^{n+1}_{r},q_{r})+\Delta t(\nabla
p^{n+1}_{r},\nabla q_{r})&=0,\quad\forall q_{r}\in Q_{r}.\end{aligned}\right.$
(4.17)

###  4.2 stability and convergence analysis of projection-POD scheme

We will carry on some numerical analysis about the newly proposed Proj-POD-
ROM. To this end, we first recall the classical conclusions concerning about
the projection error between the snapshots solution and POD-based reduced-
order solution(see [19])

######  Lemma 4.1.

(time-averaged projection error estimates) We have the following time-averaged
projection error estimates.

$\displaystyle\frac{1}{{N_{s}}}\sum^{n_{0}+M-1}_{j=n_{0}}\left\|\widetilde{\bm{u}}^{j}_{h}-\sum^{r}_{k=1}\left(\widetilde{\bm{u}}^{j}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}\right\|^{2}_{0}+\frac{1}{{N_{s}}}\sum^{n_{0}+M-1}_{j=n_{0}+1}\left\|\partial\widetilde{\bm{u}}^{j}_{h}-\sum^{r}_{k=1}\left(\partial\widetilde{\bm{u}}^{j}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}\right\|^{2}_{0}$
$\displaystyle=\sum^{d_{\widetilde{u}}}_{k=r+1}\lambda_{k},$ (4.18)
$\displaystyle\frac{1}{{N_{s}}}\sum^{n_{0}+M-1}_{j=n_{0}}\left\|p^{j}_{h}-\sum^{r}_{k=1}\left(p^{j}_{h},\psi_{k}\right)\phi_{k}\right\|^{2}_{0}+\frac{1}{{N_{s}}}\sum^{n_{0}+M-1}_{j=n_{0}+1}\left\|\partial
p^{j}_{h}-\sum^{r}_{k=1}\left(\partial
p^{j}_{h},\psi_{k}\right)\phi_{k}\right\|^{2}_{0}$
$\displaystyle=\sum^{d_{p}}_{k=r+1}\gamma_{k}.$

######  Remark 4.1.

We remark that choosing $L^{2}$ inner product to construct POD bases and by
virtue of Parseval’s identity makes the above three identities actually
represent $L^{2}$ POD projection errors.

Apart from the above time-averaged projection error estimates, we also need
the following up-to-data optimal in time error estimates(see, [18]), which is
presented in the form of hypothesis in the previous papers and is proved to be
valid in the presence of DQs in [18].

######  Lemma 4.2.

(POD Optimal pointwise in time error estimate) We have the following optimal
in time projection errors

$\displaystyle\max_{0\leq k\leq
N}\|\bm{\widetilde{u}}^{k}_{h}-\Pi^{v}_{r}\bm{\widetilde{u}}^{k}_{h}\|^{2}_{0}$
$\displaystyle\leq C\sum^{d_{\widetilde{u}}}_{i=r+1}\lambda_{i},$
$\displaystyle\max_{0\leq k\leq N}\|p^{k}_{h}-\Pi^{q}_{r}p^{k}_{h}\|^{2}_{0}$
$\displaystyle\leq C\sum^{d_{p}}_{i=r+1}\gamma_{i},$

where $C=6\max\\{1,T^{2}\\}$, and

$\displaystyle\max_{0\leq k\leq
N}\|\nabla(\bm{\widetilde{u}}^{k}_{h}-\Pi^{v}_{r}\bm{\widetilde{u}}^{k}_{h})\|^{2}_{0}$
$\displaystyle\leq
C\sum^{d_{\widetilde{u}}}_{i=r+1}\lambda_{i}\|\nabla\bm{\varphi}_{i}\|^{2}_{0},$
$\displaystyle\max_{0\leq k\leq
N}\|\nabla(p^{k}_{h}-\Pi^{q}_{r}p^{k}_{h})\|^{2}_{0}$ $\displaystyle\leq
C\sum^{d_{p}}_{i=r+1}\gamma_{i}\|\nabla\psi_{i}\|^{2}_{0},$

where $C=6\max\\{1,T^{2}\\}$.

We also need in the later analysis the time-averaged optimal projection error
estimate about DQs, which is derived in [19].

######  Lemma 4.3.

(POD Optimal time-averaged error estimate about DQs) For the difference
quotients we have the estimate:

$\frac{1}{M}\sum^{n_{0}+M-1}_{k=n_{0}}\|\partial\bm{\widetilde{u}}^{k}_{h}-\Pi^{v}_{r}\partial\bm{\widetilde{u}}^{k}_{h}\|^{2}_{0}\leq
C\sum^{d_{\widetilde{u}}}_{i=r+1}\lambda_{i}.$

The following analogues of finite element inverse inequalities will be used in
later analysis(see, [19])

######  Lemma 4.4.

(POD inverse estimates) Denote by
$S_{\widetilde{u}}=\left((s^{\widetilde{u}}_{i,j})\right)\in\mathbb{R}^{{N_{s}}\times{N_{s}}}$
and $S_{p}=\left((s^{p}_{i,j})\right)\in\mathbb{R}^{{N_{s}}\times{N_{s}}}$ are
stiffness matrices for POD spaces $\bm{V}_{r},Q_{r}$ respectively, with
entries
$s^{\widetilde{u}}_{i,j}=(\nabla\bm{\varphi}_{i},\nabla\bm{\varphi}_{j}),s^{p}_{i,j}=(\nabla\psi_{i},\nabla\psi_{j})$,
and by
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}$
the spectral norm of the matrix, then

$\displaystyle\|\nabla\widetilde{\bm{v}}\|_{0}$
$\displaystyle\leq\sqrt{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{S}_{\widetilde{u}}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}}\|\widetilde{\bm{v}}\|_{0},\
\text{for}\ \forall\widetilde{\bm{v}}\in\bm{V}_{r},$ (4.19)
$\displaystyle\|\nabla q\|_{0}$
$\displaystyle\leq\sqrt{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|S_{p}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}}\|q\|_{0},\
\text{for}\;\forall q\in Q_{r}.$

######  Remark 4.2.

Mass matrix $M$, for both $(\bm{\varphi}_{i},\bm{\varphi}_{j})$ and
$(\psi_{i},\psi_{j})$, equals identity $I$ since the standard orthogonality of
POD base in the sense of $L^{2}$ inner product.

We can prove the unconditional stability of the Proj-POD-ROM scheme.

######  Theorem 4.1.

(Unconditional stability of Proj-POD-ROM) We have the following unconditional
stability

$\displaystyle\left\|\widetilde{\bm{u}}_{r}^{N}\right\|_{0}^{2}+\sum_{n=n_{0}}^{N-1}\left\|\widetilde{\bm{u}}_{r}^{n+1}-\widetilde{\bm{u}}_{r}^{n}\right\|_{0}^{2}$
$\displaystyle+\Delta
t\sum_{n=n_{0}}^{N-1}\left(\nu\left\|\nabla\widetilde{\bm{u}}_{r}^{n+1}\right\|_{0}^{2}+\Delta
t\left\|\nabla p_{r}^{n+1}\right\|_{0}^{2},\right)$ (4.20) $\displaystyle\leq$
$\displaystyle
C\left(\left\|\widetilde{\bm{u}}_{r}^{0}\right\|_{0}^{2}+\frac{\Delta
t}{\nu}\sum_{n=n_{0}}^{N-1}\|\bm{f}^{n+1}\|_{-1}\right).$

###### Proof.

Taking $(\bm{v}_{r},q_{r})=(\widetilde{\bm{u}}^{n+1}_{r},p^{n}_{r})$ in (4.17)
and multiply by $2\Delta t$, adding both equations and integrating by parts,
adding and subtracting $2\Delta t^{2}\|\nabla p_{r}^{n+1}\|^{2}_{0}$, we get

$\displaystyle\left\|\widetilde{\bm{u}}_{r}^{n+1}\right\|_{0}^{2}-\left\|\widetilde{\bm{u}}_{r}^{n}\right\|_{0}^{2}$
$\displaystyle+\left\|\widetilde{\bm{u}}_{r}^{n+1}-\widetilde{\bm{u}}_{r}^{n}\right\|_{0}^{2}+2\nu\Delta
t\left\|\nabla\widetilde{\bm{u}}_{r}^{n+1}\right\|_{0}^{2}+2\Delta
t^{2}\|\nabla p_{r}^{n+1}\|^{2}_{0},$ (4.21) $\displaystyle\leq C\frac{\Delta
t}{\nu}\|\bm{f}^{n+1}\|_{-1}+\nu\Delta
t\left\|\nabla\widetilde{\bm{u}}_{r}^{n+1}\right\|_{0}^{2}+2\Delta
t^{2}\left(\nabla p_{r}^{n+1},\nabla(p_{r}^{n+1}-p_{r}^{n})\right).$

For the last term above, we utilize the second equation in (4.17) twice to
obtain

$\Delta t\left(\nabla
p_{r}^{n+1},\nabla(p_{r}^{n+1}-p_{r}^{n})\right)=-(\nabla\cdot(\widetilde{\bm{u}}^{n+1}_{r}-\widetilde{\bm{u}}^{n}_{r}),p_{r}^{n+1})=(\widetilde{\bm{u}}^{n+1}_{r}-\widetilde{\bm{u}}^{n}_{r},\nabla
p_{r}^{n+1}).$

This equality implies

$2\Delta t^{2}\left(\nabla
p_{r}^{n+1},\nabla(p_{r}^{n+1}-p_{r}^{n})\right)\leq\frac{2}{3}\|\widetilde{\bm{u}}^{n+1}_{r}-\widetilde{\bm{u}}^{n}_{r}\|^{2}_{0}+\frac{3}{2}\Delta
t^{2}\|\nabla p_{r}^{n+1}\|^{2}_{0}.$

Inserting the above equality into (4.21) and adding $n$ from $0$ to $N-1$ we
finally obtain (4.20). ∎

We have the following truncation error estimate between the POD reduced-order
solution and the snapshot data.

######  Lemma 4.5.

(POD truncation error) For $n=1,2,\cdots,N$,
$(\widetilde{\bm{u}}^{n}_{h},p^{n}_{h})$ is the solution of finite element
projection scheme (3.10) and $(\widetilde{\bm{u}}^{n}_{r},p^{n}_{r})$ is the
solution of POD projection scheme (4.17), then we have the following error
estimate:

$\displaystyle\|\widetilde{\bm{u}}^{N}_{h}-\widetilde{\bm{u}}^{N}_{r}\|^{2}_{0}$
$\displaystyle+\nu\Delta
t\sum^{N-1}_{n=n_{0}}\|\nabla(\widetilde{\bm{u}}^{n+1}_{h}-\widetilde{\bm{u}}^{n+1}_{r})\|^{2}_{0}+\Delta
t^{2}\sum^{N-1}_{n=n_{0}}\|\nabla(p^{n+1}_{h}-p^{n+1}_{r})\|^{2}_{0},$ (4.22)
$\displaystyle\leq
C\left(\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}\Big{(}1+(\nu+1)\|\nabla\varphi_{i}\|^{2}_{0}\Big{)}+(\Delta
t+1)\sum^{d_{p}}_{i=r+1}\epsilon_{i}\|\nabla\psi_{i}\|^{2}_{0}+\Delta
t^{2}\right).$

###### Proof.

Finite element projection scheme (3.10) subtract POD projection scheme (4.17)
to get:

$\left\\{\begin{aligned}
\left(\frac{\widetilde{\bm{u}}^{n+1}_{h}-\widetilde{\bm{u}}^{n}_{h}}{\Delta
t}-\frac{\widetilde{\bm{u}}^{n+1}_{r}-\widetilde{\bm{u}}^{n}_{r}}{\Delta
t},\bm{v}_{r}\right)+\nu(\nabla(\widetilde{\bm{u}}^{n+1}_{h}-\widetilde{\bm{u}}^{n+1}_{r}),\nabla\bm{v}_{r})+(\nabla(p^{n}_{h}-p^{n}_{r}),\bm{v}_{r})&=0,\\\
(\nabla\cdot(\widetilde{\bm{u}}^{n+1}_{h}-\widetilde{\bm{u}}^{n+1}_{r}),q_{r})+\Delta
t(\nabla(p^{n+1}_{h}-p^{n+1}_{r}),\nabla q_{r})&=0.\end{aligned}\right.$
(4.23)

Denoting

$\displaystyle\widetilde{\bm{u}}^{n+1}_{h}-\widetilde{\bm{u}}^{n+1}_{r}=\widetilde{\bm{u}}^{n+1}_{h}-\Pi^{v}_{r}\widetilde{\bm{u}}^{n+1}_{h}+\Pi^{v}_{r}\widetilde{\bm{u}}^{n+1}_{h}-\widetilde{\bm{u}}^{n+1}_{r}$
$\displaystyle:=\widetilde{\bm{\eta}}^{n+1}_{u,r}+\widetilde{\bm{w}}^{n+1}_{u,r},$
$\displaystyle
p^{n+1}_{h}-p^{n+1}_{r}=p^{n+1}_{h}-\Pi^{q}_{r}p^{n+1}_{h}+\Pi^{q}_{r}p^{n+1}_{h}-p^{n+1}_{r}$
$\displaystyle:=\eta^{n+1}_{p,r}+w^{n+1}_{p,r}.$

Testing $(\bm{v}_{r},q_{r})=(\widetilde{\bm{w}}^{n+1}_{u,r},w^{n+1}_{p,r})$,
we obtain

$\left\\{\begin{aligned}
\left(\frac{\widetilde{\bm{w}}^{n+1}_{u,r}-\widetilde{\bm{w}}^{n}_{u,r}}{\Delta
t},\widetilde{\bm{w}}^{n+1}_{u,r}\right)&+\nu\|\nabla\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0}+(\nabla
w^{n+1}_{p,r},\widetilde{\bm{w}}^{n+1}_{u,r})=\left(\frac{\widetilde{\bm{\eta}}^{n+1}_{u,r}-\widetilde{\bm{\eta}}^{n}_{u,r}}{\Delta
t},\widetilde{\bm{w}}^{n+1}_{u,r}\right)+\nu(\nabla\widetilde{\bm{\eta}}^{n+1}_{u,r},\nabla\widetilde{\bm{w}}^{n+1}_{u,r}),\\\
&\quad+(\nabla(p^{n+1}_{h}-p^{n}_{h}),\widetilde{\bm{w}}^{n+1}_{u,r})+(\nabla(p^{n}_{r}-p^{n+1}_{r}),\widetilde{\bm{w}}^{n+1}_{u,r})+(\nabla\eta^{n+1}_{p,r},\widetilde{\bm{w}}^{n+1}_{u,r}),\\\
(\nabla\cdot\widetilde{\bm{w}}^{n+1}_{u,r},w^{n+1}_{p,r})&+\Delta t\|\nabla
w^{n+1}_{p,r}\|^{2}_{0}=(\nabla\cdot\widetilde{\bm{\eta}}^{n+1}_{u,r},\widetilde{\bm{w}}^{n+1}_{u,r})+\Delta
t(\nabla\eta^{n+1}_{p,r},\nabla w^{n+1}_{p,r}).\end{aligned}\right.$

Adding both equations we have the following error equation

$\displaystyle\left(\frac{\widetilde{\bm{w}}^{n+1}_{u,r}-\widetilde{\bm{w}}^{n}_{u,r}}{\Delta
t},\widetilde{\bm{w}}^{n+1}_{u,r}\right)$
$\displaystyle+\nu\|\nabla\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0}+\Delta
t\|\nabla
w^{n+1}_{p,r}\|^{2}_{0}=\left(\frac{\widetilde{\bm{\eta}}^{n+1}_{u,r}-\widetilde{\bm{\eta}}^{n}_{u,r}}{\Delta
t},\widetilde{\bm{w}}^{n+1}_{u,r}\right),$ (4.24)
$\displaystyle+\nu(\nabla\widetilde{\bm{\eta}}^{n+1}_{u,r},\nabla\widetilde{\bm{w}}^{n+1}_{u,r})+(\nabla(p^{n+1}_{h}-p^{n}_{h}),\widetilde{\bm{w}}^{n+1}_{u,r})+(\nabla(p^{n}_{r}-p^{n+1}_{r}),\widetilde{\bm{w}}^{n+1}_{u,r}),$
$\displaystyle+(\nabla\eta^{n+1}_{p,r},\widetilde{\bm{w}}^{n+1}_{u,r})+(\nabla\cdot\widetilde{\bm{\eta}}^{n+1}_{u,r},\widetilde{\bm{w}}^{n+1}_{u,r})+\Delta
t(\nabla\eta^{n+1}_{p,r},\nabla w^{n+1}_{p,r}),$
$\displaystyle:=A_{1}+A_{2}+A_{3}+A_{4}+A_{5}+A_{6}+A_{7}.$

We will bound the seven terms above separately.

$\displaystyle|A_{1}|$
$\displaystyle=\left|\left(\frac{\widetilde{\bm{\eta}}^{n+1}_{u,r}-\widetilde{\bm{\eta}}^{n}_{u,r}}{\Delta
t},\widetilde{\bm{w}}^{n+1}_{u,r}\right)\right|\leq
C\left\|\frac{\widetilde{\bm{\eta}}^{n+1}_{u,r}-\widetilde{\bm{\eta}}^{n}_{u,r}}{\Delta
t}\right\|^{2}_{0}+\frac{1}{10}\|\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0}.$

For the first term above, we utilize the expression of
$\widetilde{\bm{u}}^{n}_{h}$ and $\Pi^{v}_{r}\widetilde{\bm{u}}^{n}_{h}$ in
the form of POD bases $\bm{\varphi}_{i}$ to get

$\displaystyle\widetilde{\bm{\eta}}^{n+1}_{u,r}=\widetilde{\bm{u}}^{n+1}_{h}-\Pi^{v}_{r}\widetilde{\bm{u}}^{n+1}_{h}=\sum^{d_{\widetilde{u}}}_{k=1}\left(\widetilde{\bm{u}}^{n+1}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}-\sum^{r}_{k=1}\left(\widetilde{\bm{u}}^{n+1}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}=\sum^{d_{\widetilde{u}}}_{k=r+1}\left(\widetilde{\bm{u}}^{n+1}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}.$

Then

$\displaystyle\left\|\frac{\widetilde{\bm{\eta}}^{n+1}_{u,r}-\widetilde{\bm{\eta}}^{n}_{u,r}}{\Delta
t}\right\|^{2}_{0}$ $\displaystyle=\left\|\frac{1}{\Delta
t}\left(\sum^{d_{\widetilde{u}}}_{k=r+1}\left(\widetilde{\bm{u}}^{n+1}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}-\sum^{d_{\widetilde{u}}}_{k=r+1}\left(\widetilde{\bm{u}}^{n}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}\right)\right\|^{2}_{0},$
$\displaystyle=\left\|\sum^{d_{\widetilde{u}}}_{k=r+1}\left(\frac{\widetilde{\bm{u}}^{n}_{h}-\widetilde{\bm{u}}^{n}_{h}}{\Delta
t},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}\right\|^{2}_{0},$
$\displaystyle=\left\|\sum^{d_{\widetilde{u}}}_{k=r+1}\left(\partial\widetilde{\bm{u}}^{n+1}_{h},\bm{\varphi}_{k}\right)\bm{\varphi}_{k}\right\|^{2}_{0},$
$\displaystyle=\left\|\partial\widetilde{\bm{u}}^{n+1}_{h}-\Pi^{v}_{r}\partial\widetilde{\bm{u}}^{n+1}_{h}\right\|^{2}_{0}.$

For other terms on the right-side hand of (4.24),

$\displaystyle|A_{2}|$
$\displaystyle=\left|\nu(\nabla\widetilde{\bm{\eta}}^{n+1}_{u,r},\nabla\widetilde{\bm{w}}^{n+1}_{u,r})\right|\leq
C\nu\|\nabla\widetilde{\bm{\eta}}^{n+1}_{u,r}\|^{2}_{0}+\frac{1}{2}\nu\|\nabla\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0}\leq
C\nu\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}\|\nabla\varphi_{i}\|^{2}_{0}+\frac{1}{2}\nu\|\nabla\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0},$
$\displaystyle|A_{3}|$
$\displaystyle=\left|(\nabla(p^{n+1}_{h}-p^{n}_{h}),\widetilde{\bm{w}}^{n+1}_{u,r})\right|\leq
C\|\nabla(p^{n+1}_{h}-p^{n}_{h})\|^{2}_{0}+\frac{1}{10}\|\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0}\leq
C\Delta t\int^{t_{n+1}}_{t_{n}}\|\nabla(\partial_{t}p_{h})\|^{2}_{0}\,{\rm
d}t+\frac{1}{10}\|\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0},$
$\displaystyle|A_{4}|$
$\displaystyle=\left|(\nabla(p^{n}_{r}-p^{n+1}_{r}),\widetilde{\bm{w}}^{n+1}_{u,r})\right|\leq
C\|\nabla(p^{n+1}_{r}-p^{n}_{r})\|^{2}_{0}+\frac{1}{10}\|\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0}\leq
C\Delta t\int^{t_{n+1}}_{t_{n}}\|\nabla(\partial_{t}p_{r})\|^{2}_{0}\,{\rm
d}t+\frac{1}{10}\|\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0},$
$\displaystyle|A_{5}|$
$\displaystyle=\left|(\nabla\eta^{n+1}_{p,r},\widetilde{\bm{w}}^{n+1}_{u,r})\right|\leq
C\|\nabla\eta^{n+1}_{p,r}\|^{2}_{0}+\frac{1}{10}\|\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0}\leq
C\sum^{d_{p}}_{i=r+1}\epsilon_{i}\|\nabla\psi_{i}\|^{2}_{0}+\frac{1}{10}\|\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0},$
$\displaystyle|A_{6}|$
$\displaystyle=\left|(\nabla\cdot\widetilde{\bm{\eta}}^{n+1}_{u,r},\widetilde{\bm{w}}^{n+1}_{u,r})\right|\leq
C|(\nabla\widetilde{\bm{\eta}}^{n+1}_{u,r},\widetilde{\bm{w}}^{n+1}_{u,r})|\leq
C\|\nabla\widetilde{\bm{\eta}}^{n+1}_{u,r}\|^{2}_{0}+\frac{1}{10}\|\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0}\leq
C\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}\|\nabla\varphi\|^{2}_{0}+\frac{1}{10}\|\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0},$
$\displaystyle|A_{7}|$ $\displaystyle=\left|\Delta
t(\nabla\eta^{n+1}_{p,r},\nabla w^{n+1}_{p,r})\right|\leq\frac{1}{2}\Delta
t\|\nabla\eta^{n+1}_{p,r}\|^{2}_{0}+\frac{1}{2}\Delta t\|\nabla
w^{n+1}_{p,r}\|^{2}_{0}\leq C\Delta
t\sum^{d_{p}}_{i=r+1}\epsilon_{i}\|\nabla\psi_{i}\|^{2}_{0}+\frac{1}{2}\Delta
t\|\nabla w^{n+1}_{p,r}\|^{2}_{0}.$

Combining with all the seven terms’ results, we have

$\displaystyle\left(\frac{\widetilde{\bm{w}}^{n+1}_{u,r}-\widetilde{\bm{w}}^{n}_{u,r}}{\Delta
t},\widetilde{\bm{w}}^{n+1}_{u,r}\right)$
$\displaystyle+\frac{1}{2}\nu\|\nabla\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0}+\frac{1}{2}\Delta
t\|\nabla
w^{n+1}_{p,r}\|^{2}_{0}\leq\left\|\partial\widetilde{\bm{u}}^{n+1}_{h}-\Pi^{v}_{r}\partial\widetilde{\bm{u}}^{n+1}_{h}\right\|^{2}_{0},$
$\displaystyle\quad+C\Delta
t\int^{t_{n+1}}_{t_{n}}\Big{(}\|\partial_{t}p_{h}\|^{2}_{0}+\|\partial_{t}p_{r}\|^{2}_{0}\Big{)}\,{\rm
d}t+C(\nu+1)\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}\|\nabla\varphi_{i}\|^{2}_{0},$
$\displaystyle\quad+C(\Delta
t+1)\sum^{d_{p}}_{i=r+1}\epsilon_{i}\|\nabla\psi_{i}\|^{2}_{0}+\frac{1}{2}\|\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0}.$

Multiplying by $\Delta t$, adding with $n$ from $n_{0}$ to $N-1$ and
rearranging, we get

$\displaystyle\|\widetilde{\bm{w}}^{N}_{u,r}\|^{2}_{0}$ $\displaystyle+\Delta
t\sum^{N-1}_{n=n_{0}}\|\widetilde{\bm{w}}^{n+1}_{u,r}-\widetilde{\bm{w}}^{n}_{u,r}\|^{2}_{0}+\nu\Delta
t\sum^{N-1}_{n=n_{0}}\|\nabla\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0}+\Delta
t^{2}\sum^{N-1}_{n=n_{0}}\|\nabla w^{n+1}_{p,r}\|^{2}_{0},$
$\displaystyle\quad\leq\Delta
t\sum^{N-1}_{n=n_{0}}\left\|\partial\widetilde{\bm{u}}^{n+1}_{h}-\Pi^{v}_{r}\partial\widetilde{\bm{u}}^{n+1}_{h}\right\|^{2}_{0}+C(\nu+1)\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}\|\nabla\varphi_{i}\|^{2}_{0}+\frac{1}{2}\Delta
t\sum^{N-1}_{n=n_{0}}\|\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0},$
$\displaystyle\quad+C(\Delta
t+1)\sum^{d_{p}}_{i=r+1}\epsilon_{i}\|\nabla\psi_{i}\|^{2}_{0}+C\Delta
t^{2}\int^{T}_{t_{n_{0}}}\Big{(}\|\partial_{t}p_{h}\|^{2}_{0}+\|\partial_{t}p_{r}\|^{2}_{0}\Big{)}\,{\rm
d}t.$

For the first term on the right-hand side above, we use Lemma 4.3 to have

$\displaystyle\Delta
t\sum^{N-1}_{n=n_{0}}\left\|\partial\widetilde{\bm{u}}^{n+1}_{h}-\Pi^{v}_{r}\partial\widetilde{\bm{u}}^{n+1}_{h}\right\|^{2}_{0}\leq
C\frac{1}{M}\sum^{M}_{k=1}\|\partial\widetilde{\bm{u}}^{k}_{h}-\Pi^{v}_{r}\partial\widetilde{\bm{u}}^{k}_{h}\|^{2}_{0}\leq
C\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}.$

By discrete Gronwall inequality and using stability result to bound the last
term above, we get

$\displaystyle\|\widetilde{\bm{w}}^{N}_{u,r}\|^{2}_{0}$ $\displaystyle+\Delta
t\sum^{N-1}_{n=n_{0}}\|\widetilde{\bm{w}}^{n+1}_{u,r}-\widetilde{\bm{w}}^{n}_{u,r}\|^{2}_{0}+\nu\Delta
t\sum^{N-1}_{n=n_{0}}\|\nabla\widetilde{\bm{w}}^{n+1}_{u,r}\|^{2}_{0}+\Delta
t^{2}\sum^{N-1}_{n=n_{0}}\|\nabla w^{n+1}_{p,r}\|^{2}_{0},$ (4.25)
$\displaystyle\quad\leq
C\left(\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}\Big{(}1+(\nu+1)\|\nabla\varphi_{i}\|^{2}_{0}\Big{)}+(\Delta
t+1)\sum^{d_{p}}_{i=r+1}\epsilon_{i}\|\nabla\psi_{i}\|^{2}_{0}+\Delta
t^{2}\right).$

∎

Finally, by triangle inequality, Lemma 3.2 and Lemma 4.5, we obtain the
following theorem which states the convergence between continuous variational
form and POD projection scheme.

######  Theorem 4.2.

(error estimate for Proj-POD-ROM) For $n=n_{0},n_{0}+1,\cdots,N$, let
$(\bm{u}^{n},p^{n})$ is the solution of (2.3) at $t=t_{n}$,
$(\widetilde{\bm{u}}^{n}_{r},p^{n}_{r})$ denotes the POD-based reduced-order
solutions obtained in (4.17), then we have the following convergence estimate:

$\displaystyle\left\|\bm{u}^{N}-\widetilde{\bm{u}}^{N}_{r}\right\|^{2}_{0}$
$\displaystyle+\nu\Delta
t\sum^{N-1}_{n=n_{0}}\left\|\nabla(\widetilde{\bm{u}}^{n}-\widetilde{\bm{u}}^{n}_{r})\right\|^{2}_{L^{2}}+\Delta
t^{2}\sum^{N-1}_{n=n_{0}}\left\|\nabla(p-p_{r})\right\|^{2}_{L^{2}},$ (4.26)
$\displaystyle\leq
C\left(\sum^{d_{\widetilde{u}}}_{i=r+1}\gamma_{i}\Big{(}1+(\nu+1)\left\|\nabla\varphi_{i}\right\|^{2}_{0}\Big{)}+(\Delta
t+1)\sum^{d_{p}}_{i=r+1}\epsilon_{i}\left\|\nabla\psi_{i}\right\|^{2}_{0}+\Delta
t^{2}+h^{2}\right).$

##  5 Numerical Tests

In this section, we present some numerical experiments to confirm the
$a\,priori$ error estimates derived in Theorem 4.2 for the POD projection
scheme 4 - Proj-POD-ROM. To this end, we first determine some necessary
parameters used before and verify the convergence in Lemma 3.2 and Lemma 3.3
of the finite element projection scheme, since with which we formed our POD
snapshot spaces, and also the convergence result in Theorem 4.2 is closely
related to the one of the finite element scheme. The open-source finite
element package iFEM [6] has been used to run the numerical experiments.

###  5.1 Problem setting

We take $\Omega=[0,1]\times[0,1]\subset\mathbb{R}^{2}$ and the time interval
$(0,1]$ with the viscosity coefficient $\nu=1$ and the prescribed solution

$\displaystyle\bm{u}$
$\displaystyle=\cos(t)\cdot\left(\begin{array}[]{c}\pi\sin(\pi x)^{2}\sin(2\pi
y)\\\ -\pi\sin(2\pi x)\sin(\pi y)^{2}\\\ \end{array}\right),$ $\displaystyle
p$ $\displaystyle=\cos(t)\cdot 10\cos(\pi x)\cos(\pi y).$

The right-hand side, the Dirichlet boundary condition and the initial
condition are chosen in accordance to the above prescribed solution.
We set $\Delta t=0.1h^{2}$ and all grids are regular $N\times N$ triangular
grids with SWNE diagonals for different $N$(i.e., diagonals coming from
connecting the southwest and the northeast vertexes on all rectangles), and we
take $N=4,8,16,32,64$ sequentially. For simplicity, we use $P^{1}$-$P^{1}$
element pair for spatial discretization, and we take the snapshots on the
finest computational mesh, i.e., h=1/64, and show this mesh in Figure 1.
When constructing the snapshots spaces, we take into account the fact that the
finite element projection scheme involves the initial pressure $p^{0}_{h}$
which is not the part of the definition of the problem, so we take
$p^{0}_{h}=0$ and it might be good to think of the numerical discrete errors
of pressure $\|p^{n}-p^{n}_{h}\|_{L^{2}}$ in the previous steps $n$ are so
large that it is not suitable to take finite element solution on those steps
into snapshot spaces, thus we choose the finite element solution
$(\bm{u}^{n}_{h},p^{n}_{h})$ and its difference quotients
$(\partial\bm{u}^{n}_{h},\partial p^{n}_{h})$ from $n=n_{0}=6$ to $n=25$ to
formulate snapshots spaces, which means we take $M=20$ and thus the number of
snapshots is $N_{s}=2M-1=39$.
The finite element solution
$(\bm{u}^{n}_{h},\widetilde{\bm{u}}^{n}_{j},p^{n}_{h})$ and exact solution
$(\bm{u}^{n},p^{n})$ at discrete termination time $t=\Delta t\cdot N$ on
finest mesh are shown in Figure 2, and following POD bases are formed from
snapshots via $L^{2}$ inner product.

Figure 1: Computational mesh with h=1/64. Figure 2: The exact solution and
finite element solution at $n=N$ on mesh of $h=1/64$. On upper row, from left
to right: ($u^{n}_{1},\widetilde{u}^{n}_{1,h},u^{n}_{1,h},p^{n}$); on lower
row, from left to right:
($u^{n}_{2},\widetilde{u}^{n}_{2,h},u^{n}_{2,h},p^{n}_{h}$).

###  5.2 Convergence of the finite element projection scheme

In this subsection, numerical tests will be used to numerically verify two
theoretical results, that is Lemma 3.2 and Lemma 3.3, with which we utilize to
get convergence analysis Theorem 4.2 of POD scheme 4. Although Lemma 3.2 has
been analyzed in [10], it lacks the corresponding numerical verifications, so
we supplement the numerical tests about Lemma 3.2 and also provide that of the
Lemma 3.3, which is newly analyzed in this paper. Since $\Delta
t=\mathcal{O}(h^{2})$, then from the Table 1 and Table 2 we can see

$\displaystyle\|\nabla(\bm{u}-\bm{\widetilde{u}}_{h})\|_{l^{2}(L^{2})}+\|p-p_{h}\|_{l^{2}(L^{2})}$
$\displaystyle=\mathcal{O}(h),$
$\displaystyle\|\bm{u}^{n}-\widetilde{\bm{u}}^{n}_{h}\|_{0}+h\|p^{n}-p^{n}_{h}\|_{0}+h\sqrt{\Delta
t}\|\nabla(p^{n}-p^{n}_{h})\|_{0}$ $\displaystyle=\mathcal{O}(h^{2}).$

which are what the theoretical results in Lemma 3.2 shows. Especially, we can
see from Table 1,
$\left\|\bm{u}^{n}-\bm{u}^{n}_{h}\right\|_{L^{2}}=\mathcal{O}(h^{2})$, which
is in accordance of Lemma 3.3.

Table 1: Errors and convergence orders of finite element solution $\widetilde{\bm{u}}^{n}_{h},\bm{u}^{n}_{h}$ with $P^{1}$-$P^{1}$ pair. 1/h | $\max\limits_{n}\left\|\bm{u}^{n}-\widetilde{\bm{u}}^{n}_{h}\right\|_{L^{2}}$ | $\max\limits_{n}\left\|\bm{u}^{n}-\bm{u}^{n}_{h}\right\|_{L^{2}}$ | $\left\|\nabla(\bm{u}-\widetilde{\bm{u}}_{h})\right\|_{\ell^{2}(L^{2})}$
---|---|---|---
error | rate | error | rate | error | rate
4 |  5.3013e$-$01 |  - |  4.6509e$-$01 |  - |  4.8187e$+$00 |  -
8 |  1.6490e$-$01 | 1.7078 | 1.4931e$-$01 | 1.6392 | 2.6626e$+$00 | 0.85582
16 | 4.3368e$-$02 | 1.9259 | 4.0108e$-$02 | 1.8963 | 1.3785e$+$00 | 0.94976
32 | 1.0969e$-$02 | 1.9828 | 1.0527e$-$02 | 1.9298 | 7.1098e$-$01 | 0.95516
64 | 2.7499e$-$03 | 1.9960 | 2.8273e$-$03 | 1.8966 | 3.7409e$-$01 | 0.92642
Table 2: Errors and convergence orders of finite element pressure $p^{n}_{h}$ with $P^{1}$-$P^{1}$ pair. 1/h | $\max\limits_{n}\|p^{n}-p^{n}_{h}\|_{L^{2}}$ | $\left\|p-p_{h}\right\|_{\ell^{2}(L^{2})}$ | $\sqrt{\Delta t}\left\|\nabla(p-p_{h})\right\|_{\ell^{2}(L^{2})}$
---|---|---|---
error | rate | error | rate | error | rate
4 |  2.7636e$+$00 |  - |  2.2987e$+$00 |  - |  1.4286e$+$00 |  -
8 | 1.1144e$+$00 | 1.3103 | 8.8892e$-$01 | 1.3707 | 4.6814e$-$01 | 1.6096
16 | 3.6664e$-$01 | 1.6038 | 2.7275e$-$01 | 1.7045 | 1.3827e$-$01 | 1.7595
32 | 1.2463e$-$01 | 1.5567 | 8.1260e$-$02 | 1.7470 | 4.5158e$-$01 | 1.6144
64 | 4.6335e$-$02 | 1.4275 | 2.5152e$-$02 | 1.6919 | 1.5553e$-$01 | 1.5378

###  5.3 Convergence of projection POD

After having the snapshots spaces, we determine
$d_{\widetilde{u}}=\text{rank}(\widetilde{\mathcal{U}})=20$(for which
$\gamma_{i}<10^{-15},i>d_{\widetilde{u}}$),
$d_{p}=\text{rank}(\mathcal{P})=20$(for which
$\epsilon_{i}<10^{-13},i>d_{p}$). Figure 3 shows the decay of POD
eigenvalues(left) of velocity $\gamma_{i},i=1,\cdots,d_{\widetilde{u}}$ and
pressure $\epsilon_{i},i=1,\cdots,d_{p}$, together with the corresponding
captured system’s energy(right) in the form of
$100\sum^{r}_{i=1}\gamma_{i}/\sum^{d_{\widetilde{u}}}_{i=1}\gamma_{i}$ for
velocity and $100\sum^{r}_{i=1}\epsilon_{i}/\sum^{d_{p}}_{i=1}\epsilon_{i}$
for pressure. We note that the first $r=4$ POD modes already capture more than
$99.99\%$ of the system’s velocity-pressure energy.
In order to confirm the computational efficiency advantage of the POD-ROM
scheme over the finite element FOM scheme, we take the number of POD modes is
$r=4$, whereas the number of degree of freedoms(DOFs) for velocity and
pressure in finite element scheme are 8450 and 4225. We show in Table 3 the
comparison of cumulative time spent by the two schemes to run to some certain
time in the process of time-stepping iterative, and the $L^{2}$ numerical
spatial discrete error of velocity and pressure in two schemes at that time.
We can see from Table 3 that, compared with the finite element FOM scheme,
POD-ROM scheme can effectively improve computational efficiency in the context
of only consuming approximately 1/6 times to get the numerical solution with
even higher accuracy.
Figure 4 plots the temporal evolution of the discrete $L^{2}$ relative
error(in semilogarithmic scale) of the reduced-order velocity and pressure
with respect to the full order ones:
$\|\widetilde{\bm{u}}^{n}_{h}-\widetilde{\bm{u}}^{n}_{r}\|_{L^{2}}/\|\widetilde{\bm{u}}^{n}_{h}\|_{L^{2}}$
and $\|p^{n}_{h}–p^{n}_{r}\|_{L^{2}}/\|p^{n}_{h}\|_{L^{2}}$ for different
number of POD modes. As expected in Theorem 4.2, the errors would decrease as
the number of POD modes $r$ increasing. For velocity, the left figure in
Figure 4 numerically demonstrate this result, and for pressure, the error
would increase slightly when $r$ from $6$ to $8$, but when we continue to
increase $r$, the error will decrease, which is consistent with the Theorem
4.2 as a whole.

Figure 3: POD velocity-pressure eigenvalues(left) and captured system’s velocity-pressure energy(right). Table 3: Comparison of errors and time consumed by time-stepping iterative between finite element scheme with $P^{1}$-$P^{1}$ pair and POD scheme with $r=4$ on the mesh of $h=1/64$ and $\Delta t=0.1h^{2}$. n | finite element scheme | POD scheme
---|---|---
$\|\bm{u}^{n}\\!-\\!\widetilde{\bm{u}}^{n}_{h}\|_{L^{2}}$ | $\|p^{n}\\!-\\!p^{n}_{h}\|_{L^{2}}$ | CPU run time(s) | $\|\bm{u}^{n}\\!-\\!\widetilde{\bm{u}}^{n}_{r}\|_{L^{2}}$ | $\|p^{n}\\!-\\!p^{n}_{r}\|_{L^{2}}$ | CPU run time(s)
2500 | 2.3789e$-$03 | 2.9458e$-$02 | 689 | 2.2860e$-$03 | 2.7823e$-$02 | 113
5000 | 2.3929e$-$03 | 2.9253e$-$02 | 1338 | 2.3007e$-$03 | 2.7642e$-$02 | 227
7500 | 2.3740e$-$03 | 2.8975e$-$02 | 2005 | 2.2826e$-$03 | 2.7379e$-$02 | 339
10000 | 2.3452e$-$03 | 2.8591e$-$02 | 2668 | 2.2549e$-$03 | 2.7015e$-$02 | 452
20000 | 2.1443e$-$03 | 2.6010e$-$02 | 5339 | 2.0618e$-$03 | 2.4573e$-$02 | 914
30000 | 1.8163e$-$03 | 2.1886e$-$02 | 8060 | 1.7464e$-$03 | 2.0673e$-$02 | 1399
40000 | 1.3805e$-$03 | 1.6464e$-$02 | 10746 | 1.3274e$-$03 | 1.5547e$-$02 | 1900
Figure 4: Temporal evolution of $L^{2}$ relative error of velocity
$\bm{u}^{n}_{r}$ with respect to $\bm{u}^{n}_{h}$(left) and velocity
$p^{n}_{r}$ with respect to $p^{n}_{h}$(right) .

##  6 Conclusions

In this paper, we proposed an efficient projection POD-ROM, which combined the
advantages of classical projection method and POD technique.
The main contribution of the present paper consisted of two aspects: the first
one was high computational efficiency. Through auxiliary intermediate velocity
variable, the classical projection method decoupled the velocity variable and
pressure variable, meanwhile decoupled the saddle-point system arose from
Stokes equations, so one strength of projection method lied in its high
computational efficiency, or low computational costs; Furthermore, POD
technique was utilized to get the ROM, which have made the newly proposed
projection POD-ROM had high computational efficiency. The second contribution
was based on the fact that, in the fully discrete scheme of the classical
projection method, the original scheme could be rewritten into a PSPG-type
pressure stabilization scheme by eliminating the end-of-step velocity, where
the pressure stabilized term $\Delta t(\nabla p^{n+1}_{h},\nabla
q^{n+1}_{h})$, $\Delta t=\mathcal{O}(h^{2})$, was inherent, so that some
flexible mixed finite element spaces pairs(for example, $P^{1}$-$P^{1}$ pair)
could be used without considering the classical LBB/inf-sup condition for
mixed POD spaces, which was different from other stabilized FE-POD-ROM to
overcome LBB/inf-sup condition by adding extra stabilization terms.
Numerical experiments have been conducted to confirm the convergence for both
projection finite element scheme and projection POD scheme. We first
numerically confirmed the PSPG-type classical projection owned the desired
convergence orders consistent with theoretical analysis, and then, finite
element solutions were taken as the snapshots to formulate POD bases/modes
which are used to get reduced-order velocity and pressure variables, so in
theory, the discrete error between exact solutions and reduced-order solutions
should converge to zero as increasing the number of POD modes, we numerically
verify this result. Apart from the projection POD-ROM convergence, we also
conduct the experiment to compare the error and computational time between
projection finite element FOM and projection POD-ROM, and the results revealed
the fact that projection POD-ROM not only had less discretization error, but
also less computational costs, compared with projection finite element FOM.
One future research direction will be the applied of projection POD on
nonlinear nonstationary Navier-Stokes equations. Another investigation is
validation of fulfillment of LBB/inf-sup condition for mixed POD spaces.

## Acknowledgment

The authors thank the anonymous referees for their constructive comments and
suggestions which improved the manuscript.

## References

  * [1] Mejdi Azaïez, Tomás Chacón Rebollo, and Samuele Rubino.  A cure for instabilities due to advection-dominance in POD solution to advection-diffusion-reaction equations.  J. Comput. Phys., 425:Paper No. 109916, 27, 2021.
  * [2] Francesco Ballarin, Andrea Manzoni, Alfio Quarteroni, and Gianluigi Rozza.  Supremizer stabilization of POD-Galerkin approximation of parametrized steady incompressible Navier-Stokes equations.  Internat. J. Numer. Methods Engrg., 102(5):1136–1161, 2015.
  * [3] Susanne C. Brenner and L. Ridgway Scott.  The mathematical theory of finite element methods, volume 15 of Texts in Applied Mathematics.  Springer, New York, third edition, 2008.
  * [4] Erik Burman and Miguel A. Fernández.  Analysis of the PSPG method for the transient Stokes’ problem.  Comput. Methods Appl. Mech. Engrg., 200(41-44):2882–2890, 2011\.
  * [5] Alfonso Caiazzo, Traian Iliescu, Volker John, and Swetlana Schyschlowa.  A numerical investigation of velocity-pressure reduced order models for incompressible flows.  J. Comput. Phys., 259:598–616, 2014.
  * [6] Long Chen.  ifem: an innovative finite element methods package in matlab.  Preprint, University of Maryland, 2008.
  * [7] Alexandre Joel Chorin.  Numerical solution of the Navier-Stokes equations.  Math. Comp., 22:745–762, 1968.
  * [8] Alexandre Joel Chorin.  On the convergence of discrete approximations to the Navier-Stokes equations.  Math. Comp., 23:341–353, 1969.
  * [9] Philippe G. Ciarlet.  The finite element method for elliptic problems, volume 40 of Classics in Applied Mathematics.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002.
  * [10] Javier de Frutos, Bosco García-Archilla, and Julia Novo.  Error analysis of projection methods for non inf-sup stable mixed finite elements. The transient Stokes problem.  Appl. Math. Comput., 322:154–173, 2018.
  * [11] Victor DeCaria, Traian Iliescu, William Layton, Michael McLaughlin, and Michael Schneier.  An artificial compression reduced order model.  SIAM J. Numer. Anal., 58(1):565–589, 2020.
  * [12] Guosheng Fu and Zhu Wang.  POD-(H)DG method for incompressible flow simulations.  J. Sci. Comput., 85(2):Paper No. 24, 20, 2020.
  * [13] Vivette Girault and Pierre-Arnaud Raviart.  Finite element methods for Navier-Stokes equations - Theory and algorithms, volume 5 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, 1986.
  * [14] J. L. Guermond, P. Minev, and Jie Shen.  An overview of projection methods for incompressible flows.  Comput. Methods Appl. Mech. Engrg., 195(44-47):6011–6045, 2006\.
  * [15] J.-L. Guermond and L. Quartapelle.  On the approximation of the unsteady Navier-Stokes equations by finite element projection methods.  Numer. Math., 80(2):207–238, 1998.
  * [16] John G. Heywood and Rolf Rannacher.  Finite element approximation of the nonstationary Navier-Stokes problem. I. Regularity of solutions and second-order error estimates for spatial discretization.  SIAM J. Numer. Anal., 19(2):275–311, 1982.
  * [17] Volker John and Julia Novo.  Analysis of the pressure stabilized Petrov-Galerkin method for the evolutionary Stokes equations avoiding time step restrictions.  SIAM J. Numer. Anal., 53(2):1005–1031, 2015.
  * [18] Birgul Koc, Samuele Rubino, Michael Schneier, John Singler, and Traian Iliescu.  On optimal pointwise in time error bounds and difference quotients for the proper orthogonal decomposition.  SIAM J. Numer. Anal., 59(4):2163–2196, 2021.
  * [19] K. Kunisch and S. Volkwein.  Galerkin proper orthogonal decomposition methods for parabolic problems.  Numer. Math., 90(1):117–148, 2001.
  * [20] K. Kunisch and S. Volkwein.  Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics.  SIAM J. Numer. Anal., 40(2):492–515, 2002.
  * [21] Zhendong Luo and Goong Chen.  Proper orthogonal decomposition methods for partial differential equations.  Mathematics in Science and Engineering. Elsevier/Academic Press, London, 2019.
  * [22] Zhendong Luo, Jing Chen, I. M. Navon, and Xiaozhong Yang.  Mixed finite element formulation and error estimates based on proper orthogonal decomposition for the nonstationary Navier-Stokes equations.  SIAM J. Numer. Anal., 47(1):1–19, 2008/09.
  * [23] Bernd R. Noack, Paul Papas, and Peter A. Monkewitz.  The need for a pressure-term representation in empirical galerkin models of incompressible shear flows.  J. Fluid Mech., 523:339–365, 2005.
  * [24] Julia Novo and Samuele Rubino.  Error analysis of proper orthogonal decomposition stabilized methods for incompressible flows.  SIAM J. Numer. Anal., 59(1):334–369, 2021.
  * [25] Rolf Rannacher.  On Chorin’s projection method for the incompressible Navier-Stokes equations.  In The Navier-Stokes equations II—theory and numerical methods (Oberwolfach, 1991), volume 1530 of Lecture Notes in Math., pages 167–183. Springer, Berlin, 1992.
  * [26] S. S. Ravindran.  Error analysis for Galerkin POD approximation of the nonstationary Boussinesq equations.  Numer. Methods Partial Differential Equations, 27(6):1639–1665, 2011.
  * [27] Samuele Rubino.  Numerical analysis of a projection-based stabilized POD-ROM for incompressible flows.  SIAM J. Numer. Anal., 58(4):2019–2058, 2020.
  * [28] Jie Shen.  On error estimates of projection methods for Navier-Stokes equations: first-order schemes.  SIAM J. Numer. Anal., 29(1):57–77, 1992.
  * [29] Jie Shen.  On error estimates of some higher order projection and penalty-projection methods for Navier-Stokes equations.  Numer. Math., 62(1):49–73, 1992.
  * [30] John R. Singler.  New POD error expressions, error bounds, and asymptotic results for reduced order models of parabolic PDEs.  SIAM J. Numer. Anal., 52(2):852–876, 2014.
  * [31] Roger Temam.  Une méthode d’approximation de la solution des équations de Navier-Stokes.  Bull. Soc. Math. France, 96:115–152, 1968.
  * [32] S. Volkwein.  Model reduction using proper orthogonal decomposition.  Lecture Notes, Faculty of Mathematics and Statistics, University of Konstanz, 2011.