Reduced-order modeling using the frequency-domain method for parabolic partial differential equations

1.   Introduction

We consider the following parabolic problem to be the target of reduced-order modeling applied in the frequency-domain method. For a given $u_{0}\in H_{0}^{1}(\Omega)$,

where $\Omega$ represents an open convex polygon in $\mathbb{R}^{2}$; $\kappa\in L^{2}(\Omega)$ and $\sigma\in W^{1, \infty}(\Omega)$ denote positive functions of $x\in\Omega$ satisfying $\kappa_{*}\leqslant\kappa\leqslant\kappa^{*}$, $\sigma_{*}\leqslant\sigma\leqslant\sigma^{*}$, and $|\nabla\sigma|\leqslant\sigma^{*}$ with positive constants $\kappa_{*}$, $\kappa^{*}$, $\sigma_{*}$, and $\sigma^{*}$; $f(\cdot, t) \in L^{2}(\Omega)$; and $\partial\Omega$ denotes the boundary of $\Omega$.

Recently, many articles have addressed reduced-order modeling to determine low-order models of dynamical systems, such as very complex turbulence flows and problems of optimization or feedback control problems (see, e.g., ^[2,3]). As a tool for deriving low-order models for the given problems, many researchers have used the proper orthogonal decomposition method (POD) (see, e.g., ^[7,8,9,16]). The combination of the isogeometric analysis and POD was investigated for parabolic problems in ^[20] and unsteady convection-dominated convection-diffusion-reaction problems in ^[15]. The POD method provides a reduced-order basis for the modal decomposition of an ensemble of functions, such as data obtained during the course of experiments or numerical simulations. For example, suppose a finite series exists with a time step $\Delta T$ of finite element numerical solutions, so-called snapshots, of a time-dependent partial differential equation, in which the solutions are approximated using the general nodal basis of a high-dimensional finite element space in a Hilbert space. The space spanned by the snapshots is called the snapshot space. Then, using the Galerkin POD of the snapshot space, an appropriate low-order orthonormal basis, so-called a POD basis, can be employed for a low-dimensional subspace of the snapshot space. Such an orthonormal basis can be easily computed using the spectral decomposition of the correlation matrix of the snapshots in the Hilbert space. Note that the number of POD basis functions is much less than the dimension of the snapshot space in general. Once a low-order POD basis is determined, we can quickly compute approximate solutions of the time-dependent partial differential equation with time step $\Delta t$ that are much less than $\Delta T$. Thus, we employ the reduced-order modeling of Galerkin POD.

In this paper, by applying the frequency-domain method to the time-dependent parabolic equation (1.1), we provide reduced-order modeling of the Galerkin POD to determine approximate solutions of frequency-dependent elliptic equations quickly. We first transform a parabolic equation to the frequency-variable elliptic equations using the Fourier integral transform in time. Such a frequency-domain method enables easily implementing a parallel computation algorithm to approximate the frequency-variable solutions because the frequency-variable elliptic equations have independent frequencies (see, e.g., ^{[4,5,6,10,11,17]}; see ^{[12,13,18,19]} for the case of using the Laplace transformation).

Applying the Fourier transformation for the space-time problem (1.1), we have the following set of complex-valued elliptic equations depending on the frequency $\omega$: for all $\omega\in\mathbb{R}$,

where $f$ and $u$ are extended by zero for $t < 0$ and $t > T$ for the Fourier transformation. The Fourier transform $\hat{u}(\cdot, \omega)$ of a function $u(\cdot, t)$ in time and the Fourier inversion are given by

This paper investigates the combination of the frequency-domain method and the reduced-order modeling. We apply Galerkin POD method to determine approximate solutions of the frequency-variable elliptic equations (1.2) instead of the time-dependent parabolic equation (1.1). A set of snapshots consists of the finite element solutions of the frequency-variable elliptic equations with some sampled frequencies, in which the solutions are approximated using the general high-dimensional nodal basis of the finite element space. Then, from the spectral decomposition of the correlation matrix of the snapshots, we determine a low-order Galerkin POD basis for a subspace of the snapshot space spanned by the snapshots. Using a low-order Galerkin POD basis, we compute approximate solutions (POD-solutions) of the frequency-variable elliptic equations for sufficiently many frequencies to determine accurate inverse Fourier transforms for the solutions in the time variable. We use the Gaussian quadrature rule based on Legendre–Gauss–Lobatto (LGL) points for the accurate numerical integration of the inverse Fourier transformation. Thus, the selected frequencies for POD-solutions are sufficiently many LGL points on an appropriate interval. The number of sample snapshots must be much less than the number of the POD-solutions to reduce the total computational cost because the snapshots are approximated using full-dimensional basis functions, but the POD solutions are computed using low-order POD basis functions. Regarding numerical computation, a fast-solving parallel computation can be run to determine the snapshots to reduce the total computational time, which is merit of the Galerkin POD method applied to the frequency-variable equations (see, e.g., ^[7,8,14]).

