A New Reduced Basis Method for Parabolic Equations Based on Single-Eigenvalue Acceleration

The Reduced Basis (RB) method is a type of model order reduction approach for numerical approximation of problems involving repeated solution of differential equations generated in engineering and applied sciences. It was first proposed in [2] for the analysis of nonlinear structures in the 1970s, and later extended to many other problems such as partial differential equations (PDEs) with multiple parameters or different initial conditions [9, 4, 41, 10, 11, 37, 12, 39], PDE-constrained parametric optimization and control problems [16, 34, 17, 18], and inverse problems [24, 31]. It should be mentioned that the work in [38, 48] has led to a decisive improvement in the computational aspects of RB methods, owing to an efficient criterion for the selection of the basis functions, a systematic splitting of the computational procedure into an offline (parameter-independent) phase and an online (parameter-dependent) phase, and the use of posteriori error estimates that guarantee certified numerical solutions for the reduced problems. The RB methods that include the offline and online phases as their essential constituents have become the most widely used ones.

The Proper Orthogonal Decomposition (POD) method, combined with the Galerkin projection method, is a typical RB method. It uses a set of orthonormal bases, which can represent the known data in the sense of least squares, to linearly approximate the target variables so as to obtain a low-dimensional approximate model. Since it is optimal in the least square sense, the POD method has the property of completely relying, without making any a priori assumption, on the data.

The POD method, whose predecessor was initially presented by K.Pearson[36] in 1901 as an eigenvector analysis method and was originally conceived in the framework of continuous second-order processes by Berkooz[5], has been widely used in many fields with different appellations. In the singular value analysis and sample identification, the method is called the Karhumen-Loeve expansion [21]. In statistics it is named as the principal component analysis (PCA) [19]. In geophysical fluid dynamics and meteorological sciences, it is called the empirical orthogonal function method (EOF)[33, 42]. We can also see the widespread applications of POD method in fluid dynamics and viscous structures[3, 8, 15, 25, 35, 40, 43, 44, 45], optimal fluid control problems [20, 30], numerical analysis of PDEs [26, 22, 23, 1, 6, 27, 29, 28] and machine learning [32, 47, 13, 7, 46].

In this paper, we develop a new RB/POD method, named as Single Eigenvalue Acceleration Method (SEAM), for a full discretization, using backward Euler scheme for temporal discretization and continuous simplicial finite elements for spatial discretization, of second order parabolic equations with homogeneous Dirichlet boundary conditions. The idea of SEAM is inspired by a POD numerical experiment for a one-dimensional heat conduction problem:

When we focus on the numerical solution matrix $U$ of (1), where each column consists of the value of the numerical solution at a node, and the number of columns is the same as the number of time steps (We divide $(0,T]$ into 1000 equal parts and $D$ into 99 equal parts to get the high-fidelity numerical solution; see Figure 1), we observe something interesting: the principal singular value is much larger than other singular values of the numerical solution matrix! Here we recall that the singular values of $U$ are the arithmetic square roots of the eigenvalues of $U^{\top}U$. Notice that when we use the POD method the choice of basis functions is often empirical. The observation tells us that we only need to choose one basis function when solving this equation with POD.

We show in Figure 1 , Figure 2 and Table 1 the SEAM (POD with one basis function) solution and computational details, respectively. We can see that when SEAM is used, the computational time can be saved greatly, and the relative error between the SEAM solution and the high-fidelity solution is very small.

Based on the above observation from the numerical experiment of 1d case, we investigate SEAM for d-dimensional (d=2,3) second order parabolic problems with homogeneous Dirichlet boundary conditions. Under the assumption that the time step size is sufficiently small and the time steps are not very large, we show that the singular value distribution of the high-fidelity solution matrix is very close to a reference matrix of rank one. Then we select the eigenfunction associated to the principal eigenvalue as the POD basis so as to obtain SEAM.

The rest of the paper is arranged as follows. Section 2 introduces the model problem and its full discretization. Section 3 gives the SEAM algorithm and its parallel version. Section 4 discusses the singular value distribution of numerical solution matrix, and Section 5 reports some numerical experiments. Finally, Section 6 gives concluding remarks.

Let $D\subset\mathcal{R}^{d}$ ($d=1,2,3$) be an open, bounded domain, and set $D_{T}:=D\times(0,T]$ for constant $T>0$. We consider the following second-order parabolic problem: find $u=u(x,t)$ such that

where $\alpha(x)=\left(\alpha_{k,j}(x)\right)_{d\times d}$ with $\alpha_{k,j}(x)=\alpha_{j,k}(x)\in W^{1,\infty}(D)$, $c(x)\in L^{\infty}(D)$, and $f=f(x,t):D_{T}\rightarrow\mathcal{R}$ and $u_{0}(x):D\rightarrow\mathcal{R}$ are given functions.

Let $\mathcal{T}_{h}=\bigcup\{T\}$ be a conforming shape-regular simplicial decomposition of the domain $D$, and let $\mathbb{V}_{h}\subset H_{0}^{1}(D)$ be a space of standard finite elements of arbitrary order with respect to $\mathcal{T}_{h}$. Denote by $\left\{x_{k}:\ k=1,2,\cdots,M\right\}$ the set of all interior nodes of $\mathcal{T}_{h}$, and by $\{\varphi_{k}(x):\ k=1,2,\cdots,M\}$ the nodal basis of $\mathbb{V}_{h}$ with

Then we have

