Applications of a space-time FOSLS formulation
for parabolic PDEs

We consider a reduced basis method, where in a potentially expensive offline phase a basis of (highly accurate approximations of the) solutions corresponding to a suitable finite subset of the parameter space is computed so that in the subsequent online phase the solution for arbitrary given parameters can be efficiently approximated from the spanned reduced basis space.

Compared to reduced basis methods based on time stepping and proper orthogonal decomposition (POD), e.g., [HO08, EKP11, Haa17], reduced basis methods based on simultaneous space-time discretization, e.g., [UP14, Yan14, YPU14], yield not only a dimension reduction in space but in space-time; see the comparison in [GMU17]. Consequently, in any case the resulting online phase can be expected to be more efficient.

As already mentioned, compared to other space-time discretization approaches, the FOSLS formulation has the advantage that it results in a coercive bilinear form. In this setting the strongest theoretical results for the reduced basis method are available (see, e.g., [Haa17]), and the construction of the reduced basis does not need to be accompanied by the construction of a corresponding test space that yields inf-sup stability.

Finally, the reliable and efficient a posteriori error estimator of the FOSLS formulation is beneficial for both building the reduced basis and for certification of the obtained approximate solution in the online phase.

We consider a large class of optimal control problems constrained by second-order parabolic PDEs. This includes controls of the source term in $L_{2}(I\times\Omega)$ or in $Y^{\prime}$ and of the initial data in $L_{2}(\Omega)$ as well as desired observations of the PDE solution in $L_{2}(I\times\Omega)$, its spatial gradient in $L_{2}(I\times\Omega)^{d}$, or its restriction to the final time point in $L_{2}(\Omega)$.

It is well-known that optimal control problems can be equivalently formulated as saddle-point problem; see, e.g., [Lio68]. The saddle-point formulation is a coupled system of the primal equation and some adjoint equation. When solving it iteratively, e.g., by a gradient descent method [Trö10, HPUU08], the approximate solutions are required over the whole space-time cylinder $I\times\Omega$ and not just at the final time point. This nullifies the potential advantage of time-stepping methods such as [MV07, MV08], which in the case of a forward problem only need to store the approximation at the current time step.

Instead space-time methods such as [GHZ12, LSTY21b, LSTY21a, LS21] aim to solve the saddle-point problem at once on some mesh of the space-time cylinder $I\times\Omega$. All of the mentioned works consider optimal control problems constrained by the heat equation with control of the source term and desired observation of the PDE solution itself. In contrast to time-stepping methods, all of them allow for locally adapted meshes because the involved bilinear forms are coercive or inf-sup stable. However, only [GHZ12], which reformulates the problem as a system of fourth order in space and second order in time requiring additional regularity, yields quasi-optimal approximations in the usual sense with the same norm on both sides, as for the others the norms in which the bilinear forms are continuous differ from the norms in which they are coercive or inf-sup stable.

In our approach, we consider the saddle-point problem when the parabolic PDE is replaced by the equivalent FOSLS. Coercivity of the latter implies that this saddle-point problem is for arbitrary discrete subspaces uniformly inf-sup stable and continuous w.r.t. the same norm so that quasi-optimality in the usual sense is valid. We also provide one reliable and one efficient a posteriori error estimator for our method. Other estimators have been derived in [GHZ12, LS21] for the respective methods.

Finally, we consider second-order parabolic PDEs on time-dependent domains. While this poses a technical difficulty for time-stepping methods [GS17, FR17, LO19, SG20, FS22], in principle, space-time methods [Leh15, HLZ16, LMN16, Moo18, BHT22] can be readily applied because the overall space-time domain remains fixed.

We show that this is indeed the case for the FOSLS method [FK21, GS21] verifying its well-posedness. The analysis requires minimal regularity assumptions on the mapping describing the motion of the spatial domain throughout time. However, we emphasize that our method ultimately only requires a mesh of the space-time domain. This is in contrast to the so-called arbitrary Lagrangian-Eulerian (ALE) time-stepping methods [SHD01, GS17, SG20], which solve the problem on the transformed time-independent domain and thus explicitly involve the mapping. Moreover, our method again allows for arbitrary discrete trial spaces and immediately provides an a reliable and efficient a posteriori error estimator.