The paper is organized as follows. Section 2 provides an overview of the Galerkin POD, and Section 3 presents reduced-order modeling for the frequency-domain method to approximate the parabolic equation. Finally, Sections 4 presents some numerical experiments.

2.   Galerkin proper orthogonal decomposition

The POD of order $\ell$ is to determine a set of ordered orthonormal basis functions, such that the snapshots can be expressed optimally using the selected first $\ell$ basis functions, where $\ell$ is a positive integer ^[7,8,16]. We briefly review the Galerkin POD in the context of the finite element method. Let $X_h$ be a finite-dimensional subspace of a given Hilbert space $X$ endowed with the inner product $(\cdot, \cdot)_{X}$ and norm $\| \, \cdot\, \|_X$, and let $\{\, \phi_p \, \}_{p = 1}^n$ be a basis of the space $X_h$. For example, we consider a nodal basis $\{\, \phi_p \, \}_{p = 1}^n$ for a Galerkin finite element subspace $X_h$ consisting of piecewise linear functions of the Sobolev space $X = H^1(\Omega)$, where $\Omega$ is a given domain. For a set of snapshots $S = \{ \, y_{1}, \, \cdots, \, y_{m} \, \} \subset X_h$, we define a snapshot subspace

and consider the orthonormal basis $\{\psi_{p}\}^{d}_{p = 1}$ of $X_{S}$, where $d = \text{dim}\, \, X_{S}$ is the dimension of $X_S$. Then, the method of Galerkin POD of order $\ell \leq d$ consists of choosing the orthonormal basis such that the mean square $X$-norm error between the snapshots $y_q$ and the corresponding $\ell$-th partial sums is minimized:

Using the orthogonality of $\psi_p$'s yields that

Hence, the minimization problem (2.1) is equivalent to the following maximization problem:

Let $\hat{X} = \big(\, (\phi_p, \phi_q)_X \, \big) \in {\mathbb{R}}^{n\times n}$ be the positive definite finite element matrix. The matrices containing the coefficients of $y_q$ and $\psi_q$ are denoted by $Y \in \mathbb{R}^{n \times m}$ and $\Psi \in \mathbb{R}^{n \times d}$, respectively, in the expansion with respect to the basis functions $\phi_p$, that is,

In addition, each snapshot can be expressed as a linear combination of the orthonormal basis functions such that

Using (2.3) and (2.4) yields

and

so that we obtain the following identities

and

Next, $Q^\ell$ denotes the first $\ell$ columns and $Q_\ell$ the first $\ell$ rows of a given matrix $Q$. Then, the problem (2.2) is equivalent to the following problem of in matrix form:

where $\| \, \cdot \, \|_F$ denotes the Frobenius norm of the matrix and $I_\ell$ denotes the identity $\ell \times \ell$ matrix.

Let $\hat Y = \hat X^{1/2}\, Y$, let $U$ and $V$ be the left and right singular vectors, respectively, in the singular value decomposition (SVD) of $\hat Y$:

For (2.7), using the Fritz John necessary conditions or Karush-Kuhn-Tucker conditions for the optimality of $\Psi^\ell$ (see ^[1] for details) yields

and the optimal solution for (2.7) is given by

Hence, the optimal solution for the problem (2.1) is given by

In contrast, from the singular value decomposition of $\hat Y$ that

Thus, the optimal solution for the POD-basis problem (2.1) can be more easily given by using the right singular vectors of $\hat Y$:

The error between the snapshots and their POD solutions is given by (see ^[8,16])

Let $K$ be the correlation matrix corresponding to the snapshots $\{y_q\}_{q = 1}^m$ in the Hilbert space $X$:

Then, the right singular vectors $V$ holds the following spectral decomposition of $K$:

In the context of the POD-basis approach for the finite element methods, the eigenvalue problem for $m\times m$ matrix $\hat Y^T \, \hat Y$ is more practical to solve than the eigenvalue problem for $n\times n$ matrix $\hat Y \, \hat Y^T$ in cases where the size of the input collection $m$ is significantly smaller than the number of coefficients $n$ needed to represent each function for the general basis functions of finite element space.