We divide the time interval $(0,T)$ into a uniform grid with nodes $t_{n}=n\tau$ and step size $\tau=T/N$ for $n=0,1,\ldots,N$. For convenience we set $g^{n}:=g(t_{n})$ for any function $g(t)$.

Let $\pi_{h}u_{0}\in\mathbb{V}_{h}$ be the piecewise interpolation of $u_{0}$. We consider the following fully discrete scheme for (2):

For $n=1,2,...,N$, find $u_{n,h}\in\mathbb{V}_{h}$ such that

with initial data $u_{0,h}=\pi_{h}u^{0}$. Here $(\cdot,\cdot)_{D}$ denotes the $L^{2}$-inner product on $D$, and $a\left(u,v\right):=\left(\alpha\nabla u,\nabla v\right)_{D}+\left(cu_{n,h},v\right)_{D}$.

Introduce the numerical solution vector $U_{n}$, the mass matrix $\mathcal{M}$, the total stiffness matrix $\mathcal{S}$, and the load vector $F$, defined respectively by

Notice that $\mathcal{M}$ and $\mathcal{S}$ are both symmetric and positive definite. Since $u_{n,h}(x)=\sum\limits_{k=1}^{M}u_{n,h}(x_{k})\varphi_{k}(x),$ for a given $U_{0}$ the system (3) is of the matrix form

In this section, we give the SEAM algorithm, based on POD, and a parallel version. Let $U_{n}$ be the high-fidelity numerical solution vector defined by (4) for each $n$, and let

be the numerical solution matrix of (2). Denote

and let

be the eigenvalues of $X$, then the singular values of $U$ are

Set

Without loss of generality, we assume that $1\leq\hat{d}<<M$. The POD method is to find the standard orthogonal basis functions $\beta_{j}$ $(j=1,2,\cdots,p)$ of $\widehat{\mathcal{U}}$ such that for any $p$, $1\leq$ $p\leq\hat{d}$, the mean square error between $U_{n}$ $(n=0,1,\cdots)$ and its orthogonal projection onto the space $\text{span}\{\beta_{1},\cdots,\beta_{p}\}$ is minimized on average, i.e.

A solution $\left\{\beta_{j}\right\}_{j=1}^{p}$ to (7) is called a POD-basis of rank $p$.

We will show later (cf. Theorem 4.1 and Corollary 4.1) that the non-zero eigenvalues of $X$ except for the principal eigenvalue are close to zero if $n$ is not very large and the time step size $\tau$ is small enough. Based upon the theoretical result, we choose $p=1$ in the POD method (7). This means that we only need to use the biggest eigenvalue $\lambda_{0}(X)$ and the associated eigenvector $\mathbf{b}_{0}=(b_{0},b_{1},\cdots,b_{n})^{T}$ to obtain the rank 1 POD basis (cf. [43])

For convenience, we name the POD method with $p=1$ as Single Eigenvalue Acceleration Method (SEAM), and the corresponding SEAM scheme for the full discretization (4) is as follows:

Given $\alpha_{0}=\beta_{1}^{\top}U_{0}\in\Re,$ for $k=1,2,\cdots,n,$ the SEAM approximation solution $U^{seam}_{k}$ at $k$-th time step is obtained by

Note that the coefficients $\beta_{1}^{\top}(\mathcal{M}+\tau\mathcal{S})\beta_{1}$ and $\beta_{1}^{\top}\mathcal{M}\beta_{1}$ in (8) are positive constants.

By Theorem 4.1 the model-order-reduction error of SEAM (8) is (cf. [39])

where $C$ is a positive constant independent of $n$ and $\tau$. This means that the SEAM algorithm may be inaccurate when $n^{2}\tau$ is not small enough. In other words, we can not expect the desired accuracy if $n$ is large. To overcome such a weakness, we will give a parallel SEAM below (cf. Remark 4.3).

For simplicity of notation, we assume $N+1=(\tilde{n}+1)(n+1)$. We divide the total numerical solution matrix as

with the submatrices

Then the parallel SEAM can be described as following:

  Algorithm 1 Parallel SEAM

In this section, we apply a perturbation technique to investigate the singular value distribution of the numerical solution matrix $U$ in (5). We need the following result due to Hoffman and Wielandt [14]:

Let us turn back to the matrix equation (4). Since

we rewrite (4) as

with $\bm{\mathcal{A}}=\mathcal{M}^{-1/2}\mathcal{S}\mathcal{M}^{-1/2}$, for a given $U_{0}$ and $n=1,2,...,N$.

In what follows we assume that the time step $\tau\in(0,1)$ is small enough such that

By noticing that $\bm{\mathcal{A}}$ is SPD, this assumption is reasonable.

Let $\rho(\cdot)$ represent the spectral radius of a matrix. As (11) implies $\rho(\tau\bm{\mathcal{A}})<1$, we have

In the case of $F=0$, (10) is reduced to

Inspired by (12), for $n=0,1,2,\cdots$ we define $U_{n}^{*}$ by

and define $\bar{U}_{n}^{*}$ by

Denote

In the sequel we use $a\lesssim b$ to denote $a\leq Cb$ for simplicity, where $C$ is a generic positive constant independent of $n$ and $\tau$ and may be different at its each occurrence.

The following lemma shows that $\bar{X}_{*}$ can be viewed as a perturbation of $\widetilde{X}$ in some sense.

This section provides some numerical experiments to verify the performance of SEAM for the fully discrete scheme (10) using continuous linear finite elements. All the tests are performed on a 12th-Gen Intel 3.20 GHz Core i9 computer.