In Section 2, we fix the general notation (Section 2.1) and recall the formulation of general second-order parabolic PDEs as a first-order system (Section 2.2). In Section 3, we introduce a corresponding reduced basis method for parameter-dependent problems (Section 3.1), discuss the generation of a suitable basis (Section 3.2), and provide a numerical experiment (Section 3.3). In Section 4, we show that the formulation of optimal control problems constrained by the first-order system as a saddle-point problem is uniformly stable for arbitrary trial spaces (Section 4.1), derive an optimality system of PDEs for a certain class of optimal control problems (Section 4.2), which is then used to derive a posteriori error estimators (Section 4.3), and provide two numerical experiments (Section 4.4), where we also demonstrate optimal convergence rate for piecewise linear trial functions if both the optimal control and the corresponding state are sufficiently smooth. In Section 5, we show that the formulation of general second-order PDEs on time-dependent domains as first-order system induces, as for time-independent domains, a linear isomorphism and provide two numerical experiments (Section 5.2). Finally, we draw some conclusion in Section 6.

In this work, by $C\lesssim D$ we will mean that $C\geq 0$ can be bounded by a multiple of $D\geq 0$, independently of parameters on which $C$ and $D$ may depend. Obviously, $C\gtrsim D$ is defined as $C\lesssim D$, and $C\eqsim D$ as $C\lesssim D$ and $C\gtrsim D$.

For normed linear spaces $E$ and $F$, we will denote by $\mathcal{L}(E,F)$ the normed linear space of bounded linear mappings $E\rightarrow F$. For simplicity only, we exclusively consider linear spaces over the scalar field $\mathbb{R}$.

Let $\Omega\subset\mathbb{R}^{d}$, $d\geq 1$, be a Lipschitz domain with boundary $\Gamma:=\partial\Omega$, and $T>0$ a given end time point with corresponding time interval $I:=(0,T)$. We abbreviate the space-time cylinder $Q:=I\times\Omega$ with lateral boundary $\Sigma:=I\times\Gamma$. We consider the following parabolic PDE with homogeneous Dirichlet boundary conditions

Here, we will require that ${\bf A}={\bf A}^{\top}\in L_{\infty}(Q)^{d\times d}$ is uniformly positive, ${\bf b}\in L_{\infty}(Q)^{d}$, and $c\in L_{\infty}(Q)$. For unique existence of a (weak) solution $u$ of (2.1), we recall the following theorem from, e.g., [SS09, Theorem 5.1].

Any $f\in Y^{\prime}=L_{2}(I;H^{-1}(\Omega))$ can be written as

for some $f_{1}\in L_{2}(Q)$ and ${\bf f}_{2}\in L_{2}(Q)^{d}$, where $\int_{Q}\operatorname{div}_{\bf x}{\bf f}_{2}\,v\,{\rm d}{\bf x}\,{\rm d}t:=-\int_{Q}{\bf f}_{2}\cdot\nabla_{\bf x}v\,{\rm d}{\bf x}\,{\rm d}t$ for $v\in Y$. With such a decomposition, and ${\bf u}=(u_{1},{\bf u}_{2})\colon Q\rightarrow\mathbb{R}\times\mathbb{R}^{d}$, (2.2) is equivalent to the first-order system¹¹1The system $\left(\begin{array}[]{@{}c@{}}\operatorname{div}{\bf u}-{\bf b}\cdot{\bf A}^{-1}{\bf u}_{2}+cu_{1}\\ {\bf u}_{2}+{\bf A}\nabla_{\bf x}u_{1}\\ u_{1}(0,\cdot)\end{array}\right)=\left(\begin{array}[]{@{}c@{}}f_{1}+{\bf b}\cdot{\bf A}^{-1}{\bf f}_{2}\\ -{\bf f}_{2}\\ u_{0}\end{array}\right)$ studied in [GS21] leads to (2.10) by pre-multiplication with $\left(\begin{array}[]{@{}ccc@{}}I&{\bf b}\cdot{\bf A}^{-1}&0\\ 0&-I&0\\ 0&0&I\end{array}\right)$, which is a linear isomorphism on the space $L$ defined in (2.12).

as is shown in the next theorem.

Notice that $G{\bf u}={\bf f}$ is equivalent to the variational problem

Since the bilinear form at the left-hand side is bounded, symmetric and coercive, it provides the ideal setting for the application of Galerkin discretizations. The Galerkin solution from the employed trial space is a quasi-best approximation to ${\bf u}$ w.r.t. $\|\cdot\|_{U}$. For any approximation $\widetilde{\bf u}=(\widetilde{u}_{1},\widetilde{\bf u}_{2})$ of ${\bf u}$ we have the computable a posteriori error estimator $\|{\bf u}-\widetilde{\bf u}\|_{U}\eqsim\|{\bf f}-G\widetilde{\bf u}\|_{L}$. Moreover, the following lemma shows that $\|u_{1}-\widetilde{u}_{1}\|_{X}\lesssim\|{\bf u}-\widetilde{\bf u}\|_{U}$.