Finally, we report on reduced-order modeling for the finite element method. If we have the following linear system of a Galerkin finite element discretization in the space $X_h \subset X$:

then using the linear transformation $\mathbf{{y}} = \Psi^\ell \, \mathbf{{x}} \in {\mathbb{R}}^n$ for $\mathbf{{x}} \in {\mathbb{R}}^\ell$, we obtain the following reduced-order modeling of order $\ell$ to approximate the solution for the above linear system:

3.   Reduced-order modeling for the frequency-domain method

3.1.   Finite element approximation based on the frequency-domain method

This section provides reduced-order modeling using the Galerkin POD for the following complex-valued elliptic equations depending on $\omega$: for all $\omega\in\mathbb{R}$

This paper applies the standard notation and definitions for the real-valued Sobolev spaces $H^{s}(\Omega)$, associated with the scalar product $(\cdot, \, \cdot)_s$ and norm $\|\cdot\|_{s}$, $s\geq 0$. Nevertheless, $H^{0}(\Omega)$ coincides with $L^{2}(\Omega)$, in which the associated inner product and norm are denoted by $(\cdot, \cdot)$ and $\|\cdot\|$, respectively. The real and imaginary parts of a complex-valued vector or scalar function $\hat v$ are denoted by $\hat v_{r}$ and $\hat v_{i}$ respectively. Then, the $L^{2}({\Omega})$ inner product and norm for complex-valued functions $\hat u = \hat u_{r}+i \hat u_{i}$ and $\hat v = \hat v_{r}+i \hat v_{i}$ are given by

From now on, we denote by $H^{s}_{c}(\Omega): = H^{s}(\Omega)\times H^{s}(\Omega)$, $L^{2}_{c}(\Omega): = L^{2}(\Omega)\times L^{2}(\Omega)$, and $H^{1}_{0, c}(\Omega): = H^{1}_{0}(\Omega)\times H^{1}_{0}(\Omega)$ where $H^1_0({\Omega})$ is the subspace of $H^1_0({\Omega})$ vanishing on the boundary of ${\Omega}$. We identify $H^{1}_{0, c}(\Omega)$ with $V$ and define the sesquilinear form $a_{\omega}(\cdot, \cdot):V\times V\rightarrow \mathbb{C}$ for ${\omega} \in \mathbb{R}$ as

Then, the variational formulation of the equation (3.1) can determine ${\hat u}(\cdot, \omega)\in V$ such that

For the finite element approximation of (3.1), Let $T_{h}$ be a quasi-regular partition of $\Omega$ into triangles with a diameter bounded by $h < 1$. We take a standard finite element subspace $V_{h}\subset V$ such that

where the positive constant $C$ is independent of $h$ and $v$, and $|\cdot|_{1, c}$ and $|\cdot|_{2, c}$ denote the seminorms of $H^1_c({\Omega})$ and $H^2_c({\Omega})$, respectively. In this paper, we take the standard piecewise linear finite element space for the space $V_h$. We assume that $V_h = P_h \oplus i \, P_h$ where $P_h$ is the standard piecewise linear real value finite element subspace of $H^1_0({\Omega})$ over $T_h$ with $\dim(P_h) = n$ and $\{ \, \phi_j \, \}_{j = 1}^n$ is the nodal basis of $P_h$. Then, the Galerkin finite element approximation is to find $\hat{u}_{h}(\omega) = \hat{u}_{h}(\cdot, \omega) \in V_{h}$ such that

In ^[10], the authors provided the existence of the solution in Theorem 2.1, the stability in Theorem 2.2 and the error estimation in Theorem 3.1. The error estimation is given as follows. A generic positive constant denoted by $C$ may differ from place to place.

Theorem 3.1. Let $\hat{u}(\omega)$ be the approximate solutions of the Eq (3.2) and let $\hat{u}_{h}(\omega)$ be the approximate solutions of the Eq (3.3). Then the following estimations hold:

and

For the approximate Fourier inversion of frequency-variable solutions to the time variable solutions, we apply the Gaussian quadrature rule based on the LGL-points on an appropriate interval $[0, \, {\omega}^*]$ with a sufficiently large ${\omega}^* > 0$ such that $\hat u({\omega}) = \hat u(\cdot, {\omega})$ is negligible for $|{\omega}| > {\omega}^*$. Let $G_{{\omega}^*} = \{{\omega}_{j}\}^{N_{\omega}}_{j = 1}$ be the set of LGL-points on the interval $[0, \omega^{*}]$. Then, the time variable approximate solution $u_{h}(x, t)$ for the real-valued solution $u(x, t)$ of the problem (1.1) is approximated by

where $w_j$ denote the Gaussian quadrature weights corresponding the LGL-points ${\omega}_j$. The error estimation between the approximate solution $u_{h, \Delta {\omega}}$ and $u(x, t)$ is given in ^[10] where the approximate Fourier inversion $u_{h, \Delta {\omega}}$ is given by using the composite mid-point rule. This paper applies the Gaussian quadrature rule for the approximate Fourier inversion to reduce the computational cost. We do not provide an error estimation but we focus on the performance of reduced-order modeling using Galerkin POD. The error estimation can be proved following the similar arguments given in ^[10].

3.2.   Reduced-order modeling for the frequency-domain method

To determine a Galerkin POD basis, we must construct a set of snapshots consisting of Galerkin finite element solutions $\hat u_h(\xi_p) \in V_h$ for the problem (3.3) with some selected frequencies, for example, $\{\, \xi_p \, \}_{p = 1}^{N_s}$, where $N_s$ denotes the number of snapshots. For the efficiency of reduced-order modeling of the frequency-domain problem, the number of samples $N_s$ must be less than the number of approximate solutions $N_{\omega}$ to be used for the Fourier inversion. The sets of the real and imaginary parts of the snapshots are denoted as

We define two correlation matrices for two Galerkin POD-bases:

and

This paper employs the $L^2({\Omega})$-inner product for the correlation matrices. The $H^1({\Omega})$-inner product could be used, but we obtained similar results in the experiments. In addition the $\ell^2$-Euclidean inner product for the coefficient vectors of the snapshots could also be applied (see ^[8,9,17]).

We let $M_h \in \mathbb{R}^{n \times n}$ be the mass matrix over the finite element space $P_h$ such that $M_h(p, q) = \big(\, \phi_q, \, \phi_p \, \big)_{L^2({\Omega})}$, and $Y_R \in \mathbb{R}^{n \times N_s}$ and $Y_I \in \mathbb{R}^{n \times N_s}$ represent the matrices containing the coefficients of $Re(\hat{u}_{h}(\xi_{p}))$ and $Im(\hat{u}_{h}(\xi_{p}))$ with respect to the basis functions $\phi_p$ of $P_h$, respectively. Following the argument in (2.11), suppose that we have the following spectral decompositions for the correlation matrices $K_R$ and $K_I$:

where $d_R = \dim(X_R)$ and $d_I = \dim(X_I)$, and $D_R = {\rm diag}(\, {\sigma}_1^R, \cdots, {\sigma}_{d_R}^R \, )$ and $D_I = {\rm diag}(\, {\sigma}_1^I, \cdots, {\sigma}_{d_R}^I \, )$. Then, the Galerkin POD bases of order $\ell$ for $X^R$ and $X^I$ are given by, for $q = 1, \cdots, \ell$,

where $\ell \leq \min\{ \, d_R, \, d_I \, \}$, and

Now we have the reduced-order subspace $V_h^\ell$ of $V_h$ of order $\ell$ such that

where

For any ${\omega} \in {\mathbb{R}}$, a function $\hat u_h({\omega}) \in V_h = P_h \oplus i \, P_h$ can be represented by

Then, from the problem (3.3) we have the following linear system of order $2n$

where $\mathbf{{A}}_R({\omega}), \; \mathbf{{A}}_I({\omega}) \in {\mathbb{R}}^{2n \times n}$ are given by

and $\mathbf{{f}}^R, \; \mathbf{{f}}^I \in {\mathbb{R}}^{n \times 1}$ are given by

Following the arguments of reduced-order modeling of order $\ell$ given in (2.12), we obtain the following reduced system of order $2 \, \ell$:

where

and

Using the reduced-order system (3.12), we can quickly compute POD approximate solutions $\hat u_h^\ell(x, {\omega}_j)$ for the LGL-points $\{{\omega}_{j}\}^{N_{\omega}}_{j = 1}$, in which the approximate solutions are used for the Fourier inversion $u_{h}^\ell(x, t)$, the POD solution of order $\ell$, like the approximate formulation (3.6):

4.   Numerical experiments

This section presents computational experiments with three examples focused on the error discretization and elapsed time for the efficient performances of using the POD-basis. For the measure of approximate errors, we define several error types:

where $u_h(x, t)$ and $u_h^\ell(x, t)$ denote the approximate solution using the general full-basis of $V_h$ and the approximate POD-solution using the Galerkin POD-basis of order $\ell$, respectively, and the errors are approximately computed where $T = 1$ and $\Delta t = 1/100$ are used. The order $\ell$ of the POD-basis is selected from the error estimation between the snapshots and their POD-approximations given in (2.10). The marks 'Time(Full)' and 'Time(POD)' indicate the total CPU elapsed times to find approximate solutions using the full and POD bases, respectively.

For every example, we simply take $\Omega = (0, 1)^{2}$, and we set $\kappa = 1$ and $\sigma = 1$ in the problem (1.1):

All simulations are computed by using the Matlab program on a 3 GHz Intel Xeon dual-core 64-bit CPU processor with 4.00 GB RAM, without parallel computation, but parallel computation can be used for large-scale real-world problems.

Example 1. For the first test problem, we take an elementary separable example:

In this example, taking $\omega^{*} = 20$ and $N_{\omega} = 12$ is sufficient to illustrate good performance in the approximate results. For snapshots, we set $N_s = 2$ (i.e., two snapshots) with two frequencies ($\xi_{j} = (j-\frac{1}{2})\Delta{\omega}$, $j = 1, 2$ with $\Delta\omega = \omega^{*}/2$). From the snapshot space, we employed only one POD-basis (i.e., $\ell = 1$). Although the order is too small, the solution is very well approximated because of the simplicity of the problem (see Table 1). The error discretization is similar for the two cases, but the elapsed CPU times are very different, and that of the POD basis is much less than that of full basis. The outcome is because $N_{\omega} = 12$ LGL frequency approximate solutions using the full-basis are needed for the full-basis case. But two frequency snapshots using the full-basis and $N_{\omega} = 12$ LGL frequency approximate solutions of POD basis solution are needed for the POD basis case, where the cost of the spectral decomposition of the correlation matrix must be added, but it is relatively cheap.

Table 2 compares the numerical results for the trapezoidal midpoint quadrature rule and Gaussian quadrature rule for the Fourier inversion. Although we employed 240 trapezoidal midpoints, the performance of this case is much worse than using the Gaussian quadrature rule on the error discretization and elapsed CPU time. Thus, we apply the Gaussian quadrature rule for the Fourier inversion.

Example 2. In this example, we take the exact solution that has two separable components.

In this example, we set ${\omega}^* = 20$ and $N_{\omega} = 20$. For the snapshots, we set $N_s = 4$ with four frequencies $\xi_{j} = (j-\frac{1}{2})\Delta{\omega}$, $j = 1, \cdots, 4$ where $\Delta{\omega} = {\omega}^{*}/4$. We employ POD basis of order $\ell = 2$. Table 3 lists the numerical results, which are very similar to those in Example 1.

Example 3. In this example, we take an exact solution that is not separable.

In this example, we set ${\omega}^* = 25$ and $N_{\omega} = 24$. For the snapshots, we set $N_s = 20$ with four frequencies: $\xi_{j} = (j-\frac{1}{2})\Delta{\omega}$, $j = 1, \cdots, 20$ where $\Delta{\omega} = {\omega}^{*}/20$. We take POD basis of order $\ell = 15$. Table 4 presents the numerical results, which are very similar to those in Examples 1 and 2. However, we require more POD-basis functions than in the cases of Examples 1 and 2 to achieve better results.

5.   Conclusions

We investigate reduced-order modeling using Galerkin POD to determine approximate solutions of time-dependent parabolic problems. We apply the frequency-domain method to the time-domain problem using the Fourier transformation for a more efficient and fast approximation. This approach of the frequency-domain method enables easily and efficiently implementing the parallel computation. Hence, we transformed the time-dependent parabolic problem into a frequency-dependent elliptic problems. However, these elliptic problems are independent. Next, we provided reduced-order modeling using Galerkin POD to determine a very small-dimensional subspace with the orthonormal POD basis. Using such a POD basis, we quickly approximated the elliptic problems for selected frequencies, in which the frequencies are LGL-points on a given interval.

The approximate solution of the time-dependent parabolic problem was approximated by the inverse Fourier transformation using the fast Gaussian quadrature rule. In the experiments, we demonstrated that the proposed reduced-order modeling performs very well in error discretization and reducing computational costs.

Acknowledgments

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. 2023R1A2C1006588).

Conflict of interest

The authors declare that there are no conflicts of interest.