Finally in this section, we provide more information about the splitting (2.3).

Given some parameter ${\bm{\mu}}$, the idea of the reduced basis method is to compute an approximation ${\bf u}^{N}[{\bm{\mu}}]$ of ${\bf u}[{\bm{\mu}}]$ from the (low-dimensional) span of some snapshots $\{{\bf u}[{\bm{\mu}}^{(1)}],\dots,{\bf u}[{\bm{\mu}}^{(N)}]\}$.

Instead of ${\bf u}[{\bm{\mu}}^{(i)}]$, very accurate approximations ${\bf u}^{\delta}[{\bm{\mu}}^{(i)}]$ thereof are computed in the offline phase. We will choose ${\bf u}^{\delta}[{\bm{\mu}}^{(i)}]$ as the best approximation to ${\bf u}[{\bm{\mu}}^{(i)}]$ w.r.t. $\|G[{\bm{\mu}}^{(i)}](\cdot)\|_{L}$ from some high-dimensional subspace $U^{\delta}\subset U$, i.e., as the solution of the Galerkin system $\langle G[{\bm{\mu}}^{(i)}]{\bf u}^{\delta}[{\bm{\mu}}^{(i)}]-{\bf f}[{\bm{\mu}}^{(i)}],G[{\bm{\mu}}^{(i)}]{\bf v}^{\delta}\rangle_{L}=0$ for all ${\bf v}^{\delta}\in U^{\delta}$. One easily checks that the parameter-separability (3.1) of the coefficients and the right-hand sides implies parameter-separability of the bilinear form

the linear form

and the squared norm

where again $\theta_{q}^{(\cdot)}:\mathcal{P}\to\mathbb{R}$, and each $b_{q}:U\times U\to\mathbb{R}$ is a continuous bilinear form and each $l_{q}:U\to\mathbb{R}$ is a continuous linear form. Besides the functions ${\bf u}^{\delta}[{\bm{\mu}}^{(i)}]$, which form the reduced basis, in the offline phase we further compute the matrices and vectors

In the online phase, we compute ${\bf u}^{N}[{\bm{\mu}}]$ as the best approximation to ${\bf u}[{\bm{\mu}}]$ w.r.t. $\|G[{\bm{\mu}}](\cdot)\|_{L}$ from the span of the reduced basis, i.e., as the solution of a low-dimensional Galerkin system. Assuming that the snapshots ${\bf u}^{\delta}[{\bm{\mu}}^{(1)}],$ $\dots$, ${\bf u}^{\delta}[{\bm{\mu}}^{(N)}]$ are linearly independent, the corresponding coefficient vector ${\bf c}^{N}[{\bm{\mu}}]$ of ${\bf u}^{N}[{\bm{\mu}}]$ with respect to this basis is just given by

We stress that the computational effort to solve this linear system depends only on $N$ and not on ${\rm dim}(U^{\delta})$ as the involved matrices and vectors have been already computed in the offline phase. Then, the quantity of interest is given by

Again, the involved vectors $F_{q}({\bf u}^{\delta}[{\bm{\mu}}^{(i)}])\big{)}_{i=1}^{N}$ have been precomputed in the offline phase. Finally, we can even instantly estimate the discretization error in the online phase by

The constants hidden in the $\eqsim$-symbol depend only on the $L_{\infty}$-norms of the coefficients ${\bf A}[{\bm{\mu}}],{\bf b}[{\bm{\mu}}],c[{\bm{\mu}}]$ and the smallest eigenvalue of ${\bf A}[{\bm{\mu}}]$.

It remains to explain how to determine a suitable reduced basis $\{{\bf u}^{\delta}[{\bm{\mu}}^{(1)}],$ $\dots$, ${\bf u}^{\delta}[{\bm{\mu}}^{(N)}]\}$. Given some sufficiently large training set $\mathcal{P}_{\rm train}\subseteq\mathcal{P}$ and some tolerance $\epsilon_{\rm tol}>0$, we employ a greedy algorithm that starting with ${\bf u}^{0}[{\bm{\mu}}]:=0$, ${\bm{\mu}}\in\mathcal{P}$, iteratively adds the snapshot

to the reduced basis and increments $N$ by one until

Note that this procedure terminates at most after $\#\mathcal{P}_{\rm train}$ steps and provides indeed a basis if the discretization error in $U^{\delta}$ is negligible. For possible choices of the training set $\mathcal{P}_{\rm train}$, we refer, e.g., to [Haa17, Remark 2.44].

We consider the example from [GMU17]: Let $\Omega=(0,1)$ and $T:=0.3$. We consider the parameter set $\mathcal{P}:=[0.5,1.5]\times[0,1]\times[0,1]\subset\mathbb{R}^{3}$ with ${\bf A}[{\bm{\mu}}]:=\mu_{1}$, ${\bf b}[{\bm{\mu}}]:=\mu_{2}$, and $c[{\bm{\mu}}]=\mu_{3}$. Moreover, we choose the right-hand sides $f_{1}$, ${\bf f}_{2}$, and $u_{0}$ independently of the parameters with

on the space-time cylinder $Q=(0,1)^{2}$, ${\bf f}_{2}:=0$, and $u_{0}(x):=\sin(2\pi x)$ on $\Omega$, which corresponds to the solution $u[(1,0.5,0.5)](t,x)=\sin(2\pi x)\cos(4\pi t)$.

We divide both $\Omega$ and $I=(0,T)$ into $2^{6}$ subintervals and choose $U^{\delta}$ to be the $2$-fold Cartesian product of continuous piecewise bi-cubic functions³³3The reason why we consider bi-cubic instead of bi-affine elements as in [GMU17] is that we measure in a stronger norm but still want the approximation error in the space $U^{\delta}$ to be negligible., with the first coordinate space restricted by homogeneous Dirichlet boundary conditions on $\Sigma$. The training set $\mathcal{P}_{\rm train}$ is chosen as $17$ equidistantly distributed points in $\mathcal{P}$ in each direction, and the tolerance $\epsilon_{\rm tol}$ is chosen as $10^{-3}$. Figure 3.1 displays the approximation error $\max_{{\bm{\mu}}\in\mathcal{P}_{\rm train}}\|G[{\bm{\mu}}]({\bf u}[{\bm{\mu}}])-G[{\bm{\mu}}]({\bf u}^{N}[{\bm{\mu}}])\|_{L}$ throughout the greedy algorithm of the offline phase with the final $N=21$, where the $y$-axis is scaled logarithmically. As expected from [BCD⁺11], we observe exponential convergence. In Figure 3.2, we plot the discretization error $\|G[{\bm{\mu}}]({\bf u}[{\bm{\mu}}])-G[{\bm{\mu}}]({\bf u}^{N}[{\bm{\mu}}])\|_{L}$ of the reduced basis method applied for ${\bm{\mu}}$ in the test sets $\big{\{}(\mu_{1},0,0)\,:\,\mu_{1}\in[0.5,1.5]\big{\}}$ and $\big{\{}(0.5,\mu_{2},0.75)\,:\,\mu_{2}\in[0,1]\big{\}}$. For comparison, we also plot the best possible error $\|G[{\bm{\mu}}]({\bf u}[{\bm{\mu}}])-G[{\bm{\mu}}]({\bf u}^{\delta}[{\bm{\mu}}])\|_{L}$ that one could hope for with the reduced basis method. Although the greedy algorithm only guarantees that these errors are below $\epsilon_{\rm tol}=10^{-3}$ on the training set $\mathcal{P}_{\rm train}$, this bound appears to hold also true on the considered test sets.

While a fair comparison to the space-time results of [GMU17] is hard, we only mention that the considered errors, although measured in the stronger norm $\|G(\cdot)\|_{L}\simeq\|\cdot\|_{U}$, are of similar magnitude as the errors of [GMU17] measured in (an approximation of) the norm $\|\cdot\|_{X}$. The number of reduced basis functions to achieve an accuracy of $\epsilon_{\rm tol}=10^{-3}$ in [GMU17] is given by $N=12$. It is also notable that, in contrast to the estimator of [GMU17], which is based on the residual corresponding to (the first components of) ${\bf u}^{\delta}[{\bm{\mu}}]-{\bf u}^{N}[{\bm{\mu}}]$, the estimator considered here is provably equivalent to the actual norm of the error $\|{\bf u}[{\bm{\mu}}]-{\bf u}^{N}[{\bm{\mu}}]\|_{U}$ and can also be efficiently computed in the online phase.

Writing (2.10) with right-hand side ${\bf f}={\bf f}^{\star}+{\bf z}$ as $\langle G{\bf u},G{\bf v}\rangle_{L}=\langle{\bf f}^{\star}+{\bf z},G{\bf v}\rangle_{L}$ for all ${\bf v}\in U$, and introducing the bilinear forms

and the linear forms

for ${\bf u},{\bf v}\in U$ and ${\bf z},{\bf y}\in Z$, and noting that the constant term $\|w^{\star}\|_{W}^{2}$ in $\|F{\bf u}-w^{\star}\|_{W}^{2}$ can be neglected for minimization, the optimal control problem (4.1) can be rewritten as

The following lemma implies well-posedness and equivalence to a saddle-point problem.