# Dynamical model reduction method for solving parameter-dependent dynamical systems

Marie Billaud-Friess, Anthony Nouy

► To cite this version:

Marie Billaud-Friess, Anthony Nouy. Dynamical model reduction method for solving parameter-dependent dynamical systems. SIAM Journal on Scientific Computing, 2017, 39 (4), pp.A1766-A1792. 10.1137/16M1071493. hal-01387231

HAL Id: hal-01387231

https://hal.science/hal-01387231

Submitted on 5 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamical model reduction method for solving parameter-dependent dynamical systems *

Marie BILLAUD-FRIESS‡and Anthony NOUY†

December 2, 2016

## Abstract

We propose a projection-based model order reduction method for the solution of parameter-dependent dynamical systems. The proposed method relies on the construction of time-dependent reduced spaces generated from evaluations of the solution of the full-order model at some selected parameters values. The approximation obtained by Galerkin projection is the solution of a reduced dynamical system with a modified flux which takes into account the time dependency of the reduced spaces. An error estimate is derived and a greedy algorithm using this error estimate is proposed for the adaptive selection of parameters values. The resulting method can be interpreted as a dynamical low-rank approximation method with a subspace point of view and a uniform control of the error over the parameter set.

## 1 Introduction

Parameter-dependent equations are considered in many problems of scientific computing such as optimization, control or uncertainty quantification. For complex numerical models, model order reduction methods are usually required for an efficient estimation of the solution for many values of the parameters (multi-query context). Classical model order reduction methods for parameter-dependent equations are the Reduced Basis (RB) method [20], the Proper Orthogonal Decomposition (POD) method [33] or the Proper Generalized Decomposition (PGD) method [26]. These methods can be interpreted as low-rank approximation methods with different constructions of the approximation for different controls of the error over the parameter set (uniform control for RB or control in mean-square sense for POD and PGD, see, e.g., [27]).

<footnote>*This work was supported by the french research group GDR MoMaS.</footnote>

<footnote>†Ecole Centrale de Nantes, GeM, UMR CNRS 6183, France.</footnote>

<footnote>‡Corresponding author (marie.billaud-friess@ec-nantes.fr).</footnote>
# Dynamical model reduction method for solving parameter-dependent dynamical systems *

Marie BILLAUD-FRIESS‡and Anthony NOUY†

December 2, 2016

## Abstract

We propose a projection-based model order reduction method for the solution of parameter-dependent dynamical systems. The proposed method relies on the construction of time-dependent reduced spaces generated from evaluations of the solution of the full-order model at some selected parameters values. The approximation obtained by Galerkin projection is the solution of a reduced dynamical system with a modified flux which takes into account the time dependency of the reduced spaces. An a posteriori error estimate is derived and a greedy algorithm using this error estimate is proposed for the adaptive selection of parameters values. The resulting method can be interpreted as a dynamical low-rank approximation method with a subspace point of view and a uniform control of the error over the parameter set.

## 1 Introduction

Parameter-dependent equations are considered in many problems of scientific computing such as optimization, control or uncertainty quantification. For complex numerical models, model order reduction methods are usually required for an efficient estimation of the solution for many values of the parameters (multi-query context). Classical model order reduction methods for parameter-dependent equations are the Reduced Basis (RB) method [20], the Proper Orthogonal Decomposition (POD) method [33] or the Proper Generalized Decomposition (PGD) method [26]. These methods can be interpreted as low-rank approximation methods with different constructions of the approximation for different controls of the error over the parameter set (uniform control for RB or control in mean-square sense for POD and PGD, see, e.g., [27]).

<footnote>*This work was supported by the french research group GDR MoMaS.</footnote>

<footnote>†Ecole Centrale de Nantes, GeM, UMR CNRS 6183, France.</footnote>

<footnote>‡Corresponding author (marie.billaud-friess@ec-nantes.fr).</footnote>
This paper is concerned with the solution of parameter-dependent non-autonomous dynamical systems of the form

$$ \left\{ \begin{array}{l} u'(t, \xi) = f(u(t, \xi), t, \xi), \\ u(0, \xi) = u^0(\xi), \end{array} \right. \qquad (1) $$

where the flux $f$ and initial condition depend on some parameters $\xi$ with values in a parameter set $\Xi$. The solution $u(t, \xi)$ belongs to the high-dimensional state space $X = \mathbb{R}^d$. Regarding time-dependent problems, different model order reduction (MOR) methods have been considered in the literature. In the context of RB methods, an approximation is obtained by a (Petrov-)Galerkin projection of the solution onto a time-independent low-dimensional space $X_r$ (so-called reduced space) of $X$, which results in an approximation of the form

$$ u(t, \xi) \approx \sum_{i=1}^{r} v_i \alpha_i(t, \xi), \qquad (2) $$

where $\{v_i\}_{i=1}^r$ constitutes a basis of $X_r$. Different methods have been proposed for the construction of time-independent reduced spaces $X_r$ (see, e.g., [7, 2] for a review of such methods). In [16], $X_r$ is obtained as the span of snapshots $u(t^k, \xi^i)$ (in both time and parameter) of the solution of the full-order model. However, for a high-dimensional state space $X$, it is not feasible (and far from optimal) to retain a large number of snapshots. Then, one can rely on a POD of the snapshots matrix in order to extract subspaces which are optimal in a mean-square sense. A more popular approach has been considered in [14, 17, 18, 19] for linear evolution problems, with an adaptive construction of the reduced space by a POD-greedy algorithm using a posteriori error estimates. At each iteration of this algorithm, the reduced space is enriched by the dominant modes of the POD of a trajectory $u(\cdot, \xi^i)$ of the full-order model, where $\xi^i$ maximizes over the parameters set an estimate of the current approximation error. This strategy has also been considered in [28, 12, 21, 34] for the solution of nonlinear problems, including nonlinear dynamical systems. In [34], it was combined with a discrete variant of the Empirical Interpolation Method [1, 8] for the approximation of nonlinear terms and for the efficient evaluation of the error estimate which is required in the greedy algorithm.

PGD method has been considered in [25, 26] for the low-rank approximation of the solution of stochastic evolution equations, with an approximation of the form

$$ u(t, \xi) \approx \sum_{i=1}^{r} v_i(t) \alpha_i(\xi), \qquad (3) $$

which is seen as a rank-$r$ element in the tensor space $X^{[0,T]} \otimes \mathbb{R}^\Xi$. This approach adopts a variational approach in time. The resulting approximation can be seen as a projection of $u(t, \xi) \in X$ onto a time-dependent reduced space

$$ X_r(t) = \operatorname{span}\{v_1(t), \dots, v_r(t)\}, $$
This paper is concerned with the solution of parameter-dependent non-autonomous dynamical systems of the form

$$ \left\{ \begin{array}{l} u'(t, \xi) = f(u(t, \xi), t, \xi), \\ u(0, \xi) = u^0(\xi), \end{array} \right. \qquad (1) $$

where the flux $f$ and initial condition depend on some parameters $\xi$ with values in a parameter set $\Xi$. The solution $u(t, \xi)$ belongs to the high-dimensional state space $X = \mathbb{R}^d$. Regarding time-dependent problems, different model order reduction (MOR) methods have been considered in the literature. In the context of RB methods, an approximation is obtained by a (Petrov-)Galerkin projection of the solution onto a time-independent low-dimensional space $X_r$ (so-called reduced space) of $X$, which results in an approximation of the form

$$ u(t, \xi) \approx \sum_{i=1}^{r} v_i \alpha_i(t, \xi), \qquad (2) $$

where $\{v_i\}_{i=1}^r$ constitutes a basis of $X_r$. Different methods have been proposed for the construction of time-independent reduced spaces $X_r$ (see, e.g., [7, 2] for a review of such methods). In [16], $X_r$ is obtained as the span of snapshots $u(t^k, \xi^i)$ (in both time and parameter) of the solution of the full-order model. However, for a high-dimensional state space $X$, it is not feasible (and far from optimal) to retain a large number of snapshots. Then, one can rely on a POD of the snapshots matrix in order to extract subspaces which are optimal in a mean-square sense. A more popular approach has been considered in [14, 17, 18, 19] for linear evolution problems, with an adaptive construction of the reduced space by a POD-greedy algorithm using a posteriori error estimates. At each iteration of this algorithm, the reduced space is enriched by the dominant modes of the POD of a trajectory $u(\cdot, \xi^i)$ of the full-order model, where $\xi^i$ maximizes over the parameters set an estimate of the current approximation error. This strategy has also been considered in [28, 12, 21, 34] for the solution of nonlinear problems, including nonlinear dynamical systems. In [34], it was combined with a discrete variant of the Empirical Interpolation Method [1, 8] for the approximation of nonlinear terms and for the efficient evaluation of the error estimate which is required in the greedy algorithm.

PGD method has been considered in [25, 26] for the low-rank approximation of the solution of stochastic evolution equations, with an approximation of the form

$$ u(t, \xi) \approx \sum_{i=1}^{r} v_i(t) \alpha_i(\xi), \qquad (3) $$

which is seen as a rank-$r$ element in the tensor space $X^{[0,T]} \otimes \mathbb{R}^\Xi$. This approach adopts a variational approach in time. The resulting approximation can be seen as a projection of $u(t, \xi) \in X$ onto a time-dependent reduced space

$$ X_r(t) = \operatorname{span}\{v_1(t), \dots, v_r(t)\}, $$
which allows to well capture transient phenomena. However, the projection is obtained with a global in time variational principle and is not optimal at each instant $t$. In the context of parameter dependent equations, reduced basis methods based on Petrov-Galerkin space-time (PG-ST) formulations have also been introduced, see for example [31, 32]. Such an approach provides a low-rank approximation of the form (3). At the discrete level, it differs from usual reduced basis approaches, introduced e.g. in [17], since they do not rely on time-stepping scheme except in very particular cases. For a detailed comparison with standard RB method, see the recent review [15].

Dynamical low-rank approximation methods have been considered in [9, 10, 23, 24, 30], with different types of constructions but an approximation of the same form

$$u(t, \xi) \approx \sum_{i=1}^{r} v_i(t) \alpha_i(t, \xi), \quad (4)$$

which at each instant $t$ can be interpreted as a rank-$r$ approximation in the tensor space $X \otimes \mathbb{R}^{\Xi}$. Again, the approximation can be seen as a projection onto a time-dependent reduced space $X_r(t)$ but here, the projection is obtained through principles which are local in time (e.g., Dirac-Frenkel principle). This results in a reduced order model which takes the form of a dynamical system.

In this paper, we propose a new model order reduction method for solving general parameter-dependent dynamical systems of the form (1). This method provides an approximation of the form (4), so that it can be interpreted as a dynamical low-rank approximation method. However, the proposed method differs from existing dynamical low-rank approximation methods in that it adopts a subspace point of view and provides a uniform control of the error over the parameter set. The reduced space $X_r(t) \subset X$ is here obtained as the span of some selected trajectories $\{u(\cdot, \xi^i)\}_{i=1}^r$ of the full-order model (1). The approximation is obtained by solving a reduced dynamical system of size $r$ obtained by Galerkin projection (see, e.g., [4]), which has the form of the initial dynamical system with a modified flux which takes into account the time dependency of the subspace. The resulting approximation (when discarding some necessary numerical approximations) interpolates the solution map $\xi \mapsto u(\cdot, \xi)$ at points $\{\xi^1, \dots, \xi^r\}$. An a posteriori error estimate $\Delta_r(t, \xi)$ (local in time) is derived in the lines of [34] using the logarithmic Lipschitz constant associated to the flux. This error estimate is used in a greedy procedure for the adaptive selection of interpolation points, where at step $r$, the next interpolation point corresponds to a maximizer over the parameters set $\Xi$ of a certain norm of $t \mapsto \Delta_r(t, \xi)$.

The paper is organized as follows. In Lemma 2, we introduce the Galerkin method for the projection of the dynamical system (1) onto a time-dependent reduced space and we derive an a posteriori error estimate. In Lemma 3, we present strategies for the construction of reduced spaces $X_r$, including the classical POD-greedy strategy for the construction of
which allows to well capture transient phenomena. However, the projection is obtained with a global in time variational principle and is not optimal at each instant $t$. In the context of parameter dependent equations, reduced basis methods based on Petrov-Galerkin space-time (PG-ST) formulations have also been introduced, see for example [31, 32]. Such an approach provides a low-rank approximation of the form (3). At the discrete level, it differs from usual reduced basis approaches, introduced e.g. in [17], since they do not rely on time-stepping scheme except in very particular cases. For a detailed comparison with standard RB method, see the recent review [15].

Dynamical low-rank approximation methods have been considered in [9, 10, 23, 24, 30], with different types of constructions but an approximation of the same form

$$u(t, \xi) \approx \sum_{i=1}^{r} v_i(t) \alpha_i(t, \xi), \quad (4)$$

which at each instant $t$ can be interpreted as a rank-$r$ approximation in the tensor space $X \otimes \mathbb{R}^{\Xi}$. Again, the approximation can be seen as a projection onto a time-dependent reduced space $X_r(t)$ but here, the projection is obtained through principles which are local in time (e.g., Dirac-Frenkel principle). This results in a reduced order model which takes the form of a dynamical system.

In this paper, we propose a new model order reduction method for solving general parameter-dependent dynamical systems of the form (1). This method provides an approximation of the form (4), so that it can be interpreted as a dynamical low-rank approximation method. However, the proposed method differs from existing dynamical low-rank approximation methods in that it adopts a subspace point of view and provides a uniform control of the error over the parameter set. The reduced space $X_r(t) \subset X$ is here obtained as the span of some selected trajectories $\{u(\cdot, \xi^i)\}_{i=1}^r$ of the full-order model (1). The approximation is obtained by solving a reduced dynamical system of size $r$ obtained by Galerkin projection (see, e.g., [4]), which has the form of the initial dynamical system with a modified flux which takes into account the time dependency of the subspace. The resulting approximation (when discarding some necessary numerical approximations) interpolates the solution map $\xi \mapsto u(\cdot, \xi)$ at points $\{\xi^1, \dots, \xi^r\}$. An a posteriori error estimate $\Delta_r(t, \xi)$ (local in time) is derived in the lines of [34] using the logarithmic Lipschitz constant associated to the flux. This error estimate is used in a greedy procedure for the adaptive selection of interpolation points, where at step $r$, the next interpolation point corresponds to a maximizer over the parameters set $\Xi$ of a certain norm of $t \mapsto \Delta_r(t, \xi)$.

The paper is organized as follows. In Lemma 2, we introduce the Galerkin method for the projection of the dynamical system (1) onto a time-dependent reduced space and we derive an a posteriori error estimate. In Lemma 3, we present strategies for the construction of reduced spaces $X_r$, including the classical POD-greedy strategy for the construction of
time-independent reduced spaces, and the proposed greedy algorithm for the construction of time-dependent reduced spaces. In Lemma 4, we detail some practical aspects for an implementation of the proposed method in a discrete time setting, and for obtaining online computations (solution of the reduced dynamical system and evaluation of the error estimate) with a complexity independent of the dimension $d$ of the full-order model. In Lemma 5, the proposed method is illustrated through numerical experiments on several test cases and compared to the POD-greedy approach.

# 2 Reduced dynamical system

In this section, we first propose a Galerkin method for computing a projection of the solution of the dynamical system (1) onto a time-dependent subspace $X_r(t)$ of $X$. Then, we derive an a posteriori error estimate. Here, the reduced space $X_r(t)$ is supposed to be given.

## 2.1 Projection

The state space $X = \mathbb{R}^d$ is equipped with the canonical inner product $\langle x, y \rangle_X := x^T y$ and associated norm $\|x\|_X$. We assume that

$$ \dim X_r(t) = r \quad \text{for all } t > 0. \qquad (5) $$

We denote by $\{v_1(t), \dots, v_r(t)\}$ an orthonormal basis of $X_r(t)$ and we introduce the orthogonal matrix

$$ V_r(t) = [v_1(t), \dots, v_r(t)] \in \mathbb{R}^{d \times r}, $$

with $V_r^T(t)V_r(t) = I_r$, where $I_r$ is the identity matrix of size $r \times r$. We denote by $\Pi_{X_r}(t) = V_r(t)V_r^T(t) \in \mathbb{R}^{d \times d}$ the orthogonal projector onto $X_r(t)$.

We define the approximation $u_r(t, \xi)$ of $u(t, \xi)$ by projecting the equations of the dynamical system (1) onto $X_r(t)$

$$ \left\{ \begin{array}{l} \Pi_{X_r}(t)u_r'(t, \xi) = \Pi_{X_r}(t)f(u_r(t, \xi), t, \xi), \\ u_r(0, \xi) = \Pi_{X_r}(0)u_r^0(\xi). \end{array} \right. \qquad (6) $$

Note that for subspaces $X_r(t)$ generated by trajectories of the full-order model, $X_r(0)$ contains evaluations of the initial condition $u_r^0(\xi)$, so that $\Pi_{X_r}(0)u_r^0 = u_r^0$ in the case of a parameter-independent initial condition<fnref n="1" />. We define

$$ u_r(t, \xi) = V_r(t)\alpha_r(t, \xi), $$

<footnote n="1">¹Note that the dimension of $X_r(0)$ may be different from the dimension of $X_r(t)$, $t > 0$, which does not contradict the assumption (5).</footnote>
time-independent reduced spaces, and the proposed greedy algorithm for the construction of time-dependent reduced spaces. In Lemma 4, we detail some practical aspects for an implementation of the proposed method in a discrete time setting, and for obtaining online computations (solution of the reduced dynamical system and evaluation of the error estimate) with a complexity independent of the dimension $d$ of the full-order model. In Lemma 5, the proposed method is illustrated through numerical experiments on several test cases and compared to the POD-greedy approach.

# 2 Reduced dynamical system

In this section, we first propose a Galerkin method for computing a projection of the solution of the dynamical system (1) onto a time-dependent subspace $X_r(t)$ of $X$. Then, we derive an a posteriori error estimate. Here, the reduced space $X_r(t)$ is supposed to be given.

## 2.1 Projection

The state space $X = \mathbb{R}^d$ is equipped with the canonical inner product $\langle x, y \rangle_X := x^T y$ and associated norm $\|x\|_X$. We assume that

$$ \dim X_r(t) = r \quad \text{for all } t > 0. \qquad (5) $$

We denote by $\{v_1(t), \dots, v_r(t)\}$ an orthonormal basis of $X_r(t)$ and we introduce the orthogonal matrix

$$ V_r(t) = [v_1(t), \dots, v_r(t)] \in \mathbb{R}^{d \times r}, $$

with $V_r^T(t)V_r(t) = I_r$, where $I_r$ is the identity matrix of size $r \times r$. We denote by $\Pi_{X_r}(t) = V_r(t)V_r^T(t) \in \mathbb{R}^{d \times d}$ the orthogonal projector onto $X_r(t)$.

We define the approximation $u_r(t, \xi)$ of $u(t, \xi)$ by projecting the equations of the dynamical system (1) onto $X_r(t)$

$$ \left\{ \begin{array}{l} \Pi_{X_r}(t)u_r'(t, \xi) = \Pi_{X_r}(t)f(u_r(t, \xi), t, \xi), \\ u_r(0, \xi) = \Pi_{X_r}(0)u_r^0(\xi). \end{array} \right. \qquad (6) $$

Note that for subspaces $X_r(t)$ generated by trajectories of the full-order model, $X_r(0)$ contains evaluations of the initial condition $u_r^0(\xi)$, so that $\Pi_{X_r}(0)u_r^0 = u_r^0$ in the case of a parameter-independent initial condition<fnref n="1" />. We define

$$ u_r(t, \xi) = V_r(t)\alpha_r(t, \xi), $$

<footnote n="1">¹Note that the dimension of $X_r(0)$ may be different from the dimension of $X_r(t)$, $t > 0$, which does not contradict the assumption (5).</footnote>
with $\alpha_r(t, \xi) = \mathbf{V}_r(t)^T \mathbf{u}_r(t, \xi) \in \mathbb{R}^r$. Let $\Pi_{X_r^\perp}(t) = \mathbf{I}_d - \Pi_{X_r}(t)$ be the orthogonal projector onto $X_r^\perp(t)$ which is the orthogonal complement of $X_r(t)$ in $X$. From the first equation of (6) and $\Pi_{X_r^\perp}(t)\mathbf{V}_r(t) = 0$, we deduce

$$ \mathbf{u}'_r(t, \xi) = \Pi_{X_r}(t) \mathbf{f}(\mathbf{u}_r(t, \xi), t, \xi) + \Pi_{X_r^\perp}(t) \mathbf{V}'_r(t) \alpha_r(t, \xi). \quad (7) $$

A reduced dynamical system of dimension $r$ is then obtained for the components $\alpha_r$:

$$ \left\{ \begin{array}{l} \alpha'_r(t, \xi) = \mathbf{f}_r(\alpha_r(t, \xi), t, \xi), \\ \alpha_r(0, \xi) = \alpha_r^0(\xi), \end{array} \right. \qquad (8) $$

with $\alpha_r^0(\xi) = \mathbf{V}_r^T(0)\mathbf{u}^0(\xi)$ and a reduced flux $\mathbf{f}_r$ defined by

$$ \mathbf{f}_r(\alpha_r, t, \xi) = \mathbf{V}_r(t)^T \mathbf{f}(\mathbf{V}_r(t) \alpha_r, t, \xi) - \mathbf{V}_r(t)^T \mathbf{V}'_r(t) \alpha_r, \quad (9) $$

where the last term takes into account the time dependency of the reduced basis. For a time-independent reduced basis, i.e. such that $\mathbf{V}'_r(t) = 0$, we recover a classical projected dynamical system (see, e.g., [18]).

**Remark 1.** Assuming that for a fixed $\xi$, $\mathbf{v} \mapsto \mathbf{f}(\mathbf{v}, t, \xi)$ is uniformly Lipschitz continuous, $(\mathbf{v}, t) \mapsto \mathbf{f}(\mathbf{v}, t, \xi)$ is continuous, and assuming that $\mathbf{V}_r$ is continuously differentiable with $\mathbf{V}'_r$ uniformly bounded, then $\mathbf{v} \mapsto \mathbf{f}_r(\mathbf{v}, t, \xi)$ is uniformly Lipschitz continuous with respect to $\mathbf{v}$ and $(\mathbf{v}, t) \mapsto \mathbf{f}_r(\mathbf{v}, t, \xi)$ is continuous. Then, the reduced dynamical system (8) admits a unique solution $t \mapsto \alpha_r(t, \xi)$ which is continuously differentiable.

## 2.2 A posteriori error estimate

From equations (1) and (7), we deduce that the error $e_r = \mathbf{u} - \mathbf{u}_r$ satisfies

$$ \begin{aligned} e'_r(t, \xi) = & \underbrace{\mathbf{f}(\mathbf{u}(t, \xi), t, \xi) - \mathbf{f}(\mathbf{u}_r(t, \xi), t, \xi)}_{\delta_1} + \underbrace{\Pi_{X_r^\perp}(t) \mathbf{f}(\mathbf{u}_r(t, \xi), t, \xi)}_{\delta_2} \\ & - \underbrace{\Pi_{X_r^\perp}(t) \mathbf{V}'_r(t) \mathbf{V}_r(t)^T \mathbf{u}_r(t, \xi)}_{\delta_3}, \end{aligned} \quad (10) $$

with $e_r(0, \xi) = \Pi_{X_r^\perp}(0)\mathbf{u}^0(\xi)$. The time derivative of the error $e_r$ is the sum of three contributions: the error $\delta_1$ between the flux and its approximation, the error $\delta_2$ between the flux approximation and its projection onto $X_r(t)$, and an additional term $\delta_3$ taking into account the time dependency of the basis.

We first recall the definition of the local logarithmic Lipschitz constant, as defined in [34, §2.1]. This constant will provide a local information on the flux $\mathbf{f}$ around the approximation $\mathbf{u}_r$ and will allow to derive an a posteriori error estimate.
with $\alpha_r(t, \xi) = \mathbf{V}_r(t)^T \mathbf{u}_r(t, \xi) \in \mathbb{R}^r$. Let $\Pi_{X_r^\perp}(t) = \mathbf{I}_d - \Pi_{X_r}(t)$ be the orthogonal projector onto $X_r^\perp(t)$ which is the orthogonal complement of $X_r(t)$ in $X$. From the first equation of (6) and $\Pi_{X_r^\perp}(t)\mathbf{V}_r(t) = 0$, we deduce

$$ \mathbf{u}'_r(t, \xi) = \Pi_{X_r}(t) \mathbf{f}(\mathbf{u}_r(t, \xi), t, \xi) + \Pi_{X_r^\perp}(t) \mathbf{V}'_r(t) \alpha_r(t, \xi). \quad (7) $$

A reduced dynamical system of dimension $r$ is then obtained for the components $\alpha_r$:

$$ \left\{ \begin{array}{l} \alpha'_r(t, \xi) = \mathbf{f}_r(\alpha_r(t, \xi), t, \xi), \\ \alpha_r(0, \xi) = \alpha_r^0(\xi), \end{array} \right. \qquad (8) $$

with $\alpha_r^0(\xi) = \mathbf{V}_r^T(0)\mathbf{u}^0(\xi)$ and a reduced flux $\mathbf{f}_r$ defined by

$$ \mathbf{f}_r(\alpha_r, t, \xi) = \mathbf{V}_r(t)^T \mathbf{f}(\mathbf{V}_r(t)\alpha_r, t, \xi) - \mathbf{V}_r(t)^T \mathbf{V}'_r(t) \alpha_r, \quad (9) $$

where the last term takes into account the time dependency of the reduced basis. For a time-independent reduced basis, i.e. such that $\mathbf{V}'_r(t) = 0$, we recover a classical projected dynamical system (see, e.g., [18]).

**Remark 1.** Assuming that for a fixed $\xi$, $\mathbf{v} \mapsto \mathbf{f}(\mathbf{v}, t, \xi)$ is uniformly Lipschitz continuous, $(\mathbf{v}, t) \mapsto \mathbf{f}(\mathbf{v}, t, \xi)$ is continuous, and assuming that $\mathbf{V}_r$ is continuously differentiable with $\mathbf{V}'_r$ uniformly bounded, then $\mathbf{v} \mapsto \mathbf{f}_r(\mathbf{v}, t, \xi)$ is uniformly Lipschitz continuous with respect to $\mathbf{v}$ and $(\mathbf{v}, t) \mapsto \mathbf{f}_r(\mathbf{v}, t, \xi)$ is continuous. Then, the reduced dynamical system (8) admits a unique solution $t \mapsto \alpha_r(t, \xi)$ which is continuously differentiable.

## 2.2 A posteriori error estimate

From equations (1) and (7), we deduce that the error $e_r = \mathbf{u} - \mathbf{u}_r$ satisfies

$$ \begin{aligned} e'_r(t, \xi) = & \underbrace{\mathbf{f}(\mathbf{u}(t, \xi), t, \xi) - \mathbf{f}(\mathbf{u}_r(t, \xi), t, \xi)}_{\delta_1} + \underbrace{\Pi_{X_r^\perp}(t) \mathbf{f}(\mathbf{u}_r(t, \xi), t, \xi)}_{\delta_2} \\ & - \underbrace{\Pi_{X_r^\perp}(t) \mathbf{V}'_r(t) \mathbf{V}_r(t)^T \mathbf{u}_r(t, \xi)}_{\delta_3}, \end{aligned} \quad (10) $$

with $e_r(0, \xi) = \Pi_{X_r^\perp}(0)\mathbf{u}^0(\xi)$. The time derivative of the error $e_r$ is the sum of three contributions: the error $\delta_1$ between the flux and its approximation, the error $\delta_2$ between the flux approximation and its projection onto $X_r(t)$, and an additional term $\delta_3$ taking into account the time dependency of the basis.

We first recall the definition of the local logarithmic Lipschitz constant, as defined in [34, §2.1]. This constant will provide a local information on the flux $\mathbf{f}$ around the approximation $\mathbf{u}_r$ and will allow to derive an a posteriori error estimate.
**Definition 1** (Local logarithmic Lipschitz constant). For a Lipschitz continuous function $f: X \to X$, the local logarithmic Lipschitz constant of $f$ at $v \in X$ is defined by

$$L_X[f](v) = \sup_{u \in X, u \neq v} \frac{\langle u - v, f(u) - f(v) \rangle_X}{\|u - v\|_X^2}.$$

For an affine flux, the local logarithmic Lipschitz constant is constant and can be obtained by solving an eigenvalue problem. Indeed, if $f(u) = Au + g$ with $A \in \mathbb{R}^{d \times d}$ and $g \in \mathbb{R}^d$, then

$$L_X[f](v) = L_X[A] = \sup_{0 \neq u \in X} \frac{\langle u, Au \rangle_X}{\|u\|_X^2} = \lambda_{max} \left( \frac{A + A^T}{2} \right), \quad (11)$$

where $\lambda_{max}(M)$ denotes the maximum eigenvalue of the matrix $M \in \mathbb{R}^{d \times d}$.

Now, we recall a comparison lemma [34] in a form suitable for the derivation of our a posteriori error estimate.

**Lemma 1** (Comparison lemma). Let $T > 0$ and let $u, \alpha, \beta : [0, T] \to \mathbb{R}$ be integrable functions. Assume that $u$ is differentiable and $u' \le \beta u + \alpha$. Then for all $t \in [0, T]$, it holds $u(t) \le v(t)$, where $v(t)$ is the solution of the differential equation $v' = \beta v + \alpha$ with initial condition $v(0) = u(0)$.

*Proof.* Denoting $\gamma(t) = \int_0^t \beta(\tau)d\tau$, we have

$$\begin{aligned} u(t) &= e^{\gamma(t)}e^{-\gamma(t)}u(t) - e^{\gamma(t)}u(0) + e^{\gamma(t)}u(0), \\ &= e^{\gamma(t)}\int_0^t (e^{-\gamma(\tau)}u(\tau))'d\tau + e^{\gamma(t)}v(0), \\ &\le e^{\gamma(t)}\left(\int_0^t e^{-\gamma(\tau)}\alpha(\tau)d\tau + v(0)\right) = v(t), \end{aligned}$$

which ends the proof.

<figure><img src="image_1.png" /></figure>

Following [34], we now provide a bound for the error norm $\|e_r(t, \xi)\|_X$.

**Proposition 1.** The error norm $\|e_r(t, \xi)\|_X$ satisfies

$$\|e_r(t, \xi)\|_X \le \Delta_r(t, \xi) \quad (12)$$

for all $t \ge 0$, where $\Delta_r(t, \xi)$ is the solution of the ordinary differential equation

$$\left\{ \begin{array}{l} \Delta'_{r}(t, \xi) = L_{X}[f](u_{r}(t, \xi))\Delta_{r}(t, \xi) + \|r(t, \xi)\|_{X}, \\ \Delta_{r}(0, \xi) = \|e_{r}(0, \xi)\|_{X}, \end{array} \right. \quad (13)$$

with $r(t, \xi) = \Pi_{X_r^\perp}(t)(V'_r(t)V_r(t))^T u_r(t, \xi) - f(u_r(t, \xi), t, \xi))$ and $e_r(0, \xi) = u^0(\xi) - \Pi_{X_r}(0)u^0(\xi)$.
**Definition 1** (Local logarithmic Lipschitz constant). For a Lipschitz continuous function $f: X \to X$, the local logarithmic Lipschitz constant of $f$ at $v \in X$ is defined by

$$L_X[f](v) = \sup_{u \in X, u \neq v} \frac{\langle u - v, f(u) - f(v) \rangle_X}{\|u - v\|_X^2}.$$

For an affine flux, the local logarithmic Lipschitz constant is constant and can be obtained by solving an eigenvalue problem. Indeed, if $f(u) = Au + g$ with $A \in \mathbb{R}^{d \times d}$ and $g \in \mathbb{R}^d$, then

$$L_X[f](v) = L_X[A] = \sup_{0 \neq u \in X} \frac{\langle u, Au \rangle_X}{\|u\|_X^2} = \lambda_{max} \left( \frac{A + A^T}{2} \right), \quad (11)$$

where $\lambda_{max}(M)$ denotes the maximum eigenvalue of the matrix $M \in \mathbb{R}^{d \times d}$.

Now, we recall a comparison lemma [34] in a form suitable for the derivation of our a posteriori error estimate.

**Lemma 1** (Comparison lemma). Let $T > 0$ and let $u, \alpha, \beta: [0, T] \to \mathbb{R}$ be integrable functions. Assume that $u$ is differentiable and $u' \le \beta u + \alpha$. Then for all $t \in [0, T]$, it holds $u(t) \le v(t)$, where $v(t)$ is the solution of the differential equation $v' = \beta v + \alpha$ with initial condition $v(0) = u(0)$.

*Proof.* Denoting $\gamma(t) = \int_0^t \beta(\tau)d\tau$, we have

$$\begin{aligned} u(t) &= e^{\gamma(t)}e^{-\gamma(t)}u(t) - e^{\gamma(t)}u(0) + e^{\gamma(t)}u(0), \\ &= e^{\gamma(t)}\int_0^t(e^{-\gamma(\tau)}u(\tau))'d\tau + e^{\gamma(t)}v(0), \\ &\le e^{\gamma(t)}\left(\int_0^t e^{-\gamma(\tau)}\alpha(\tau)d\tau + v(0)\right) = v(t), \end{aligned}$$

which ends the proof.

<figure><img src="image_2.png" /></figure>

Following [34], we now provide a bound for the error norm $\|e_r(t, \xi)\|_X$.

**Proposition 1.** The error norm $\|e_r(t, \xi)\|_X$ satisfies

$$\|e_r(t, \xi)\|_X \le \Delta_r(t, \xi) \quad (12)$$

for all $t \ge 0$, where $\Delta_r(t, \xi)$ is the solution of the ordinary differential equation

$$\left\{ \begin{array}{l} \Delta'_{r}(t, \xi) = L_{X}[f](u_{r}(t, \xi))\Delta_{r}(t, \xi) + \|r(t, \xi)\|_{X}, \\ \Delta_{r}(0, \xi) = \|e_{r}(0, \xi)\|_{X}, \end{array} \right. \quad (13)$$

with $r(t, \xi) = \Pi_{X_r^\perp}(t)(V'_r(t)V_r(t))^T u_r(t, \xi) - f(u_r(t, \xi), t, \xi))$ and $e_r(0, \xi) = u^0(\xi) - \Pi_{X_r}(0)u^0(\xi)$.
Proof. For the sake of clarity, we omit the dependence on $\xi$ in the proof. Taking the scalar product of (10) with $e_r(t)$, it comes

$$ 
\begin{aligned}
\frac{1}{2} \frac{d}{dt} \|e_r(t)\|_X^2 &= \langle e_r(t), e'_r(t) \rangle_X, \\
&= \langle e_r(t), (\mathbf{f}(\mathbf{u}(t), t) - \mathbf{f}(\mathbf{u}_r(t), t)) + \mathbf{r}(t) \rangle_X, \\
&\le L_X[\mathbf{f}](\mathbf{u}_r(t)) \|e_r(t)\|_X^2 + \|e_r(t)\|_X \|\mathbf{r}(t)\|_X,
\end{aligned}
$$

where we have used the definition 1 of the local logarithmic Lipschitz constant at $\mathbf{u}_r(t)$. We then obtain $\frac{d}{dt} \|e_r(t)\|_X \le L_X[\mathbf{f}](\mathbf{u}_r(t)) \|e_r(t)\|_X + \|\mathbf{r}(t)\|_X$, and the result follows from Lemma 1. □

In practice, the error bound $\Delta_r(t)$ will be estimated by solving approximately the differential equation (13) using a numerical scheme (see Lemma 4).

**Remark 2.** Following the proof of Theorem 1, we easily prove that the solution of (13), with $L_X[\mathbf{f}](\mathbf{u}_r)$ replaced by the local Lipschitz constant

$$ K_X[\mathbf{f}](\mathbf{u}_r) = \sup_{v \in X, v \neq \mathbf{u}_r} \frac{\|\mathbf{f}(v) - \mathbf{f}(\mathbf{u}_r)\|_X}{\|v - \mathbf{u}_r\|_X}, $$

also provides an error bound. However, since $L_X[\mathbf{f}] \le K_X[\mathbf{f}]$, a sharper error bound is obtained when using the local logarithmic Lipschitz constant. Moreover, $L_X[\mathbf{f}](\mathbf{u}_r)$ is easier to compute than $K_X[\mathbf{f}](\mathbf{u}_r)$ in practice.

**Remark 3.** The error bound $\Delta_r(t)$ may grow exponentially with time. However, it is a certified error bound and a good candidate for piloting adaptive algorithms (see later). Moreover, if particular information is available on the structure of the problem (e.g. if it arises from spatial discretization of a PDE), it is possible to improve this error bound using energy norms [20, §3.3]. In the context of Petrov-Galerkin space-time formulations, better error estimates can be obtained for particular classes of problems, e.g. [32] for parabolic equation and [35] for Burgers' equation. For the general non linear dynamical systems considered here, a way to improve the effectivity of error bound $\Delta_r$ is to improve the local Lipschitz constant, as suggested in [34, §5].

# 3 Greedy algorithms for the construction of reduced spaces

In this section, we introduce greedy algorithms for the construction of an increasing sequence of reduced spaces $\{X_r\}_{r>0}$. These spaces are generated from successive evaluations of the solution $\mathbf{u}$ of the full-order model at parameters values $\{\xi^r\}_{r>0}$ which are selected adaptively.
Proof. For the sake of clarity, we omit the dependence on $\xi$ in the proof. Taking the scalar product of (10) with $e_r(t)$, it comes

$$ 
\begin{aligned}
\frac{1}{2} \frac{d}{dt} \|e_r(t)\|_X^2 &= \langle e_r(t), e'_r(t) \rangle_X, \\
&= \langle e_r(t), (\mathbf{f}(\mathbf{u}(t), t) - \mathbf{f}(\mathbf{u}_r(t), t)) + \mathbf{r}(t) \rangle_X, \\
&\le L_X[\mathbf{f}](\mathbf{u}_r(t)) \|e_r(t)\|_X^2 + \|e_r(t)\|_X \|\mathbf{r}(t)\|_X,
\end{aligned}
$$

where we have used the definition 1 of the local logarithmic Lipschitz constant at $\mathbf{u}_r(t)$. We then obtain $\frac{d}{dt} \|e_r(t)\|_X \le L_X[\mathbf{f}](\mathbf{u}_r(t)) \|e_r(t)\|_X + \|\mathbf{r}(t)\|_X$, and the result follows from Lemma 1. □

In practice, the error bound $\Delta_r(t)$ will be estimated by solving approximately the differential equation (13) using a numerical scheme (see Lemma 4).

**Remark 2.** Following the proof of Theorem 1, we easily prove that the solution of (13), with $L_X[\mathbf{f}](\mathbf{u}_r)$ replaced by the local Lipschitz constant

$$ K_X[\mathbf{f}](\mathbf{u}_r) = \sup_{v \in X, v \neq \mathbf{u}_r} \frac{\|\mathbf{f}(v) - \mathbf{f}(\mathbf{u}_r)\|_X}{\|v - \mathbf{u}_r\|_X}, $$

also provides an error bound. However, since $L_X[\mathbf{f}] \le K_X[\mathbf{f}]$, a sharper error bound is obtained when using the local logarithmic Lipschitz constant. Moreover, $L_X[\mathbf{f}](\mathbf{u}_r)$ is easier to compute than $K_X[\mathbf{f}](\mathbf{u}_r)$ in practice.

**Remark 3.** The error bound $\Delta_r(t)$ may grow exponentially with time. However, it is a certified error bound and a good candidate for piloting adaptive algorithms (see later). Moreover, if particular information is available on the structure of the problem (e.g. if it arises from spatial discretization of a PDE), it is possible to improve this error bound using energy norms [20, §3.3]. In the context of Petrov-Galerkin space-time formulations, better error estimates can be obtained for particular classes of problems, e.g. [32] for parabolic equation and [35] for Burgers' equation. For the general non linear dynamical systems considered here, a way to improve the effectivity of error bound $\Delta_r$ is to improve the local Lipschitz constant, as suggested in [34, §5].

# 3 Greedy algorithms for the construction of reduced spaces

In this section, we introduce greedy algorithms for the construction of an increasing sequence of reduced spaces $\{X_r\}_{r>0}$. These spaces are generated from successive evaluations of the solution $\mathbf{u}$ of the full-order model at parameters values $\{\xi^r\}_{r>0}$ which are selected adaptively.
For a given subspace $X_r$ in $X$, possibly time-dependent, we consider $u_r(\cdot, \xi)$ the solution of (6). We assume that an a posteriori estimate $\Delta_r(t, \xi)$ of the error $\|e_r(t, \xi)\|_X$ is available (see Subsection 2.2). Defining

$$X_{r+1}(t) = X_r(t) + \operatorname{span}\{u(t, \xi^{r+1}(t))\},$$

where $\xi^{r+1}(t)$ maximizes $\xi \mapsto \Delta_r(t, \xi)$ over $\Xi$, would result in a standard greedy algorithm for the approximation of the solution manifold $\{u(t, \xi) : \xi \in \Xi\}$ at time $t$ with a sequence of subspaces $\{X_r(t)\}_{r>0}$ (see [6, 5, 11] for convergence results on greedy algorithms). However, choosing time-dependent parameters values is infeasible in practice, even when working with a discretization on a time grid. Indeed, it requires the solution of the full-order model for too many parameters values.

Therefore, we rely on evaluations of the solution at time-independent parameters values whose selection requires a global in time error estimate. Here, we introduce an a posteriori estimate of the $L^2(0, T)$-norm of the error $\|e_r(\cdot, \xi)\|_X$ defined by

$$\Delta_r^{(0,T)}(\xi) := \|\Delta_r(\cdot, \xi)\|_{(0,T),2} = \left(\int_0^T \Delta_r(t, \xi)^2 dt\right)^{1/2},$$

where $\| \cdot \|_{(0,T),2}$ denotes the natural norm in $L^2(0, T)$.

**Remark 4.** Note that other norms of $t \mapsto \Delta_r(t, \xi)$ (e.g., the $L^\infty(0, T)$-norm or a weighted $L^2(0, T)$-norm) could be used for defining the error indicator $\Delta_r^{(0,T)}(\xi)$, depending on how we want to control the quality of the reduced-order model.

## 3.1 T-greedy algorithm

We first propose a very natural strategy which consists in defining

$$X_{r+1}(t) = X_r(t) + \operatorname{span}\{u(t, \xi^{r+1})\}.$$

A basic strategy would consist in choosing the sequence of parameters values $\{\xi^1, \dots, \xi^r, \dots\}$ at random. In this work, we adopt an adaptive greedy strategy called T-greedy algorithm, where $\xi^{r+1}$ maximizes over $\Xi$ the error estimate $\Delta_r^{(0,T)}(\xi)$. Note that in practice, parameters values are selected in a finite subset $\Xi_{train}$ in $\Xi$ called a training set. The T-Greedy algorithm is summarized in Algorithm 1.

In practice, we can fix a desired precision $\varepsilon$ and stop the Algorithm 1 after step 2 if $\max_{\xi \in \Xi_{train}} \Delta_r^{(0,T)}(\xi) < \varepsilon$. At iteration $r$ of the algorithm, the approximation $u_r(t)$ takes the form

$$u_r(t, \xi) = \sum_{i=1}^{r} v_i(t) \alpha_i(t, \xi), \qquad (14)$$
For a given subspace $X_r$ in $X$, possibly time-dependent, we consider $u_r(\cdot, \xi)$ the solution of (6). We assume that an a posteriori estimate $\Delta_r(t, \xi)$ of the error $\|e_r(t, \xi)\|_X$ is available (see Subsection 2.2). Defining

$$X_{r+1}(t) = X_r(t) + \operatorname{span}\{u(t, \xi^{r+1}(t))\},$$

where $\xi^{r+1}(t)$ maximizes $\xi \mapsto \Delta_r(t, \xi)$ over $\Xi$, would result in a standard greedy algorithm for the approximation of the solution manifold $\{u(t, \xi) : \xi \in \Xi\}$ at time $t$ with a sequence of subspaces $\{X_r(t)\}_{r>0}$ (see [6, 5, 11] for convergence results on greedy algorithms). However, choosing time-dependent parameters values is infeasible in practice, even when working with a discretization on a time grid. Indeed, it requires the solution of the full-order model for too many parameters values.

Therefore, we rely on evaluations of the solution at time-independent parameters values whose selection requires a global in time error estimate. Here, we introduce an a posteriori estimate of the $L^2(0, T)$-norm of the error $\|e_r(\cdot, \xi)\|_X$ defined by

$$\Delta_r^{(0,T)}(\xi) := \|\Delta_r(\cdot, \xi)\|_{(0,T),2} = \left(\int_0^T \Delta_r(t, \xi)^2 dt\right)^{1/2},$$

where $\| \cdot \|_{(0,T),2}$ denotes the natural norm in $L^2(0, T)$.

**Remark 4.** Note that other norms of $t \mapsto \Delta_r(t, \xi)$ (e.g., the $L^\infty(0, T)$-norm or a weighted $L^2(0, T)$-norm) could be used for defining the error indicator $\Delta_r^{(0,T)}(\xi)$, depending on how we want to control the quality of the reduced-order model.

## 3.1 T-greedy algorithm

We first propose a very natural strategy which consists in defining

$$X_{r+1}(t) = X_r(t) + \operatorname{span}\{u(t, \xi^{r+1})\}.$$

A basic strategy would consist in choosing the sequence of parameters values $\{\xi^1, \dots, \xi^r, \dots\}$ at random. In this work, we adopt an adaptive greedy strategy called T-greedy algorithm, where $\xi^{r+1}$ maximizes over $\Xi$ the error estimate $\Delta_r^{(0,T)}(\xi)$. Note that in practice, parameters values are selected in a finite subset $\Xi_{train}$ in $\Xi$ called a training set. The T-Greedy algorithm is summarized in Algorithm 1.

In practice, we can fix a desired precision $\varepsilon$ and stop the Algorithm 1 after step 2 if $\max_{\xi \in \Xi_{train}} \Delta_r^{(0,T)}(\xi) < \varepsilon$. At iteration $r$ of the algorithm, the approximation $u_r(t)$ takes the form

$$u_r(t, \xi) = \sum_{i=1}^{r} v_i(t) \alpha_i(t, \xi), \qquad (14)$$
**Algorithm 1 T-Greedy algorithm**

1: Set $r = 0, X_0 = \{0\}$.

2: Compute $u_r(\cdot, \xi)$ and $\Delta_r(\cdot, \xi)$ for all $\xi$ in $\Xi_{train}$.

3: Find $\xi^{r+1} \in \arg\max_{\xi \in \Xi_{train}} \Delta_r^{(0,T)}(\xi)$.

4: Compute $t \mapsto u(t, \xi^{r+1})$.

5: Set $X_{r+1}(t) = X_r(t) + \text{span}\{u(t, \xi^{r+1})\}$.

6: Set $r \leftarrow r + 1$ and go to step 2.

where $\{v_1(t), \dots, v_r(t)\}$ constitutes a basis of the space $X_r(t)$. The time-dependent sub-space $X_r(t)$ contains the solution $u(t, \xi)$ of the full-order model for parameters values $\xi^1, \dots, \xi^r$. Therefore, by the property of the Galerkin projection, the solution $u_r$ of the reduced dynamical system interpolates the solution $u$ at these parameters values.

Satisfying assumption (5) requires that for all $t > 0$, the vectors $u(t, \xi^1), \dots, u(t, \xi^r)$ are linearly independent. Under this assumption, we can define the basis of $X_r(t)$ by

$$v_i(t) = u(t, \xi^i) \quad \text{or} \quad v_i(t) = \frac{(I_d - \Pi_{X_{i-1}}(t))u(t, \xi^i)}{\|(I_d - \Pi_{X_{i-1}}(t))u(t, \xi^i)\|_X}, \quad (15)$$

the second choice resulting in an orthonormal basis which is more convenient for numerical stability issues. Moreover, assuming that the flux $f(\cdot, \cdot, \xi^i)$ is continuous and uniformly Lipschitz continuous with respect to its first variable, the trajectories $t \mapsto u(t, \xi^i)$, and therefore the functions $t \mapsto v_i(t)$, are continuously differentiable.

## 3.2 POD-Greedy algorithm

Here, we describe the POD-greedy algorithm introduced in [17] for the adaptive construction of time-independent reduced spaces $X_r$. Given an approximation space $X_r$ and the corresponding solution $u_r$ of the reduced dynamical system in $X_r$, we select a new parameters value $\xi^{r+1}$ which maximizes the error estimate $\xi \mapsto \Delta_r^{(0,T)}(\xi)$, such as for the T-greedy strategy. Then the trajectory $u(\cdot, \xi^{r+1})$ of the full-order model is computed and the space $X_r$ is enriched by the $\ell$-dimensional subspace generated by the first $\ell$ POD modes $\text{POD}_\ell(s_r(\cdot, \xi^{r+1}))$ of $s_r(\cdot, \xi^{r+1}) = u(\cdot, \xi^{r+1}) - \Pi_{X_r} u(\cdot, \xi^{r+1})$, such that (see, e.g., [33])

$$\text{span}(\text{POD}_\ell(s_r(\cdot, \xi^{r+1}))) = \arg \min_{\text{dim}(V)=\ell} \int_0^T \|s_r(t, \xi^{r+1}) - \Pi_V s_r(t, \xi^{r+1})\|_X^2 dt.$$

The algorithm is summarized in Algorithm 2.

As for Algorithm 1, we can fix a desired precision $\varepsilon$ and stop the Algorithm 2 after step 2 if $\arg\max_{\xi \in \Xi_{train}} \Delta_r^{(0,T)}(\xi) < \varepsilon$. Let us remark that choosing $\ell = 1$ allows to obtain a slow increase of the dimension of $X_r$ along the iterations. Since the reduced space is
**Algorithm 1** T-Greedy algorithm

1: Set $r = 0, X_0 = \{0\}$.

2: Compute $\mathbf{u}_r(\cdot, \boldsymbol{\xi})$ and $\Delta_r(\cdot, \boldsymbol{\xi})$ for all $\boldsymbol{\xi}$ in $\Xi_{train}$.

3: Find $\boldsymbol{\xi}^{r+1} \in \arg\max_{\boldsymbol{\xi} \in \Xi_{train}} \Delta_r^{(0,T)}(\boldsymbol{\xi})$.

4: Compute $t \mapsto \mathbf{u}(t, \boldsymbol{\xi}^{r+1})$.

5: Set $X_{r+1}(t) = X_r(t) + \text{span}\{\mathbf{u}(t, \boldsymbol{\xi}^{r+1})\}$.

6: Set $r \leftarrow r + 1$ and go to step 2.

where $\{v_1(t), \dots, v_r(t)\}$ constitutes a basis of the space $X_r(t)$. The time-dependent sub-space $X_r(t)$ contains the solution $\mathbf{u}(t, \boldsymbol{\xi})$ of the full-order model for parameters values $\boldsymbol{\xi}^1, \dots, \boldsymbol{\xi}^r$. Therefore, by the property of the Galerkin projection, the solution $\mathbf{u}_r$ of the reduced dynamical system interpolates the solution $\mathbf{u}$ at these parameters values.

Satisfying assumption (5) requires that for all $t > 0$, the vectors $\mathbf{u}(t, \boldsymbol{\xi}^1), \dots, \mathbf{u}(t, \boldsymbol{\xi}^r)$ are linearly independent. Under this assumption, we can define the basis of $X_r(t)$ by

$$ v_i(t) = \mathbf{u}(t, \boldsymbol{\xi}^i) \quad \text{or} \quad v_i(t) = \frac{(\mathbf{I}_d - \mathbf{\Pi}_{X_{i-1}}(t))\mathbf{u}(t, \boldsymbol{\xi}^i)}{\|(\mathbf{I}_d - \mathbf{\Pi}_{X_{i-1}}(t))\mathbf{u}(t, \boldsymbol{\xi}^i)\|_X}, \quad (15) $$

the second choice resulting in an orthonormal basis which is more convenient for numerical stability issues. Moreover, assuming that the flux $f(\cdot, \cdot, \boldsymbol{\xi}^i)$ is continuous and uniformly Lipschitz continuous with respect to its first variable, the trajectories $t \mapsto \mathbf{u}(t, \boldsymbol{\xi}^i)$, and therefore the functions $t \mapsto v_i(t)$, are continuously differentiable.

## 3.2 POD-Greedy algorithm

Here, we describe the POD-greedy algorithm introduced in [17] for the adaptive construction of time-independent reduced spaces $X_r$. Given an approximation space $X_r$ and the corresponding solution $\mathbf{u}_r$ of the reduced dynamical system in $X_r$, we select a new parameters value $\boldsymbol{\xi}^{r+1}$ which maximizes the error estimate $\boldsymbol{\xi} \mapsto \Delta_r^{(0,T)}(\boldsymbol{\xi})$, such as for the T-greedy strategy. Then the trajectory $\mathbf{u}(\cdot, \boldsymbol{\xi}^{r+1})$ of the full-order model is computed and the space $X_r$ is enriched by the $\ell$-dimensional subspace generated by the first $\ell$ POD modes $\text{POD}_\ell(\mathbf{s}_r(\cdot, \boldsymbol{\xi}^{r+1}))$ of $\mathbf{s}_r(\cdot, \boldsymbol{\xi}^{r+1}) = \mathbf{u}(\cdot, \boldsymbol{\xi}^{r+1}) - \mathbf{\Pi}_{X_r} \mathbf{u}(\cdot, \boldsymbol{\xi}^{r+1})$, such that (see, e.g., [33])

$$ \text{span}(\text{POD}_\ell(\mathbf{s}_r(\cdot, \boldsymbol{\xi}^{r+1}))) = \arg \min_{\text{dim}(V)=\ell} \int_0^T \|\mathbf{s}_r(t, \boldsymbol{\xi}^{r+1}) - \mathbf{\Pi}_V \mathbf{s}_r(t, \boldsymbol{\xi}^{r+1})\|_X^2 dt. $$

The algorithm is summarized in Algorithm 2.

As for Algorithm 1, we can fix a desired precision $\varepsilon$ and stop the Algorithm 2 after step 2 if $\arg\max_{\boldsymbol{\xi} \in \Xi_{train}} \Delta_r^{(0,T)}(\boldsymbol{\xi}) < \varepsilon$. Let us remark that choosing $\ell = 1$ allows to obtain a slow increase of the dimension of $X_r$ along the iterations. Since the reduced space is
<caption>**Algorithm 2** POD-greedy algorithm</caption>

1: Set $r = 0, X_0 = \{0\}$.
2: Compute $\mathbf{u}_r(\cdot, \boldsymbol{\xi})$ and $\Delta_r(\cdot, \boldsymbol{\xi})$ for all $\boldsymbol{\xi}$ in $\Xi_{train}$.
3: Find $\boldsymbol{\xi}^{r+1} \in \arg\max_{\boldsymbol{\xi} \in \Xi_{train}} \Delta_r^{(0,T)}(\boldsymbol{\xi})$.
4: Compute $\mathbf{u}(\cdot, \boldsymbol{\xi}^{r+1})$ and $\mathbf{s}_r(\cdot, \boldsymbol{\xi}^{r+1}) = \mathbf{u}(\cdot, \boldsymbol{\xi}^{r+1}) - \Pi_{X_r}\mathbf{u}(\cdot, \boldsymbol{\xi}^{r+1})$.
5: Set $X_{r+\ell} = X_r + \operatorname{span}(\mathrm{POD}_{\ell}(\mathbf{s}_r(\cdot, \boldsymbol{\xi}^{r+1})))$.
6: Set $r \leftarrow r + \ell$ and go to step 2.

enriched by only the first modes of the POD of the trajectories $t \mapsto \mathbf{u}(t, \boldsymbol{\xi}^i)$, $1 \le i \le r$, the approximation $\mathbf{u}_r$ does not in general interpolate the function $\boldsymbol{\xi} \mapsto \mathbf{u}(\cdot, \boldsymbol{\xi})$ at points $\boldsymbol{\xi}^1, \dots, \boldsymbol{\xi}^r$. As it will be illustrated in the numerical experiments, this enrichment strategy may yield very high-dimensional reduced spaces for reaching a desired accuracy. A typical example is the advection problem $\partial_t u + \boldsymbol{\xi} \partial_x u = 0$ on the torus with initial condition $u^0$, for which a very high-dimensional time-independent reduced space may be required to approximate the solution $t \mapsto u^0(x - \boldsymbol{\xi}t)$, even for one instance of the parameters $\boldsymbol{\xi}$. More generally, it is well known that POD is not well suited for solving problems with propagating fronts.

# 4 Practical implementation

In this section, we provide the practical ingredients for an efficient offline/online implementation of the proposed model reduction method in a discrete time setting. We first introduce a time integration scheme for the solution of the full and reduced order dynamical systems, and for the computation of the error estimate. Then, for the model order reduction method to be efficient, the computation of the approximation $\mathbf{u}_r$ as well as the evaluation of the error estimate $\Delta_r$ have to be performed online with a complexity independent of the dimension $d$ of the full-order model. This requires some assumptions on the dependence of the flux on the parameters and some pre-computations in an offline phase.

## 4.1 Time integration

Let $\mathbb{T} = \{t_k\}_{k=0}^K$ be a regular discretization of $[0, T]$ with $t_k = k\delta t$ and $\delta t = \frac{T}{K}$. Given $a: [0, T] \to V$, with $V$ a vector space, we denote by $a^k \approx a(t^k)$ an approximation of $a$ at time $t^k$ and $\delta a^k = a^{k+1} - a^k$. We assume that the flux $f: X \times (0, T) \times \Xi \to X$ can be decomposed as follows:

$$f(\mathbf{u}, t, \boldsymbol{\xi}) = \mathbf{A}(t, \boldsymbol{\xi})\mathbf{u} + \mathbf{h}(\mathbf{u}, t, \boldsymbol{\xi}) + \mathbf{g}(t, \boldsymbol{\xi}),$$

where $\mathbf{A}(t, \boldsymbol{\xi}) \in \mathbb{R}^{d \times d}$, $\mathbf{h}(\cdot, t, \boldsymbol{\xi}) : X \to \mathbb{R}^d$ and $\mathbf{g}(t, \boldsymbol{\xi}) \in \mathbb{R}^d$. For solving (1), we use the following semi-implicit time integration scheme

$$(\mathbf{I}_d - \delta t \mathbf{A}^{k+1}(\boldsymbol{\xi}))\mathbf{u}^{k+1}(\boldsymbol{\xi}) = \mathbf{u}^k(\boldsymbol{\xi}) + \delta t \mathbf{h}^k(\mathbf{u}^k(\boldsymbol{\xi}), \boldsymbol{\xi}) + \delta t \mathbf{g}^k(\boldsymbol{\xi}), \quad (16)$$
<caption>**Algorithm 2** POD-greedy algorithm</caption>

1: Set $r = 0, X_0 = \{0\}$.
2: Compute $\mathbf{u}_r(\cdot, \boldsymbol{\xi})$ and $\Delta_r(\cdot, \boldsymbol{\xi})$ for all $\boldsymbol{\xi}$ in $\Xi_{train}$.
3: Find $\boldsymbol{\xi}^{r+1} \in \arg\max_{\boldsymbol{\xi} \in \Xi_{train}} \Delta_r^{(0,T)}(\boldsymbol{\xi})$.
4: Compute $\mathbf{u}(\cdot, \boldsymbol{\xi}^{r+1})$ and $\mathbf{s}_r(\cdot, \boldsymbol{\xi}^{r+1}) = \mathbf{u}(\cdot, \boldsymbol{\xi}^{r+1}) - \Pi_{X_r}\mathbf{u}(\cdot, \boldsymbol{\xi}^{r+1})$.
5: Set $X_{r+\ell} = X_r + \text{span}(\text{POD}_{\ell}(\mathbf{s}_r(\cdot, \boldsymbol{\xi}^{r+1})))$.
6: Set $r \leftarrow r + \ell$ and go to step 2.

enriched by only the first modes of the POD of the trajectories $t \mapsto \mathbf{u}(t, \boldsymbol{\xi}^i)$, $1 \le i \le r$, the approximation $\mathbf{u}_r$ does not in general interpolate the function $\boldsymbol{\xi} \mapsto \mathbf{u}(\cdot, \boldsymbol{\xi})$ at points $\boldsymbol{\xi}^1, \dots, \boldsymbol{\xi}^r$. As it will be illustrated in the numerical experiments, this enrichment strategy may yield very high-dimensional reduced spaces for reaching a desired accuracy. A typical example is the advection problem $\partial_t u + \boldsymbol{\xi} \partial_x u = 0$ on the torus with initial condition $u^0$, for which a very high-dimensional time-independent reduced space may be required to approximate the solution $t \mapsto u^0(x - \boldsymbol{\xi}t)$, even for one instance of the parameters $\boldsymbol{\xi}$. More generally, it is well known that POD is not well suited for solving problems with propagating fronts.

# 4 Practical implementation

In this section, we provide the practical ingredients for an efficient offline/online implementation of the proposed model reduction method in a discrete time setting. We first introduce a time integration scheme for the solution of the full and reduced order dynamical systems, and for the computation of the error estimate. Then, for the model order reduction method to be efficient, the computation of the approximation $\mathbf{u}_r$ as well as the evaluation of the error estimate $\Delta_r$ have to be performed online with a complexity independent of the dimension $d$ of the full-order model. This requires some assumptions on the dependence of the flux on the parameters and some pre-computations in an offline phase.

## 4.1 Time integration

Let $\mathbb{T} = \{t_k\}_{k=0}^K$ be a regular discretization of $[0, T]$ with $t_k = k\delta t$ and $\delta t = \frac{T}{K}$. Given $a: [0, T] \to V$, with $V$ a vector space, we denote by $a^k \approx a(t^k)$ an approximation of $a$ at time $t^k$ and $\delta a^k = a^{k+1} - a^k$. We assume that the flux $f: X \times (0, T) \times \Xi \to X$ can be decomposed as follows:

$$f(\mathbf{u}, t, \boldsymbol{\xi}) = \mathbf{A}(t, \boldsymbol{\xi})\mathbf{u} + \mathbf{h}(\mathbf{u}, t, \boldsymbol{\xi}) + \mathbf{g}(t, \boldsymbol{\xi}),$$

where $\mathbf{A}(t, \boldsymbol{\xi}) \in \mathbb{R}^{d \times d}$, $\mathbf{h}(\cdot, t, \boldsymbol{\xi}) : X \to \mathbb{R}^d$ and $\mathbf{g}(t, \boldsymbol{\xi}) \in \mathbb{R}^d$. For solving (1), we use the following semi-implicit time integration scheme

$$(\mathbf{I}_d - \delta t \mathbf{A}^{k+1}(\boldsymbol{\xi}))\mathbf{u}^{k+1}(\boldsymbol{\xi}) = \mathbf{u}^k(\boldsymbol{\xi}) + \delta t \mathbf{h}^k(\mathbf{u}^k(\boldsymbol{\xi}), \boldsymbol{\xi}) + \delta t \mathbf{g}^k(\boldsymbol{\xi}), \quad (16)$$
The matrix $(\mathbf{I}_d - \delta t \mathbf{A}^{k+1}(\boldsymbol{\xi}))$ is assumed to be invertible.

### 4.1.1 Time integration of the reduced dynamical system

Let $X_r^k$ denote the reduced space at time $t^k$, $\{\boldsymbol{v}_1^k, \dots, \boldsymbol{v}_r^k\}$ an orthonormal basis of $X_r^k$, $\mathbf{V}_r^k = [\boldsymbol{v}_1^k, \dots, \boldsymbol{v}_r^k]$, $\Pi_{X_r^k} = \mathbf{V}_r^{kT}\mathbf{V}_r^k$ the orthogonal projector onto $X_r^k$, and $\Pi_{X_r^{k\perp}} = \mathbf{I}_d - \Pi_{X_r^k}$. The approximation $\boldsymbol{u}_r^k$ at time $t^k$ is obtained by projecting the discrete dynamical system (16) onto $X_r^{k+1}$

$$ \boldsymbol{u}_r^{k+1}(\boldsymbol{\xi}) = \Pi_{X_r^{k+1}} (\boldsymbol{u}_r^k(\boldsymbol{\xi}) + \delta t (\mathbf{A}^{k+1}(\boldsymbol{\xi})\boldsymbol{u}_r^{k+1}(\boldsymbol{\xi}) + \boldsymbol{h}^k(\boldsymbol{u}_r^k(\boldsymbol{\xi}), \boldsymbol{\xi}) + \boldsymbol{g}^k(\boldsymbol{\xi}))) , \quad (17) $$

with $\boldsymbol{u}_r^0(\boldsymbol{\xi}) = \Pi_{X_r^0} \boldsymbol{u}^0(\boldsymbol{\xi})$. Then, we obtain the following discrete reduced dynamical system:

$$ (\mathbf{I}_r - \delta t \mathbf{A}_r^{k+1}(\boldsymbol{\xi})) \boldsymbol{\alpha}_r^{k+1}(\boldsymbol{\xi}) = \boldsymbol{\alpha}_r^k(\boldsymbol{\xi}) + \delta t \mathbf{V}_r^{k+1T} \left( \boldsymbol{h}^k(\mathbf{V}_r^k \boldsymbol{\alpha}_r^k(\boldsymbol{\xi}), \boldsymbol{\xi}) + \boldsymbol{g}^k(\boldsymbol{\xi}) - \frac{\delta \mathbf{V}_r^k}{\delta t} \boldsymbol{\alpha}_r^k(\boldsymbol{\xi}) \right), \quad (18) $$

with initial condition $\boldsymbol{\alpha}_r^0(\boldsymbol{\xi}) = \mathbf{V}_r^{0T} \boldsymbol{u}^0(\boldsymbol{\xi})$, where

$$ \mathbf{A}_r^{k+1}(\boldsymbol{\xi}) = \mathbf{V}_r^{k+1T} \mathbf{A}^{k+1}(\boldsymbol{\xi}) \mathbf{V}_r^{k+1}, \quad (19) $$

and where $\frac{\delta \mathbf{V}_r^k}{\delta t}$ takes into account the possible time dependency of the reduced basis (this term is equal to zero for time-independent reduced spaces). In practice, for an efficient solution of (18), $\boldsymbol{h}^k(\boldsymbol{u}_r^k(\boldsymbol{\xi}), \boldsymbol{\xi})$ is replaced by an approximation $\tilde{\boldsymbol{h}}^k(\boldsymbol{u}_r^k(\boldsymbol{\xi}), \boldsymbol{\xi})$ using an empirical interpolation method (see Subsection 4.2.1).

### 4.1.2 Time integration for error estimation

We first note that if in (13) the constant $L_X[\boldsymbol{f}](\boldsymbol{u}_r)$ is replaced by an upper bound, then the solution of the ordinary differential equation provides an upper bound for $\Delta_r$. Noting that $L_X[\boldsymbol{f}](\boldsymbol{u}_r) \le L_X[\mathbf{A}] + L_X[\boldsymbol{h}](\boldsymbol{u}_r)$, we introduce estimations $\tilde{L}_X[\mathbf{A}]$ and $\tilde{L}_X[\boldsymbol{h}](\boldsymbol{u}_r)$ of $L_X[\mathbf{A}]$ and $L_X[\boldsymbol{h}](\boldsymbol{u}_r)$ respectively and consider the following ordinary differential equation whose solution $\tilde{\Delta}_r$ provides an estimation of an upper bound of $\Delta_r$:

$$ \left\{ \begin{array}{l} \tilde{\Delta}'_r(t, \boldsymbol{\xi}) = \tilde{L}_X[\mathbf{A}]\tilde{\Delta}_r(t, \boldsymbol{\xi}) + \tilde{L}_X[\boldsymbol{h}](\boldsymbol{u}_r)\tilde{\Delta}_r(t, \boldsymbol{\xi}) + \|\tilde{\boldsymbol{r}}(t, \boldsymbol{\xi})\|_X, \\ \tilde{\Delta}_r(0, \boldsymbol{\xi}) = \|\boldsymbol{e}_r(0, \boldsymbol{\xi})\|_X. \end{array} \right. \quad (20) $$

Here, $\|\tilde{\boldsymbol{r}}(t, \boldsymbol{\xi})\|_X$ denotes an approximation of $\|\boldsymbol{r}(t, \boldsymbol{\xi})\|_X$ where $\boldsymbol{h}(\boldsymbol{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi})$ is replaced by an approximation $\tilde{\boldsymbol{h}}(\boldsymbol{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi})$ obtained by an empirical interpolation method (see Subsection 4.2.1).
The matrix $(\mathbf{I}_d - \delta t \mathbf{A}^{k+1}(\boldsymbol{\xi}))$ is assumed to be invertible.

### 4.1.1 Time integration of the reduced dynamical system

Let $X_r^k$ denote the reduced space at time $t^k$, $\{\boldsymbol{v}_1^k, \dots, \boldsymbol{v}_r^k\}$ an orthonormal basis of $X_r^k$, $\mathbf{V}_r^k = [\boldsymbol{v}_1^k, \dots, \boldsymbol{v}_r^k]$, $\Pi_{X_r^k} = \mathbf{V}_r^{kT}\mathbf{V}_r^k$ the orthogonal projector onto $X_r^k$, and $\Pi_{X_r^{k\perp}} = \mathbf{I}_d - \Pi_{X_r^k}$. The approximation $\boldsymbol{u}_r^k$ at time $t^k$ is obtained by projecting the discrete dynamical system (16) onto $X_r^{k+1}$

$$ \boldsymbol{u}_r^{k+1}(\boldsymbol{\xi}) = \Pi_{X_r^{k+1}} (\boldsymbol{u}_r^k(\boldsymbol{\xi}) + \delta t (\mathbf{A}^{k+1}(\boldsymbol{\xi})\boldsymbol{u}_r^{k+1}(\boldsymbol{\xi}) + \boldsymbol{h}^k(\boldsymbol{u}_r^k(\boldsymbol{\xi}), \boldsymbol{\xi}) + \boldsymbol{g}^k(\boldsymbol{\xi}))) , \quad (17) $$

with $\boldsymbol{u}_r^0(\boldsymbol{\xi}) = \Pi_{X_r^0} \boldsymbol{u}^0(\boldsymbol{\xi})$. Then, we obtain the following discrete reduced dynamical system:

$$ (\mathbf{I}_r - \delta t \mathbf{A}_r^{k+1}(\boldsymbol{\xi})) \boldsymbol{\alpha}_r^{k+1}(\boldsymbol{\xi}) = \boldsymbol{\alpha}_r^k(\boldsymbol{\xi}) + \quad (18) $$

$$ \delta t \mathbf{V}_r^{k+1T} \left( \boldsymbol{h}^k(\mathbf{V}_r^k \boldsymbol{\alpha}_r^k(\boldsymbol{\xi}), \boldsymbol{\xi}) + \boldsymbol{g}^k(\boldsymbol{\xi}) - \frac{\delta \mathbf{V}_r^k}{\delta t} \boldsymbol{\alpha}_r^k(\boldsymbol{\xi}) \right), $$

with initial condition $\boldsymbol{\alpha}_r^0(\boldsymbol{\xi}) = \mathbf{V}_r^{0T} \boldsymbol{u}^0(\boldsymbol{\xi})$, where

$$ \mathbf{A}_r^{k+1}(\boldsymbol{\xi}) = \mathbf{V}_r^{k+1T} \mathbf{A}^{k+1}(\boldsymbol{\xi}) \mathbf{V}_r^{k+1}, \quad (19) $$

and where $\frac{\delta \mathbf{V}_r^k}{\delta t}$ takes into account the possible time dependency of the reduced basis (this term is equal to zero for time-independent reduced spaces). In practice, for an efficient solution of (18), $\boldsymbol{h}^k(\boldsymbol{u}_r^k(\boldsymbol{\xi}), \boldsymbol{\xi})$ is replaced by an approximation $\tilde{\boldsymbol{h}}^k(\boldsymbol{u}_r^k(\boldsymbol{\xi}), \boldsymbol{\xi})$ using an empirical interpolation method (see Subsection 4.2.1).

### 4.1.2 Time integration for error estimation

We first note that if in (13) the constant $L_X[\boldsymbol{f}](\boldsymbol{u}_r)$ is replaced by an upper bound, then the solution of the ordinary differential equation provides an upper bound for $\Delta_r$. Noting that $L_X[\boldsymbol{f}](\boldsymbol{u}_r) \le L_X[\mathbf{A}] + L_X[\boldsymbol{h}](\boldsymbol{u}_r)$, we introduce estimations $\tilde{L}_X[\mathbf{A}]$ and $\tilde{L}_X[\boldsymbol{h}](\boldsymbol{u}_r)$ of $L_X[\mathbf{A}]$ and $L_X[\boldsymbol{h}](\boldsymbol{u}_r)$ respectively and consider the following ordinary differential equation whose solution $\tilde{\Delta}_r$ provides an estimation of an upper bound of $\Delta_r$:

$$ \left\{ \begin{array}{l} \tilde{\Delta}'_r(t, \boldsymbol{\xi}) = \tilde{L}_X[\mathbf{A}]\tilde{\Delta}_r(t, \boldsymbol{\xi}) + \tilde{L}_X[\boldsymbol{h}](\boldsymbol{u}_r)\tilde{\Delta}_r(t, \boldsymbol{\xi}) + \|\tilde{\boldsymbol{r}}(t, \boldsymbol{\xi})\|_X, \\ \tilde{\Delta}_r(0, \boldsymbol{\xi}) = \|\boldsymbol{e}_r(0, \boldsymbol{\xi})\|_X. \end{array} \right. \quad (20) $$

Here, $\|\tilde{\boldsymbol{r}}(t, \boldsymbol{\xi})\|_X$ denotes an approximation of $\|\boldsymbol{r}(t, \boldsymbol{\xi})\|_X$ where $\boldsymbol{h}(\boldsymbol{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi})$ is replaced by an approximation $\tilde{\boldsymbol{h}}(\boldsymbol{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi})$ obtained by an empirical interpolation method (see Subsection 4.2.1).
Equation (20) is then solved using the following semi-implicit scheme (consistent with (16)):

$$ \tilde{\Delta}_r^{k+1}(\xi) = (1 - \delta t \tilde{L}_X[\mathbf{A}^{k+1}(\xi)])^{-1} (\tilde{\Delta}_r^k(\xi) + \delta t \tilde{L}_X[\mathbf{h}^k](\mathbf{u}_r^k(\xi)) + \delta t \|\tilde{\mathbf{r}}^k(\xi)\|_X). \quad (21) $$

with

$$ \tilde{\mathbf{r}}^k(\xi) = \Pi_{X_r^{k+1} \perp} \left( \frac{\delta \mathbf{V}_r^k}{\delta t} \alpha_r^k(\xi) - \mathbf{A}^{k+1}(\xi) \mathbf{u}_r^{k+1} - \tilde{\mathbf{h}}^k(\mathbf{u}_r^k(\xi), \xi) - \mathbf{g}^k(\xi) \right), \quad (22) $$

where $\frac{\delta \mathbf{V}_r^k}{\delta t}$ takes into account the possible time dependency of the reduced basis (this term is equal to zero for time-independent reduced spaces).

## 4.2 Offline/Online implementation

A time and parameter-dependent function $a(t, \xi)$ with values in some vector space $V$ is said to admit a time-dependent affine representation if

$$ a(t, \xi) = \sum_{i=1}^{Q_A} \theta_a^i(t, \xi) a^i(t), $$

with $\theta_a^i(t, \xi) \in \mathbb{R}$ and $a^i(t) \in V$. For a function discretized on a time grid $\{t^k\}_{k=0}^K$, that means $a^k(\xi) = \sum_{i=1}^{Q_A} \theta_a^{i,k}(\xi) a^{i,k}$, for $0 \le k \le K$. For the sake of presentation, we keep the notation $a(t, \xi)$ for time-dependent functions, even if $t$ only takes a finite set of values $t^k$, $0 \le k \le K$. Then, we assume that $\mathbf{A}(t, \xi)$ and $\mathbf{g}(t, \xi)$ admit time-dependent affine representations

$$ \mathbf{A}(t, \xi) = \sum_{i=1}^{Q_A} \theta_{\mathbf{A}}^i(t, \xi) \mathbf{A}^i(t), \quad \mathbf{g}(t, \xi) = \sum_{i=1}^{Q_g} \theta_{\mathbf{g}}^i(t, \xi) \mathbf{g}^i(t). $$

### 4.2.1 Solution of the reduced dynamical system

The solution of the reduced dynamical system (18) requires an efficient evaluation of

$$ \mathbf{A}_r(t^k, \xi) = \mathbf{V}_r(t^k)^T \mathbf{A}(t^k, \xi) \mathbf{V}_r(t^k), \quad \mathbf{g}_r(t^k, \xi) = \mathbf{V}_r(t^k)^T \mathbf{g}(t^k, \xi), $$

and

$$ \mathbf{h}_r^{k+1}(\alpha_r(t^k, \xi), t^k, \xi) = \mathbf{V}_r(t^{k+1})^T \mathbf{h}(\alpha_r(t^k, \xi), t^k, \xi). $$

In the offline phase, after computing the reduced space $X_r(t)$ and associated basis $\mathbf{V}_r(t)$, the reduced matrices $\mathbf{A}_r^i(t) = \mathbf{V}_r(t)^T \mathbf{A}^i(t) \mathbf{V}_r(t)$ and reduced vectors $\mathbf{g}_r^i(t) = \mathbf{V}_r(t)^T \mathbf{g}_r^i(t)$ can be precomputed. Then in the online phase, the reduced matrix $\mathbf{A}_r(t, \xi)$ and reduced
Equation (20) is then solved using the following semi-implicit scheme (consistent with (16)):

$$ \tilde{\Delta}_r^{k+1}(\xi) = (1 - \delta t \tilde{L}_X[\mathbf{A}^{k+1}(\xi)])^{-1} (\tilde{\Delta}_r^k(\xi) + \delta t \tilde{L}_X[\mathbf{h}^k](\mathbf{u}_r^k(\xi)) + \delta t \|\tilde{\mathbf{r}}^k(\xi)\|_X). \quad (21) $$

with

$$ \tilde{\mathbf{r}}^k(\xi) = \Pi_{X_r^{k+1}\perp} \left( \frac{\delta \mathbf{V}_r^k}{\delta t} \alpha_r^k(\xi) - \mathbf{A}^{k+1}(\xi) \mathbf{u}_r^{k+1} - \tilde{\mathbf{h}}^k(\mathbf{u}_r^k(\xi), \xi) - \mathbf{g}^k(\xi) \right), \quad (22) $$

where $\frac{\delta \mathbf{V}_r^k}{\delta t}$ takes into account the possible time dependency of the reduced basis (this term is equal to zero for time-independent reduced spaces).

## 4.2 Offline/Online implementation

A time and parameter-dependent function $a(t, \xi)$ with values in some vector space $V$ is said to admit a time-dependent affine representation if

$$ a(t, \xi) = \sum_{i=1}^{Q_A} \theta_a^i(t, \xi) a^i(t), $$

with $\theta_a^i(t, \xi) \in \mathbb{R}$ and $a^i(t) \in V$. For a function discretized on a time grid $\{t^k\}_{k=0}^K$, that means $a^k(\xi) = \sum_{i=1}^{Q_A} \theta_a^{i,k}(\xi) a^{i,k}$, for $0 \le k \le K$. For the sake of presentation, we keep the notation $a(t, \xi)$ for time-dependent functions, even if $t$ only takes a finite set of values $t^k$, $0 \le k \le K$. Then, we assume that $\mathbf{A}(t, \xi)$ and $\mathbf{g}(t, \xi)$ admit time-dependent affine representations

$$ \mathbf{A}(t, \xi) = \sum_{i=1}^{Q_A} \theta_{\mathbf{A}}^i(t, \xi) \mathbf{A}^i(t), \quad \mathbf{g}(t, \xi) = \sum_{i=1}^{Q_g} \theta_{\mathbf{g}}^i(t, \xi) \mathbf{g}^i(t). $$

### 4.2.1 Solution of the reduced dynamical system

The solution of the reduced dynamical system (18) requires an efficient evaluation of

$$ \mathbf{A}_r(t^k, \xi) = \mathbf{V}_r(t^k)^T \mathbf{A}(t^k, \xi) \mathbf{V}_r(t^k), \quad \mathbf{g}_r(t^k, \xi) = \mathbf{V}_r(t^k)^T \mathbf{g}(t^k, \xi), $$

and

$$ \mathbf{h}_r^{k+1}(\alpha_r(t^k, \xi), t^k, \xi) = \mathbf{V}_r(t^{k+1})^T \mathbf{h}(\alpha_r(t^k, \xi), t^k, \xi). $$

In the offline phase, after computing the reduced space $X_r(t)$ and associated basis $\mathbf{V}_r(t)$, the reduced matrices $\mathbf{A}_r^i(t) = \mathbf{V}_r(t)^T \mathbf{A}^i(t) \mathbf{V}_r(t)$ and reduced vectors $\mathbf{g}_r^i(t) = \mathbf{V}_r(t)^T \mathbf{g}_r^i(t)$ can be precomputed. Then in the online phase, the reduced matrix $\mathbf{A}_r(t, \xi)$ and reduced
vector $g_r(t, \xi)$ can be evaluated for any parameters value and any time with a complexity independent of $d$ using

$$A_r(t, \xi) = \sum_{i=1}^{Q_A} \theta_A^i(t, \xi) A_r^i(t), \quad g_r(t) = \sum_{i=1}^{Q_g} \theta_g^i(t, \xi) g_r^i(t).$$

An efficient way of treating the nonlinear term is to approximate $h(\mathbf{u}_r(t, \xi), t, \xi)$ using the Empirical Interpolation Method (EIM). Here, we briefly recall the principle of this method (see [3] for a detailed presentation). For a given basis $\{h_1, \dots, h_m\}$ in $\mathbb{R}^d$, with $m \le d$, the EIM approximation $\tilde{h}(\mathbf{u}_r(t, \xi), t, \xi)$ of order $m$ of $h(\mathbf{u}_r(t, \xi), t, \xi)$ is given by

$$\tilde{h}(\mathbf{u}_r(t, \xi), t, \xi) = H_m(\mathbf{P}_m^T \mathbf{H}_m)^{-1} \mathbf{P}_m^T h(\mathbf{u}_r(t, \xi), t, \xi), \quad (23)$$

where $\mathbf{P}_m = [e_{j_1}, \dots, e_{j_m}] \in \mathbb{R}^{m \times d}$, with $e_j$ the $j$-th unit vector in $\mathbb{R}^d$, and $H_m = [h_1, \dots, h_m] \in \mathbb{R}^{d \times m}$. Here, the vectors $h_1, \dots, h_m$ are chosen as linear combinations of columns $(i_1, \dots, i_m)$ of the matrix of snapshots $h(\mathbf{u}(t^k, \xi^i), t^k, \xi^i)$ already computed during the greedy selection of $\{\xi^1, \dots, \xi^r\}$. The set of indices $\mathcal{I}_m = \{(i_l, j_l), 1 \le l \le m\}$ is obtained with a greedy algorithm (see [3, §4.4]) which ensures that the matrix $\mathbf{P}_m^T \mathbf{H}_m$ is non singular and well-conditioned. An approximation of the reduced flux $V_r(t^{k+1})^T h(\mathbf{u}_r(t^k, \xi), t^k, \xi)$ is then given by

$$V_r(t^{k+1})^T U_m \mathbf{P}_m^T h(\mathbf{u}_r(t^k, \xi), t^k, \xi),$$

where $U_m = H_m(\mathbf{P}_m^T \mathbf{H}_m)^{-1} \in \mathbb{R}^{m \times d}$. The matrices $V_r(t^{k+1})^T U_m \in \mathbb{R}^{r \times m}$ can be pre-computed during the offline phase and stored to be later multiplied by the vectors $\mathbf{P}_m^T h(\mathbf{u}_r(t^k, \xi), t^k, \xi)$ which contains the components $(j_1, \dots, j_m)$ of the flux $h(\mathbf{u}_r(t^k, \xi), t^k, \xi)$. In practice, the EIM algorithm is stopped when the precision $\varepsilon$ is reached, i.e. when

$$\max_{k,i} \|h(\mathbf{u}_r(t^k, \xi^i), t^k, \xi^i) - \tilde{h}(\mathbf{u}_r(t^k, \xi^i), t^k, \xi^i)\|_X < \varepsilon. \quad (24)$$

In the following computations, the flux approximation defined by (23) is used for the solution of the discrete reduced dynamical system (18). Here, we assume that $\varepsilon$ is chosen small enough for the interpolation error between $\tilde{h}$ and $h$ to be neglected in our error analysis.

**Remark 5.** As in [34], this method can be reduced to the adaptive selection of the indices $(i_1, \dots, i_m)$ only, where the set of vectors $\{h_1, \dots, h_m\}$ is obtained by a POD of the matrix containing the snapshots of the discrete flux [8].

**Remark 6.** In the present work, we have applied an EIM to construct an approximation of the nonlinear flux $h(\mathbf{u}_r(t, \xi), t, \xi)$ on a time-independent basis $\{h_1, \dots, h_m\}$ in $\mathbb{R}^d$. An alternative approach, not considered here, would consist in approximating the flux on a time-dependent basis $\{h_1(t), \dots, h_m(t)\}$ constructed using an algorithm similar to the proposed T-greedy algorithm, using a global in time error estimate providing indices $\mathcal{I}_m$.
vector $g_r(t, \xi)$ can be evaluated for any parameters value and any time with a complexity independent of $d$ using

$$A_r(t, \xi) = \sum_{i=1}^{Q_A} \theta_A^i(t, \xi) A_r^i(t), \quad g_r(t) = \sum_{i=1}^{Q_g} \theta_g^i(t, \xi) g_r^i(t).$$

An efficient way of treating the nonlinear term is to approximate $h(\mathbf{u}_r(t, \xi), t, \xi)$ using the Empirical Interpolation Method (EIM). Here, we briefly recall the principle of this method (see [3] for a detailed presentation). For a given basis $\{\mathbf{h}_1, \dots, \mathbf{h}_m\}$ in $\mathbb{R}^d$, with $m \le d$, the EIM approximation $\tilde{h}(\mathbf{u}_r(t, \xi), t, \xi)$ of order $m$ of $h(\mathbf{u}_r(t, \xi), t, \xi)$ is given by

$$\tilde{h}(\mathbf{u}_r(t, \xi), t, \xi) = H_m(\mathbf{P}_m^T \mathbf{H}_m)^{-1} \mathbf{P}_m^T h(\mathbf{u}_r(t, \xi), t, \xi), \quad (23)$$

where $\mathbf{P}_m = [e_{j_1}, \dots, e_{j_m}] \in \mathbb{R}^{m \times d}$, with $e_j$ the $j$-th unit vector in $\mathbb{R}^d$, and $\mathbf{H}_m = [\mathbf{h}_1, \dots, \mathbf{h}_m] \in \mathbb{R}^{d \times m}$. Here, the vectors $\mathbf{h}_1, \dots, \mathbf{h}_m$ are chosen as linear combinations of columns $(i_1, \dots, i_m)$ of the matrix of snapshots $h(\mathbf{u}(t^k, \xi^i), t^k, \xi^i)$ already computed during the greedy selection of $\{\xi^1, \dots, \xi^r\}$. The set of indices $\mathcal{I}_m = \{(i_l, j_l), 1 \le l \le m\}$ is obtained with a greedy algorithm (see [3, §4.4]) which ensures that the matrix $\mathbf{P}_m^T \mathbf{H}_m$ is non singular and well-conditioned. An approximation of the reduced flux $\mathbf{V}_r(t^{k+1})^T h(\mathbf{u}_r(t^k, \xi), t^k, \xi)$ is then given by

$$\mathbf{V}_r(t^{k+1})^T \mathbf{U}_m \mathbf{P}_m^T h(\mathbf{u}_r(t^k, \xi), t^k, \xi),$$

where $\mathbf{U}_m = \mathbf{H}_m(\mathbf{P}_m^T \mathbf{H}_m)^{-1} \in \mathbb{R}^{m \times d}$. The matrices $\mathbf{V}_r(t^{k+1})^T \mathbf{U}_m \in \mathbb{R}^{r \times m}$ can be pre-computed during the offline phase and stored to be later multiplied by the vectors $\mathbf{P}_m^T h(\mathbf{u}_r(t^k, \xi), t^k, \xi)$ which contains the components $(j_1, \dots, j_m)$ of the flux $h(\mathbf{u}_r(t^k, \xi), t^k, \xi)$. In practice, the EIM algorithm is stopped when the precision $\varepsilon$ is reached, i.e. when

$$\max_{k,i} \|h(\mathbf{u}_r(t^k, \xi^i), t^k, \xi^i) - \tilde{h}(\mathbf{u}_r(t^k, \xi^i), t^k, \xi^i)\|_X < \varepsilon. \quad (24)$$

In the following computations, the flux approximation defined by (23) is used for the solution of the discrete reduced dynamical system (18). Here, we assume that $\varepsilon$ is chosen small enough for the interpolation error between $\tilde{h}$ and $h$ to be neglected in our error analysis.

**Remark 5.** As in [34], this method can be reduced to the adaptive selection of the indices $(i_1, \dots, i_m)$ only, where the set of vectors $\{\mathbf{h}_1, \dots, \mathbf{h}_m\}$ is obtained by a POD of the matrix containing the snapshots of the discrete flux [8].

**Remark 6.** In the present work, we have applied an EIM to construct an approximation of the nonlinear flux $h(\mathbf{u}_r(t, \xi), t, \xi)$ on a time-independent basis $\{\mathbf{h}_1, \dots, \mathbf{h}_m\}$ in $\mathbb{R}^d$. An alternative approach, not considered here, would consist in approximating the flux on a time-dependent basis $\{\mathbf{h}_1(t), \dots, \mathbf{h}_m(t)\}$ constructed using an algorithm similar to the proposed T-greedy algorithm, using a global in time error estimate providing indices $\mathcal{I}_m$.
### 4.2.2 Evaluation of the error estimate

Now we detail practical aspects for the evaluation of the a posteriori error estimate $\tilde{\Delta}_r$ using the integration scheme (21). For an efficient online computation of the error estimate, we need to evaluate with a complexity independent of $d$ the term $\|\tilde{\mathbf{r}}^k(\boldsymbol{\xi})\|_X$ and the estimations $\tilde{L}_X[\mathbf{A}]$ and $\tilde{L}_X[\mathbf{h}](\mathbf{u}_r(t, \boldsymbol{\xi}))$ of the logarithmic Lipschitz constants $L_X[\mathbf{A}]$ and $L_X[\mathbf{h}](\mathbf{u}_r(t, \boldsymbol{\xi}))$.

We first address the estimation of logarithmic Lipschitz constants. Using the time-dependent affine representation of $\mathbf{A}$, we have

$$L_X[\mathbf{A}] = \sup_{0 \neq \mathbf{x} \in X} \sum_{i=1}^{Q_A} \theta_{\mathbf{A}}^i(t, \boldsymbol{\xi}) \frac{\langle \mathbf{x}, \mathbf{A}^i(t) \mathbf{x} \rangle_X}{\|\mathbf{x}\|_X^2} \le \sum_{i=1}^{Q_A} |\theta_{\mathbf{A}}^i(t, \boldsymbol{\xi})| L_X[\mathbf{A}^i(t)] := \tilde{L}_X[\mathbf{A}]. \quad (25)$$

For the local logarithmic constant of $\mathbf{h}$ at $\mathbf{u}_r$, we consider the first-order linearization of the flux

$$\mathbf{h}(\mathbf{v}, t, \boldsymbol{\xi}) \approx \mathbf{h}(\mathbf{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi}) + \nabla \mathbf{h}(\mathbf{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi})(\mathbf{v} - \mathbf{u}_r(t, \boldsymbol{\xi})),$$

where $\nabla \mathbf{h}(\mathbf{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi}) \in \mathbb{R}^{d \times d}$ denotes the gradient of $\mathbf{h}(\cdot, t, \boldsymbol{\xi})$ at $\mathbf{u}_r(t, \boldsymbol{\xi})$, and the corresponding approximation

$$L_X[\mathbf{h}](\mathbf{u}_r(t, \boldsymbol{\xi})) \approx L_X[\nabla \mathbf{h}(\mathbf{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi})].$$

Then, a first order approximation of $L_X[\mathbf{h}](\mathbf{u}_r)$ is obtained by computing the largest eigenvalue of the symmetric part of $\nabla \mathbf{h}(\mathbf{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi})$, which is a problem with complexity depending on $d$. In [34], in order to obtain a complexity independent on $d$, the authors use a partial similarity transformation of the matrix $\nabla \mathbf{h}$ based on POD (combined with a matrix-DEIM approximation) which preserves the largest eigenvalue of the symmetric part of $\nabla \mathbf{h}$. Here, we adopt a simpler strategy which consists in interpolating the Lipschitz constants of the matrices $\nabla \mathbf{h}$ already computed during the solution of the full-order dynamical system in the offline step. We introduce a very simple nearest neighbor interpolation

$$L_X[\nabla \mathbf{h}(\mathbf{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi})] \approx \sum_{i=1}^{r} \gamma^i(\boldsymbol{\xi}) L_X[\nabla \mathbf{h}](\mathbf{u}(t, \boldsymbol{\xi}^i), t, \boldsymbol{\xi}^i) := \tilde{L}_X[\mathbf{h}(\mathbf{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi})], \quad (26)$$

where the $\{\boldsymbol{\xi}^i\}_{i=1}^r$ are the parameters values selected during the greedy procedure, and where

$$\gamma^i(\boldsymbol{\xi}) = \begin{cases} 1 & \text{if } d(\boldsymbol{\xi}, \boldsymbol{\xi}^i) = \min_{1 \le j \le r} d(\boldsymbol{\xi}, \boldsymbol{\xi}^j), \\ 0 & \text{otherwise,} \end{cases}$$

with $d(\cdot, \cdot)$ a metric on $\Xi$ (typically the Euclidian metric). Note that when using the interpolation method, $\tilde{L}_X[\mathbf{h}](\mathbf{u}_r)$ is not necessarily an upper bound of $L_X[\mathbf{h}](\mathbf{u}_r)$ and therefore, the solution of (13) with $L_X[\mathbf{h}](\mathbf{u}_r)$ replaced by $\tilde{L}_X[\mathbf{h}](\mathbf{u}_r)$ may not provide a certified
### 4.2.2 Evaluation of the error estimate

Now we detail practical aspects for the evaluation of the a posteriori error estimate $\tilde{\Delta}_r$ using the integration scheme (21). For an efficient online computation of the error estimate, we need to evaluate with a complexity independent of $d$ the term $\|\tilde{\mathbf{r}}^k(\boldsymbol{\xi})\|_X$ and the estimations $\tilde{L}_X[\mathbf{A}]$ and $\tilde{L}_X[\mathbf{h}](\mathbf{u}_r(t, \boldsymbol{\xi}))$ of the logarithmic Lipschitz constants $L_X[\mathbf{A}]$ and $L_X[\mathbf{h}](\mathbf{u}_r(t, \boldsymbol{\xi}))$.

We first address the estimation of logarithmic Lipschitz constants. Using the time-dependent affine representation of $\mathbf{A}$, we have

$$L_X[\mathbf{A}] = \sup_{0 \neq \mathbf{x} \in X} \sum_{i=1}^{Q_A} \theta_A^i(t, \boldsymbol{\xi}) \frac{\langle \mathbf{x}, \mathbf{A}^i(t) \mathbf{x} \rangle_X}{\|\mathbf{x}\|_X^2} \leq \sum_{i=1}^{Q_A} |\theta_A^i(t, \boldsymbol{\xi})| L_X[\mathbf{A}^i(t)] := \tilde{L}_X[\mathbf{A}]. \quad (25)$$

For the local logarithmic constant of $\mathbf{h}$ at $\mathbf{u}_r$, we consider the first-order linearization of the flux

$$\mathbf{h}(\mathbf{v}, t, \boldsymbol{\xi}) \approx \mathbf{h}(\mathbf{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi}) + \nabla \mathbf{h}(\mathbf{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi})(\mathbf{v} - \mathbf{u}_r(t, \boldsymbol{\xi})),$$

where $\nabla \mathbf{h}(\mathbf{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi}) \in \mathbb{R}^{d \times d}$ denotes the gradient of $\mathbf{h}(\cdot, t, \boldsymbol{\xi})$ at $\mathbf{u}_r(t, \boldsymbol{\xi})$, and the corresponding approximation

$$L_X[\mathbf{h}](\mathbf{u}_r(t, \boldsymbol{\xi})) \approx L_X[\nabla \mathbf{h}(\mathbf{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi})].$$

Then, a first order approximation of $L_X[\mathbf{h}](\mathbf{u}_r)$ is obtained by computing the largest eigenvalue of the symmetric part of $\nabla \mathbf{h}(\mathbf{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi})$, which is a problem with complexity depending on $d$. In [34], in order to obtain a complexity independent on $d$, the authors use a partial similarity transformation of the matrix $\nabla \mathbf{h}$ based on POD (combined with a matrix-DEIM approximation) which preserves the largest eigenvalue of the symmetric part of $\nabla \mathbf{h}$. Here, we adopt a simpler strategy which consists in interpolating the Lipschitz constants of the matrices $\nabla \mathbf{h}$ already computed during the solution of the full-order dynamical system in the offline step. We introduce a very simple nearest neighbor interpolation

$$L_X[\nabla \mathbf{h}(\mathbf{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi})] \approx \sum_{i=1}^{r} \gamma^i(\boldsymbol{\xi}) L_X[\nabla \mathbf{h}](\mathbf{u}(t, \boldsymbol{\xi}^i), t, \boldsymbol{\xi}^i) := \tilde{L}_X[\mathbf{h}(\mathbf{u}_r(t, \boldsymbol{\xi}), t, \boldsymbol{\xi})], \quad (26)$$

where the $\{\boldsymbol{\xi}^i\}_{i=1}^r$ are the parameters values selected during the greedy procedure, and where

$$\gamma^i(\boldsymbol{\xi}) = \begin{cases} 1 & \text{if } d(\boldsymbol{\xi}, \boldsymbol{\xi}^i) = \min_{1 \leq j \leq r} d(\boldsymbol{\xi}, \boldsymbol{\xi}^j), \\ 0 & \text{otherwise,} \end{cases}$$

with $d(\cdot, \cdot)$ a metric on $\Xi$ (typically the Euclidian metric). Note that when using the interpolation method, $\tilde{L}_X[\mathbf{h}](\mathbf{u}_r)$ is not necessarily an upper bound of $L_X[\mathbf{h}](\mathbf{u}_r)$ and therefore, the solution of (13) with $L_X[\mathbf{h}](\mathbf{u}_r)$ replaced by $\tilde{L}_X[\mathbf{h}](\mathbf{u}_r)$ may not provide a certified
error bound. However, as it will be seen in the numerical experiments, it provides a very good error indicator in practice.

Now we detail the computation of $||\tilde{\mathbf{r}}^k(\boldsymbol{\xi})||_X$ in an offline/online strategy, where $\tilde{\mathbf{r}}^k(\boldsymbol{\xi})$ is given by (22) with $\tilde{\mathbf{h}}$ the approximation of $\mathbf{h}$ defined by (23). In what follows, the indices $i$ and $j$ take values in $\{1, \dots, Q_A\}$ or $\{1, \dots, Q_g\}$. In the offline phase, we pre-compute the following matrices and vectors required for the online phase:

$$ 
\begin{aligned}
\mathbf{K}_{ij}^{1,k} &= \mathbf{V}_r^{kT} \mathbf{A}^{i,kT} \Pi_{X_r^{k\perp}} \mathbf{A}^{j,k} \mathbf{V}_r^k, & \mathbf{K}_i^{2,k} &= \frac{\delta \mathbf{V}_r^{kT}}{\delta t} \Pi_{X_r^{k+1\perp}} \mathbf{A}^{i,k+1} \mathbf{V}_r^{k+1}, \\
\mathbf{K}^{3,k} &= \frac{\delta \mathbf{V}_r^{kT}}{\delta t} \Pi_{X_r^{k+1\perp}} \frac{\delta \mathbf{V}_r^k}{\delta t}, & \mathbf{K}^{4,k} &= \frac{\delta \mathbf{V}_r^{kT}}{\delta t} \Pi_{X_r^{k+1\perp}} \mathbf{U}_m, \\
\mathbf{K}_i^{5,k} &= \mathbf{V}_r^{kT} \mathbf{A}^{i,kT} \Pi_{X_r^{k\perp}} \mathbf{U}_m,
\end{aligned}
$$

and

$$ 
\begin{aligned}
\mathbf{b}_i^{1,k} &= \frac{\delta \mathbf{V}_r^{kT}}{\delta t} \Pi_{X_r^{k+1\perp}} \mathbf{g}^{i,k}, & \mathbf{b}_{ij}^{2,k} &= \mathbf{V}_r^{k+1,T} \mathbf{A}^{i,k+1,T} \Pi_{X_r^{k+1\perp}} \mathbf{g}^{j,k}, \\
\mathbf{b}_i^{3,k} &= \mathbf{U}_m^T \Pi_{X_r^{k+1\perp}} \mathbf{g}^{i,k}.
\end{aligned}
$$

Then, in the online phase we compute

$$ 
\begin{aligned}
||\tilde{\mathbf{r}}^k(\boldsymbol{\xi})||_X^2 &= \langle \boldsymbol{\alpha}_r^{k+1}(\boldsymbol{\xi}), \mathbf{M}_1^{k+1}(\boldsymbol{\xi})\boldsymbol{\alpha}_r^{k+1}(\boldsymbol{\xi}) \rangle_X + \langle \boldsymbol{\alpha}_r^k(\boldsymbol{\xi}), \mathbf{M}_2^k(\boldsymbol{\xi})\boldsymbol{\alpha}_r^{k+1}(\boldsymbol{\xi}) \rangle_X \\
&\quad + \langle \boldsymbol{\alpha}_r^k(\boldsymbol{\xi}), \mathbf{M}_3^k(\boldsymbol{\xi})\boldsymbol{\alpha}_r^k(\boldsymbol{\xi}) \rangle_X + \langle \boldsymbol{\alpha}_r^k(\boldsymbol{\xi}), \mathbf{M}_4^k(\boldsymbol{\xi})\mathbf{v}_1^k(\boldsymbol{\xi}) \rangle_X \\
&\quad + \langle \boldsymbol{\alpha}_r^{k+1}(\boldsymbol{\xi}), \mathbf{M}_5^{k+1}(\boldsymbol{\xi})\mathbf{v}_1^k(\boldsymbol{\xi}) \rangle_X + \langle \mathbf{v}_1^k(\boldsymbol{\xi}), \mathbf{M}_6^k\mathbf{v}_1^k(\boldsymbol{\xi}) \rangle_X \\
&\quad + \langle \mathbf{v}_1^k(\boldsymbol{\xi}), \mathbf{v}_2^k(\boldsymbol{\xi}) \rangle_X + \langle \boldsymbol{\alpha}_r^k(\boldsymbol{\xi}), \mathbf{v}_3(\boldsymbol{\xi}) \rangle_X + \langle \boldsymbol{\alpha}_r^{k+1}(\boldsymbol{\xi}), \mathbf{v}_4^k(\boldsymbol{\xi}) \rangle_X + b^k
\end{aligned}
$$

where the reduced quantities are defined by

$$ 
\begin{aligned}
\mathbf{M}_1^k(\boldsymbol{\xi}) &= \sum_{i,j=1}^{Q_A} \theta_A^{i,k}(\boldsymbol{\xi}) \theta_A^{j,k}(\boldsymbol{\xi}) \mathbf{K}_{ij}^{1,k}, & \mathbf{M}_2^k(\boldsymbol{\xi}) &= -2 \sum_{i=1}^{Q_A} \theta_A^{i,k+1}(\boldsymbol{\xi}) \mathbf{K}_i^{2,k}, & \mathbf{M}_3^k(\boldsymbol{\xi}) &= \mathbf{K}^{3,k}, \\
\mathbf{M}_4^k(\boldsymbol{\xi}) &= -2\mathbf{K}^{4,k}, & \mathbf{M}_5^k(\boldsymbol{\xi}) &= 2 \sum_{i=1}^{Q_A} \theta_A^{i,k}(\boldsymbol{\xi}) \mathbf{K}_i^{5,k}, & \mathbf{M}_6^k &= \mathbf{U}_m^T \Pi_{X_r^{k\perp}} \mathbf{U}_m,
\end{aligned}
$$

$$ 
\begin{aligned}
\mathbf{v}_1^k(\boldsymbol{\xi}) &= \mathbf{P}_m^T \mathbf{h}^k(\mathbf{u}_r^k(\boldsymbol{\xi}), \boldsymbol{\xi}), & \mathbf{v}_2^k(\boldsymbol{\xi}) &= 2 \sum_{i=1}^{Q_g} \theta_g^{i,k}(\boldsymbol{\xi}) \mathbf{b}_i^{3,k}, \\
\mathbf{v}_3^k(\boldsymbol{\xi}) &= -2 \sum_{i=1}^{Q_g} \theta_g^{i,k}(\boldsymbol{\xi}) \mathbf{b}_i^{1,k}, & \mathbf{v}_4^k(\boldsymbol{\xi}) &= 2 \sum_{i=1}^{Q_A} \sum_{j=1}^{Q_g} \theta_A^{i,k+1}(\boldsymbol{\xi}) \theta_g^{j,k}(\boldsymbol{\xi}) \mathbf{b}_{ij}^{2,k},
\end{aligned}
$$
error bound. However, as it will be seen in the numerical experiments, it provides a very good error indicator in practice.

Now we detail the computation of $\|\tilde{\mathbf{r}}^k(\boldsymbol{\xi})\|_X$ in an offline/online strategy, where $\tilde{\mathbf{r}}^k(\boldsymbol{\xi})$ is given by (22) with $\tilde{\mathbf{h}}$ the approximation of $\mathbf{h}$ defined by (23). In what follows, the indices $i$ and $j$ take values in $\{1, \dots, Q_A\}$ or $\{1, \dots, Q_g\}$. In the offline phase, we pre-compute the following matrices and vectors required for the online phase:

$$ 
\begin{aligned}
\mathbf{K}_{ij}^{1,k} &= \mathbf{V}_r^{kT} \mathbf{A}^{i,kT} \Pi_{X_r^{k\perp}} \mathbf{A}^{j,k} \mathbf{V}_r^k, & \mathbf{K}_i^{2,k} &= \frac{\delta \mathbf{V}_r^{kT}}{\delta t} \Pi_{X_r^{k+1\perp}} \mathbf{A}^{i,k+1} \mathbf{V}_r^{k+1}, \\
\mathbf{K}^{3,k} &= \frac{\delta \mathbf{V}_r^{kT}}{\delta t} \Pi_{X_r^{k+1\perp}} \frac{\delta \mathbf{V}_r^k}{\delta t}, & \mathbf{K}^{4,k} &= \frac{\delta \mathbf{V}_r^{kT}}{\delta t} \Pi_{X_r^{k+1\perp}} \mathbf{U}_m, \\
\mathbf{K}_i^{5,k} &= \mathbf{V}_r^{kT} \mathbf{A}^{i,kT} \Pi_{X_r^{k\perp}} \mathbf{U}_m,
\end{aligned}
$$

and

$$ 
\begin{aligned}
\mathbf{b}_i^{1,k} &= \frac{\delta \mathbf{V}_r^{kT}}{\delta t} \Pi_{X_r^{k+1\perp}} \mathbf{g}^{i,k}, & \mathbf{b}_{ij}^{2,k} &= \mathbf{V}_r^{k+1,T} \mathbf{A}^{i,k+1,T} \Pi_{X_r^{k+1\perp}} \mathbf{g}^{j,k}, \\
\mathbf{b}_i^{3,k} &= \mathbf{U}_m^T \Pi_{X_r^{k+1\perp}} \mathbf{g}^{i,k}.
\end{aligned}
$$

Then, in the online phase we compute

$$ 
\begin{aligned}
\|\tilde{\mathbf{r}}^k(\boldsymbol{\xi})\|_X^2 &= \langle \boldsymbol{\alpha}_r^{k+1}(\boldsymbol{\xi}), \mathbf{M}_1^{k+1}(\boldsymbol{\xi})\boldsymbol{\alpha}_r^{k+1}(\boldsymbol{\xi}) \rangle_X + \langle \boldsymbol{\alpha}_r^k(\boldsymbol{\xi}), \mathbf{M}_2^k(\boldsymbol{\xi})\boldsymbol{\alpha}_r^{k+1}(\boldsymbol{\xi}) \rangle_X \\
&\quad + \langle \boldsymbol{\alpha}_r^k(\boldsymbol{\xi}), \mathbf{M}_3^k(\boldsymbol{\xi})\boldsymbol{\alpha}_r^k(\boldsymbol{\xi}) \rangle_X + \langle \boldsymbol{\alpha}_r^k(\boldsymbol{\xi}), \mathbf{M}_4^k(\boldsymbol{\xi})\mathbf{v}_1^k(\boldsymbol{\xi}) \rangle_X \\
&\quad + \langle \boldsymbol{\alpha}_r^{k+1}(\boldsymbol{\xi}), \mathbf{M}_5^{k+1}(\boldsymbol{\xi})\mathbf{v}_1^k(\boldsymbol{\xi}) \rangle_X + \langle \mathbf{v}_1^k(\boldsymbol{\xi}), \mathbf{M}_6^k\mathbf{v}_1^k(\boldsymbol{\xi}) \rangle_X \\
&\quad + \langle \mathbf{v}_1^k(\boldsymbol{\xi}), \mathbf{v}_2^k(\boldsymbol{\xi}) \rangle_X + \langle \boldsymbol{\alpha}_r^k(\boldsymbol{\xi}), \mathbf{v}_3(\boldsymbol{\xi}) \rangle_X + \langle \boldsymbol{\alpha}_r^{k+1}(\boldsymbol{\xi}), \mathbf{v}_4^k(\boldsymbol{\xi}) \rangle_X + b^k
\end{aligned}
$$

where the reduced quantities are defined by

$$ 
\begin{aligned}
\mathbf{M}_1^k(\boldsymbol{\xi}) &= \sum_{i,j=1}^{Q_A} \theta_A^{i,k}(\boldsymbol{\xi}) \theta_A^{j,k}(\boldsymbol{\xi}) \mathbf{K}_{ij}^{1,k}, & \mathbf{M}_2^k(\boldsymbol{\xi}) &= -2 \sum_{i=1}^{Q_A} \theta_A^{i,k+1}(\boldsymbol{\xi}) \mathbf{K}_i^{2,k}, & \mathbf{M}_3^k(\boldsymbol{\xi}) &= \mathbf{K}^{3,k}, \\
\mathbf{M}_4^k(\boldsymbol{\xi}) &= -2\mathbf{K}^{4,k}, & \mathbf{M}_5^k(\boldsymbol{\xi}) &= 2 \sum_{i=1}^{Q_A} \theta_A^{i,k}(\boldsymbol{\xi}) \mathbf{K}_i^{5,k}, & \mathbf{M}_6^k &= \mathbf{U}_m^T \Pi_{X_r^{k\perp}} \mathbf{U}_m,
\end{aligned}
$$

$$ 
\begin{aligned}
\mathbf{v}_1^k(\boldsymbol{\xi}) &= \mathbf{P}_m^T \mathbf{h}^k(\mathbf{u}_r^k(\boldsymbol{\xi}), \boldsymbol{\xi}), & \mathbf{v}_2^k(\boldsymbol{\xi}) &= 2 \sum_{i=1}^{Q_g} \theta_g^{i,k}(\boldsymbol{\xi}) \mathbf{b}_i^{3,k}, \\
\mathbf{v}_3^k(\boldsymbol{\xi}) &= -2 \sum_{i=1}^{Q_g} \theta_g^{i,k}(\boldsymbol{\xi}) \mathbf{b}_i^{1,k}, & \mathbf{v}_4^k(\boldsymbol{\xi}) &= 2 \sum_{i=1}^{Q_A} \sum_{j=1}^{Q_g} \theta_A^{i,k+1}(\boldsymbol{\xi}) \theta_g^{j,k}(\boldsymbol{\xi}) \mathbf{b}_{ij}^{2,k},
\end{aligned}
$$
and

$$b^k = \sum_{i,j=1}^{Q_g} \theta_g^{i,k}(\xi) \theta_g^{i,k}(\xi) \langle g^{i,k}, g^{j,k} \rangle_X.$$

# 5 Numerical experiments

In this section, we present numerical results where we compare model order reduction methods using either time-independent or time-dependent reduced spaces. Also, we evaluate the effectiveness of the a posteriori error estimate derived in Subsection 2.2. We denote by MTI (resp. MTD) the model order reduction method using time-independent (resp. time-dependent) reduced spaces. If not mentioned, reduced spaces for MTI (resp. MTD) are constructed with a POD-greedy algorithm (resp. T-greedy algorithm) using a training set $\Xi_{train}$ whose size will be specified in each case.

We consider three test cases which are particular cases of the following nonlinear partial differential equation defined on a domain $\Omega$ which is a domain in $\mathbb{R}^s$ (with $s = 1, 2$) and a time interval $I = (0, T)$:

$$\frac{\partial}{\partial t}u + \nabla \cdot (c(u, \xi)u) + \mathbf{a}(x, \xi) \cdot \nabla u - \mu(\xi)\Delta u = g(x, t, \xi), \quad \text{on } \Omega \times I, \quad (27)$$

with appropriate boundary conditions and a parameter-independent initial condition $u^0(x)$. We denote $x = (x_1, \dots, x_s)^T$ the spatial variable and the corresponding spatial differential operators $\nabla = (\frac{\partial}{\partial x_1}, \dots, \frac{\partial}{\partial x_s})^T$ and $\Delta = \sum_{i=1}^s \frac{\partial^2}{\partial x_i^2}$. Here $\xi = (\xi^1, \dots, \xi^p)$ denotes a random vector with values in $\Xi$ and with independent components. The functions $c: X \times \Xi \to \mathbb{R}$ and $\mathbf{a}: \Omega \times \Xi \to \mathbb{R}^s$ will be specified in each test case. Finally, $\mu: \Xi \to \mathbb{R}$ is a parameter-dependent coefficient and $g: \Omega \times I \times \Xi \to \mathbb{R}$ is a given source term. We consider an approximation of $u$ obtained with an appropriate scheme (e.g. finite differences, finite element) depending on the test case. This yields a system of $d$ ordinary differential equations of the form (1). Here, $d = d(n)$ corresponds to the dimension of the discrete problem taking into account the boundary conditions (e.g. when using finite differences for one-dimensional problem, $d = n - 2$ for Dirichlet conditions, or $d = n - 1$ for periodic conditions). Finally, we denote the spatial approximation of the solution as follows $u(t) = (u_i(t))_{i=1}^d \in \mathbb{R}^d$, with approximations of both initial condition $u^0 \in \mathbb{R}^d$, and source term $g(t, \xi) \in \mathbb{R}^d$ that will be specified later. We assume that the time integration scheme is accurate enough for the error due to this approximation to be neglected.

Numerical experiments were conducted with an in-house code written in Matlab®.
and

$$b^k = \sum_{i,j=1}^{Q_g} \theta_g^{i,k}(\xi) \theta_g^{i,k}(\xi) \langle g^{i,k}, g^{j,k} \rangle_X.$$

# 5 Numerical experiments

In this section, we present numerical results where we compare model order reduction methods using either time-independent or time-dependent reduced spaces. Also, we evaluate the effectiveness of the a posteriori error estimate derived in Subsection 2.2. We denote by MTI (resp. MTD) the model order reduction method using time-independent (resp. time-dependent) reduced spaces. If not mentioned, reduced spaces for MTI (resp. MTD) are constructed with a POD-greedy algorithm (resp. T-greedy algorithm) using a training set $\Xi_{train}$ whose size will be specified in each case.

We consider three test cases which are particular cases of the following nonlinear partial differential equation defined on a domain $\Omega$ which is a domain in $\mathbb{R}^s$ (with $s = 1, 2$) and a time interval $I = (0, T)$:

$$\frac{\partial}{\partial t}u + \nabla \cdot (c(u, \xi)u) + a(x, \xi) \cdot \nabla u - \mu(\xi)\Delta u = g(x, t, \xi), \quad \text{on } \Omega \times I, \quad (27)$$

with appropriate boundary conditions and a parameter-independent initial condition $u^0(x)$. We denote $x = (x_1, \dots, x_s)^T$ the spatial variable and the corresponding spatial differential operators $\nabla = (\frac{\partial}{\partial x_1}, \dots, \frac{\partial}{\partial x_s})^T$ and $\Delta = \sum_{i=1}^s \frac{\partial^2}{\partial x_i^2}$. Here $\xi = (\xi^1, \dots, \xi^p)$ denotes a random vector with values in $\Xi$ and with independent components. The functions $c: X \times \Xi \to \mathbb{R}$ and $a: \Omega \times \Xi \to \mathbb{R}^s$ will be specified in each test case. Finally, $\mu: \Xi \to \mathbb{R}$ is a parameter-dependent coefficient and $g: \Omega \times I \times \Xi \to \mathbb{R}$ is a given source term. We consider an approximation of $u$ obtained with an appropriate scheme (e.g. finite differences, finite element) depending on the test case. This yields a system of $d$ ordinary differential equations of the form (1). Here, $d = d(n)$ corresponds to the dimension of the discrete problem taking into account the boundary conditions (e.g. when using finite differences for one-dimensional problem, $d = n - 2$ for Dirichlet conditions, or $d = n - 1$ for periodic conditions). Finally, we denote the spatial approximation of the solution as follows $u(t) = (u_i(t))_{i=1}^d \in \mathbb{R}^d$, with approximations of both initial condition $u^0 \in \mathbb{R}^d$, and source term $g(t, \xi) \in \mathbb{R}^d$ that will be specified later. We assume that the time integration scheme is accurate enough for the error due to this approximation to be neglected.

Numerical experiments were conducted with an in-house code written in Matlab®.
## 5.1 Test case 1

Let $\Omega = (0,1)$ and $I = (0,0.2)$. We consider an advection equation with $\mu(\xi) = 0$, $g(x, t, \xi) = 0$, $c(\mathbf{u}, \xi) = 0$ and $a(x, \xi) = a(\xi) = a_0 + a_1\xi^1$, with $\xi^1 \sim U(-1, 1)$, $a_0 = 1$ and $a_1 = 0.5$. The initial condition is a smooth function given by $u_{cont}^0(x) = \frac{1}{\sqrt{2\pi}} \exp(-\frac{x-0.6}{0.05}^2)$. We consider an approximation of $u$ obtained with an appropriate finite difference scheme over a uniform discretization $\{x_i\}_{i=1}^n$ of $\Omega$, where $x_i = a + i\delta x$, $\delta x = \frac{b-a}{n-1}$. This yields a system of $d$ ordinary differential equations of the form (1). We impose periodic boundary conditions and consider a finite difference upwind scheme on a uniform discretization with $n = 2001$ points. We use an explicit Euler time integration scheme

$$ \mathbf{u}^{k+1}(\xi) = (\mathbf{I}_d + \delta t \mathbf{C}(\xi)) \mathbf{u}^k(\xi), \qquad (28) $$

with initial condition $\mathbf{u}^0 = (u^0(x_i))_{i=1}^d$. Here $\delta t = 0.5 \frac{\delta x}{c_0+c_1}$, where $\mathbf{C}(\xi) \in \mathbb{R}^{d \times d}$ corresponds to the discrete advection operator with periodic boundary conditions.

**Deterministic case** We consider a deterministic problem with a fixed value $\xi = \xi_0 = 0.65$. Here, we compare the approximations obtained by projections on reduced spaces constructed in two different ways. In the first method (MTI), the reduced space $X_r$ is time-independent and generated by the first $r$ modes of the POD of the trajectory $t \mapsto \mathbf{u}(t, \xi_0)$, with $r \in \{1, 2, 5, 10, 20, 50, 100, 200\}$. In the second method (MTD), we consider the one-dimensional time-dependent space $X_r(t) = \text{span}\{\mathbf{u}(t, \xi_0)\}, r = 1$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Test case 1: comparison for discontinuous (left) and continuous (right) initial conditions at final times of the exact solution $\mathbf{u}$ and the approximation $\mathbf{u}_r$ computed with MTI for $r \in \{1, 5, 10, 20, 50, 100\}$.</figcaption></figure>
## 5.1 Test case 1

Let $\Omega = (0,1)$ and $I = (0,0.2)$. We consider an advection equation with $\mu(\xi) = 0$, $g(x, t, \xi) = 0$, $c(\mathbf{u}, \xi) = 0$ and $a(x, \xi) = a(\xi) = a_0 + a_1\xi^1$, with $\xi^1 \sim U(-1, 1)$, $a_0 = 1$ and $a_1 = 0.5$. The initial condition is a smooth function given by $u_{cont}^0(x) = \frac{1}{\sqrt{2\pi}} \exp(-\left(\frac{x-0.6}{0.05}\right)^2)$. We consider an approximation of $u$ obtained with an appropriate finite difference scheme over a uniform discretization $\{x_i\}_{i=1}^n$ of $\Omega$, where $x_i = a + i\delta x$, $\delta x = \frac{b-a}{n-1}$. This yields a system of $d$ ordinary differential equations of the form (1). We impose periodic boundary conditions and consider a finite difference upwind scheme on a uniform discretization with $n = 2001$ points. We use an explicit Euler time integration scheme

$$ \mathbf{u}^{k+1}(\xi) = (\mathbf{I}_d + \delta t \mathbf{C}(\xi)) \mathbf{u}^k(\xi), \qquad (28) $$

with initial condition $\mathbf{u}^0 = (u^0(x_i))_{i=1}^d$. Here $\delta t = 0.5 \frac{\delta x}{c_0+c_1}$, where $\mathbf{C}(\xi) \in \mathbb{R}^{d \times d}$ corresponds to the discrete advection operator with periodic boundary conditions.

**Deterministic case** We consider a deterministic problem with a fixed value $\xi = \xi_0 = 0.65$. Here, we compare the approximations obtained by projections on reduced spaces constructed in two different ways. In the first method (MTI), the reduced space $X_r$ is time-independent and generated by the first $r$ modes of the POD of the trajectory $t \mapsto \mathbf{u}(t, \xi_0)$, with $r \in \{1, 2, 5, 10, 20, 50, 100, 200\}$. In the second method (MTD), we consider the one-dimensional time-dependent space $X_r(t) = \text{span}\{\mathbf{u}(t, \xi_0)\}, r = 1$.

<figure><img src="image_3.png" /><figcaption>Figure 1: Test case 1: comparison for discontinuous (left) and continuous (right) initial conditions at final times of the exact solution $\mathbf{u}$ and the approximation $\mathbf{u}_r$ computed with MTI for $r \in \{1, 5, 10, 20, 50, 100\}$.</figcaption></figure>
<table><caption>Table 1: Test case 1: relative errors $E_2$ and $E_\infty$ obtained with the MTI with respect to $dim(X_r)$ for both continuous (left) and discontinuous (right) initial conditions.</caption><thead><tr><th colspan="3">Relative errors for u<sup>0</sup><sub>cont</sub></th><th colspan="3">Relative errors for u<sup>0</sup><sub>disc</sub></th></tr><tr><th>dim(X<sub>r</sub>)</th><th>E<sub>2</sub></th><th>E<sub>∞</sub></th><th>dim(X<sub>r</sub>)</th><th>E<sub>2</sub></th><th>E<sub>∞</sub></th></tr></thead><tbody><tr><td>1</td><td>0.027831</td><td>0.76058</td><td>1</td><td>0.022116</td><td>0.76058</td></tr><tr><td>2</td><td>0.024075</td><td>0.95713</td><td>2</td><td>0.0135</td><td>0.46427</td></tr><tr><td>5</td><td>0.0078461</td><td>0.26984</td><td>5</td><td>0.0060424</td><td>0.2078</td></tr><tr><td>10</td><td>0.00021853</td><td>0.0075156</td><td>10</td><td>0.0040044</td><td>0.13771</td></tr><tr><td>20</td><td>4.8616e-09</td><td>1.6719e-07</td><td>20</td><td>0.002043</td><td>0.07026</td></tr><tr><td>50</td><td>4.0924e-17</td><td>1.4074e-15</td><td>50</td><td>4.7139e-05</td><td>0.0016211</td></tr><tr><td>100</td><td>4.4353e-17</td><td>1.5253e-15</td><td>100</td><td>1.3891e-12</td><td>4.7771e-11</td></tr><tr><td>200</td><td>4.5623e-17</td><td>1.569e-15</td><td>200</td><td>6.3966e-17</td><td>2.1998e-15</td></tr></tbody></table>

As shown on Figure 1, the approximations obtained with MTI are satisfactory when the dimension of the reduced space $X_r$ is greater than 20. This is confirmed by the values of the relative errors $E_q = ||u_r - u||_{I,q}/||u||_{I,q}$ given in Table 1 for $q = 2$ and $q = \infty$, where $||u||_{I,q}$ is the natural norm on $L^q(I; X)$. In particular, the relative error in $L^2$-norm (resp $L^\infty$-norm) is of order $10^{-17}$ (respectively $10^{-15}$) for a space $X_r$ with a reduced dimension of 50. As expected, MTD, with a reduced space $X_1(t)$ of dimension 1 containing the exact solution, gives relative errors at the machine precision ($E_2 = 5.7 \ 10^{-17}$ and $E_\infty = 1.3 \ 10^{-15}$). For obtaining a very accurate precision, the method MTI requires a reduced space with rather high dimension, thus highlighting the limits of the POD method for transport equations. Now, we illustrate the impact of the smoothness of the solution. For this purpose, we consider the same test case with a discontinuous initial condition given by $u_{disc}^0(x) = \mathbf{1}_{[0.1,0.9]}(x) \cdot (\lfloor 3x \rfloor + \sin(10x))^2$. The approximate solution is plotted on Figure 1 and the relative errors between the reduced approximation and the exact solution are summarized in the right table of Table 1. Here, we clearly observe that for non smooth initial condition $u_{disc}^0$ a reduced space with higher dimension is needed for well approximating the solution (e.g. dimension 200 against 50 for reaching the machine precision). Concerning the MTD, we still obtain relative errors up to the machine precision with a one-dimensional time dependent reduced space ($E_2 = 3.9 \ 10^{-17}$ and $E_\infty = 1.3 \ 10^{-15}$).

**General case** We now consider the parameter-dependent problem and compare the approximations obtained with MTI and MTD for a subspace of dimension $r = 20$ generated from a training set of size 30. Figure 2 plots the evolution with time of the approximation obtained with MTD evaluated at $\xi = 0.65$. In Figure 3, the exact solution is compared to the approximations obtained by MTI and MTD, at final time and for $\xi = 0.65$. The

<footnote n="2">² Here, $\mathbf{1}_A$ is the characteristic function of a subset $A \subset \mathbb{R}^d$ and $\lfloor x \rfloor$ denotes the integer part of $x$.</footnote>
<table><caption>Table 1: Test case 1: relative errors $E_2$ and $E_\infty$ obtained with the MTI with respect to $dim(X_r)$ for both continuous (left) and discontinuous (right) initial conditions.</caption><thead><tr><th colspan="3">Relative errors for u<sup>0</sup><sub>cont</sub></th><th colspan="3">Relative errors for u<sup>0</sup><sub>disc</sub></th></tr><tr><th>dim(X<sub>r</sub>)</th><th>E<sub>2</sub></th><th>E<sub>∞</sub></th><th>dim(X<sub>r</sub>)</th><th>E<sub>2</sub></th><th>E<sub>∞</sub></th></tr></thead><tbody><tr><td>1</td><td>0.027831</td><td>0.76058</td><td>1</td><td>0.022116</td><td>0.76058</td></tr><tr><td>2</td><td>0.024075</td><td>0.95713</td><td>2</td><td>0.0135</td><td>0.46427</td></tr><tr><td>5</td><td>0.0078461</td><td>0.26984</td><td>5</td><td>0.0060424</td><td>0.2078</td></tr><tr><td>10</td><td>0.00021853</td><td>0.0075156</td><td>10</td><td>0.0040044</td><td>0.13771</td></tr><tr><td>20</td><td>4.8616e-09</td><td>1.6719e-07</td><td>20</td><td>0.002043</td><td>0.07026</td></tr><tr><td>50</td><td>4.0924e-17</td><td>1.4074e-15</td><td>50</td><td>4.7139e-05</td><td>0.0016211</td></tr><tr><td>100</td><td>4.4353e-17</td><td>1.5253e-15</td><td>100</td><td>1.3891e-12</td><td>4.7771e-11</td></tr><tr><td>200</td><td>4.5623e-17</td><td>1.569e-15</td><td>200</td><td>6.3966e-17</td><td>2.1998e-15</td></tr></tbody></table>

As shown on Figure 1, the approximations obtained with MTI are satisfactory when the dimension of the reduced space $X_r$ is greater than 20. This is confirmed by the values of the relative errors $E_q = ||u_r - u||_{I,q} / ||u||_{I,q}$ given in Table 1 for $q = 2$ and $q = \infty$, where $||u||_{I,q}$ is the natural norm on $L^q(I; X)$. In particular, the relative error in $L^2$-norm (resp $L^\infty$-norm) is of order $10^{-17}$ (respectively $10^{-15}$) for a space $X_r$ with a reduced dimension of 50. As expected, MTD, with a reduced space $X_1(t)$ of dimension 1 containing the exact solution, gives relative errors at the machine precision ($E_2 = 5.7 \ 10^{-17}$ and $E_\infty = 1.3 \ 10^{-15}$). For obtaining a very accurate precision, the method MTI requires a reduced space with rather high dimension, thus highlighting the limits of the POD method for transport equations. Now, we illustrate the impact of the smoothness of the solution. For this purpose, we consider the same test case with a discontinuous initial condition given by $u_{disc}^0(x) = \mathbf{1}_{[0.1,0.9]}(x) \cdot (\lfloor 3x \rfloor + \sin(10x))^2$. The approximate solution is plotted on Figure 1 and the relative errors between the reduced approximation and the exact solution are summarized in the right table of Table 1. Here, we clearly observe that for non smooth initial condition $u_{disc}^0$ a reduced space with higher dimension is needed for well approximating the solution (e.g. dimension 200 against 50 for reaching the machine precision). Concerning the MTD, we still obtain relative errors up to the machine precision with a one-dimensional time dependent reduced space ($E_2 = 3.9 \ 10^{-17}$ and $E_\infty = 1.3 \ 10^{-15}$).

**General case** We now consider the parameter-dependent problem and compare the approximations obtained with MTI and MTD for a subspace of dimension $r = 20$ generated from a training set of size 30. Figure 2 plots the evolution with time of the approximation obtained with MTD evaluated at $\xi = 0.65$. In Figure 3, the exact solution is compared to the approximations obtained by MTI and MTD, at final time and for $\xi = 0.65$. The

<footnote n="2">² Here, $\mathbf{1}_A$ is the characteristic function of a subset $A \subset \mathbb{R}^d$ and $\lfloor x \rfloor$ denotes the integer part of $x$.</footnote>
solutions of the reduced order models are very close to the exact solution. Nevertheless, we notice that the approximation computed with MTI present small oscillations whereas the one obtained with MTD matches very well the exact solution (see right plot of Figure 3).

<figure><img src="image_3.png" /><figcaption>Figure 2: Test case 1: Time evolution of $u_r$ computed with MTD for $r = 20$ and $\xi = 0.65$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Test case 1: Reduced approximations $u_r$ provided by MTD and MTI compared to the exact solution $u$ for both continuous (left) and discontinuous (right) initial conditions, at final time $T = 1$ for $r = 20$ and $\xi = 0.65$.</figcaption></figure>

Then, we estimate the expectation $\mathbb{E}(E_q(\xi))$ and maximum $\max_{\xi \in \Xi}(E_q(\xi))$ of the relative errors $E_q$ respectively by the empirical mean and by the maximum of the values of $E_q$ taken at 50 randomly chosen values of the parameters. These quantities are depicted on Figure 4 for different values of $r$ and for both continuous and discontinuous initial conditions. For the same dimension of the reduced spaces, MTD clearly provides a more accurate approximation than MTI in particular when considering discontinuous initial condition.
solutions of the reduced order models are very close to the exact solution. Nevertheless, we notice that the approximation computed with MTI present small oscillations whereas the one obtained with MTD matches very well the exact solution (see right plot of Figure 3).

<figure><img src="image_4.png" /><figcaption>Figure 2: Test case 1: Time evolution of $u_r$ computed with MTD for $r = 20$ and $\xi = 0.65$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Test case 1: Reduced approximations $u_r$ provided by MTD and MTI compared to the exact solution $u$ for both continuous (left) and discontinuous (right) initial conditions, at final time $T = 1$ for $r = 20$ and $\xi = 0.65$.</figcaption></figure>

Then, we estimate the expectation $\mathbb{E}(E_q(\xi))$ and maximum $\max_{\xi \in \Xi}(E_q(\xi))$ of the relative errors $E_q$ respectively by the empirical mean and by the maximum of the values of $E_q$ taken at 50 randomly chosen values of the parameters. These quantities are depicted on Figure 4 for different values of $r$ and for both continuous and discontinuous initial conditions. For the same dimension of the reduced spaces, MTD clearly provides a more accurate approximation than MTI in particular when considering discontinuous initial condition.
<figure><img src="image_5.png" /><figcaption>Figure 4: Test case 1: Statistical estimations of the expectation and maximum of the relative errors $E_2$ and $E_\infty$ with respect to the reduced dimension $r$, for MTD (left) and MTI (right) for both continuous (top) and discontinuous (bottom) initial conditions.</figcaption></figure>

**Validation of the error estimate** We now investigate the efficiency of the proposed error estimate $\tilde{\Delta}_r$. Figure 5 shows the evolution with time of $\tilde{\Delta}_r(t, \boldsymbol{\xi})$ and the exact error $\|e_r\|_X$ for MTD and MTI, evaluated at $\boldsymbol{\xi} = 0.65$. As expected, we observe that the error estimates have an exponential behavior in time. The error estimate $\tilde{\Delta}_r$ is much sharper for MTD than for MTI. For MTI, $\tilde{\Delta}_r$ is a very pessimistic upper bound of the error. This observation is confirmed by computing the effectivity index $\kappa(t, \boldsymbol{\xi}) = \tilde{\Delta}_r(t, \boldsymbol{\xi})/\|e_r(t, \boldsymbol{\xi})\|_X$. Figure 11 plots the evolution with time of a statistical estimation of the mean of the effectivity index $\mathbb{E}(\kappa(t, \boldsymbol{\xi}))$, which remains close to 1 with MTD and takes high values with MTI.
<figure><img src="image_6.png" /><figcaption>Figure 4: Test case 1: Statistical estimations of the expectation and maximum of the relative errors $E_2$ and $E_\infty$ with respect to the reduced dimension $r$, for MTD (left) and MTI (right) for both continuous (top) and discontinuous (bottom) initial conditions.</figcaption></figure>

**Validation of the error estimate** We now investigate the efficiency of the proposed error estimate $\tilde{\Delta}_r$. Figure 5 shows the evolution with time of $\tilde{\Delta}_r(t, \boldsymbol{\xi})$ and the exact error $\|e_r\|_X$ for MTD and MTI, evaluated at $\boldsymbol{\xi} = 0.65$. As expected, we observe that the error estimates have an exponential behavior in time. The error estimate $\tilde{\Delta}_r$ is much sharper for MTD than for MTI. For MTI, $\tilde{\Delta}_r$ is a very pessimistic upper bound of the error. This observation is confirmed by computing the effectivity index $\kappa(t, \boldsymbol{\xi}) = \tilde{\Delta}_r(t, \boldsymbol{\xi})/\|e_r(t, \boldsymbol{\xi})\|_X$. Figure 11 plots the evolution with time of a statistical estimation of the mean of the effectivity index $\mathbb{E}(\kappa(t, \boldsymbol{\xi}))$, which remains close to 1 with MTD and takes high values with MTI.
<figure><img src="image_6.png" /><figcaption>Figure 5: Test case 1: Evolution with time of $\tilde{\Delta}_r(t, \xi)$ and $\|e_r(t, \xi)\|_X$ for MTD (left) and MTI (right) for both continuous (top) and discontinuous (bottom) initial conditions, for $r = 20$ and $\xi = 0.65$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Test case 1: Evolution with time of a statistical estimation of $\mathbb{E}(\kappa(t, \xi))$ for MTD (left) and MTI (right) for both continuous and discontinuous initial conditions.</figcaption></figure>

**CPU times** Now, we briefly discuss the CPU computational costs of both MTI and MTD methods with POD-greedy and T-greedy algorithms respectively. We restrict the presentation to the case of the initial condition $u_{disc}^0$ (similar results are observed for $u_{cont}^0$). The dimension of reduced spaces is fixed to $r = 20$ for both methods. As summarized in Table 2, both offline and online costs are roughly of the same order but the MTD provides a reduced order model which is by 2 orders of magnitude (see Figure 4) better than the reduced order model provided by MTI. Yet, we notice that for the MTI, CPU-times are slightly higher. This is probably due to additional computations involved by time-dependent reduced quantities in MTD. We may guess that the CPU times for both methods would be similar for non autonomous problems with time affine dependent coefficients.
<figure><img src="image_7.png" /><figcaption>Figure 5: Test case 1: Evolution with time of $\tilde{\Delta}_r(t, \xi)$ and $\|e_r(t, \xi)\|_X$ for MTD (left) and MTI (right) for both continuous (top) and discontinuous (bottom) initial conditions, for $r = 20$ and $\xi = 0.65$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Test case 1: Evolution with time of a statistical estimation of $\mathbb{E}(\kappa(t, \xi))$ for MTD (left) and MTI (right) for both continuous and discontinuous initial conditions.</figcaption></figure>

**CPU times** Now, we briefly discuss the CPU computational costs of both MTI and MTD methods with POD-greedy and T-greedy algorithms respectively. We restrict the presentation to the case of the initial condition $u_{disc}^0$ (similar results are observed for $u_{cont}^0$). The dimension of reduced spaces is fixed to $r = 20$ for both methods. As summarized in Table 2, both offline and online costs are roughly of the same order but the MTD provides a reduced order model which is by 2 orders of magnitude (see Figure 4) better than the reduced order model provided by MTI. Yet, we notice that for the MTI, CPU-times are slightly higher. This is probably due to additional computations involved by time-dependent reduced quantities in MTD. We may guess that the CPU times for both methods would be similar for non autonomous problems with time affine dependent coefficients.
<table><caption>Table 2: Test case 1: Offline and online CPU-times for both MTD and MTI, with $u_{disc}^0$.</caption><thead><tr><th></th><th>Offline</th><th>Online</th></tr></thead><tbody><tr><th>MTD</th><td>1955.1</td><td>11.24</td></tr><tr><th>MTI</th><td>1517.5</td><td>8.94</td></tr></tbody></table>

## 5.2 Test case 2

We now consider a two-dimensional advection-diffusion equation with source term $g = 0$, $c(u, \xi) = 0$, $\mu(\xi) = \mu^0(2+\cos(\pi\xi^1)^2)$, with $\mu^0 = 0.5$, and $a(x, \xi) = a(\xi)(x_2-0.5, 0.5-x_1)^T$, with $a(\xi) = a_0 \sin(\pi\xi^2)$ and $a_0 = 0.1$. Here, we choose $\xi^1, \xi^2 \sim U(-1, 1)$ as independent uniform random variables. The spatial domain is $\Omega = (0, 1)^2$ and the time interval $I = (0, 0.2)$. The initial condition is given by $u^0(x) = e^{-(x_1-\frac{2}{3})^2-(x_2-\frac{2}{3})^2} \sin(2\pi x_1) \sin(2\pi x_2)$ and we impose homogeneous Dirichlet boundary conditions. We consider a $\mathbb{P}_1$-Lagrange finite element discretization with $n = 1681$ nodes$^3$ together with an implicit Euler scheme with $K = 400$ time steps, yielding the following scheme

$$ (\mathbf{I}_X - \delta t \mathbf{A}(\xi)) \mathbf{u}^{k+1}(\xi) = \mathbf{u}^k(\xi), \qquad (29) $$

with $\mathbf{A}(\xi) = \mu(\xi)\mathbf{A}_D + a(\xi)\mathbf{A}_C \in \mathbb{R}^{d \times d}$ ($d = 1521$), where $\mathbf{A}_C$ and $\mathbf{A}_D$ are the product of the mass matrix inverse with the discrete two dimensional diffusion and convection operators, obtained by finite element approximation, respectively.

We first present the reduced approximations computed by MTI and MTD with reduced spaces of dimension $r = 30$ selected with greedy algorithms using a training set of size 60. Figure 7 represents the approximation obtained with MTD evaluated for different values of the parameter $\xi \in \{(0.2, -1), (0.2, 0), (0.2, 0.5)\}$.

<figure><img src="image_8.png" /><figcaption>Figure 7: Test case 2: reduced approximation $\mathbf{u}_r$ computed with MTD for $r = 30$ at final time and $\xi \in \{(0.2, -1), (0.2, 0), (0.2, 0.5)\}$.</figcaption></figure>

For the same parameter values, the absolute distance to the exact solution at final time is given on Figure 7 for MTI and MTD. It shows that the approximations obtained by both

<footnote>$^3$The mesh is chosen fine enough to ensure a Péclet number smaller than 1.</footnote>
<table><caption>Table 2: Test case 1: Offline and online CPU-times for both MTD and MTI, with $u_{disc}^0$.</caption><thead><tr><th></th><th>Offline</th><th>Online</th></tr></thead><tbody><tr><th>MTD</th><td>1955.1</td><td>11.24</td></tr><tr><th>MTI</th><td>1517.5</td><td>8.94</td></tr></tbody></table>

## 5.2 Test case 2

We now consider a two-dimensional advection-diffusion equation with source term $g = 0$, $c(u, \xi) = 0$, $\mu(\xi) = \mu^0(2+\cos(\pi\xi^1)^2)$, with $\mu^0 = 0.5$, and $a(x, \xi) = a(\xi)(x_2-0.5, 0.5-x_1)^T$, with $a(\xi) = a_0 \sin(\pi\xi^2)$ and $a_0 = 0.1$. Here, we choose $\xi^1, \xi^2 \sim U(-1, 1)$ as independent uniform random variables. The spatial domain is $\Omega = (0, 1)^2$ and the time interval $I = (0, 0.2)$. The initial condition is given by $u^0(x) = e^{-(x_1-\frac{2}{3})^2-(x_2-\frac{2}{3})^2} \sin(2\pi x_1) \sin(2\pi x_2)$ and we impose homogeneous Dirichlet boundary conditions. We consider a $\mathbb{P}_1$-Lagrange finite element discretization with $n = 1681$ nodes$^3$ together with an implicit Euler scheme with $K = 400$ time steps, yielding the following scheme

$$ (\mathbf{I}_X - \delta t \mathbf{A}(\xi)) \mathbf{u}^{k+1}(\xi) = \mathbf{u}^k(\xi), \qquad (29) $$

with $\mathbf{A}(\xi) = \mu(\xi)\mathbf{A}_D + a(\xi)\mathbf{A}_C \in \mathbb{R}^{d \times d}$ ($d = 1521$), where $\mathbf{A}_C$ and $\mathbf{A}_D$ are the product of the mass matrix inverse with the discrete two dimensional diffusion and convection operators, obtained by finite element approximation, respectively.

We first present the reduced approximations computed by MTI and MTD with reduced spaces of dimension $r = 30$ selected with greedy algorithms using a training set of size 60. Figure 7 represents the approximation obtained with MTD evaluated for different values of the parameter $\xi \in \{(0.2, -1), (0.2, 0), (0.2, 0.5)\}$.

<figure><img src="image_9.png" /><figcaption>Figure 7: Test case 2: reduced approximation $\mathbf{u}_r$ computed with MTD for $r = 30$ at final time and $\xi \in \{(0.2, -1), (0.2, 0), (0.2, 0.5)\}$.</figcaption></figure>

For the same parameter values, the absolute distance to the exact solution at final time is given on Figure 7 for MTI and MTD. It shows that the approximations obtained by both

<footnote>$^3$The mesh is chosen fine enough to ensure a Péclet number smaller than 1.</footnote>
order reduction methods are in good agreement with the exact solution. Once again, we observe that MTD provides a better approximation with an error of order $10^{-9}$ against $10^{-5}$ for MTI.

<figure><img src="image_9.png" /><figcaption>Figure 8: Test case 2: absolute error between the exact solution and reduced approximations $|u - u_r|$ provided by MTD (top) and MTI (bottom) at final time for $r = 30$ and $\xi \in \{(0.2, -1), (0.2, 0), (0.2, 0.5)\}$.</figcaption></figure>

Estimations of the expectations $\mathbb{E}(E_q(\xi))$ and maxima $\max_{\xi \in \Xi}(E_q(\xi))$ of the relative errors $E_q$ (using a random sample of size 50 in $\Xi$) are plotted on Figure 9 for different values of $r$. For this advection-diffusion problem, we first observe that both order reduction methods provide accurate approximations with low-dimensional spaces. However, for the same dimension of the reduced spaces, the approximation obtained by MTD is more accurate than the approximation obtained by MTI. Indeed, in order to reach a relative error of $10^{-8}$, MTD requires $r = 10$ while MTI requires $r = 30$.
order reduction methods are in good agreement with the exact solution. Once again, we observe that MTD provides a better approximation with an error of order $10^{-9}$ against $10^{-5}$ for MTI.

<figure><img src="image_10.png" /><figcaption>Figure 8: Test case 2: absolute error between the exact solution and reduced approximations $|u - u_r|$ provided by MTD (top) and MTI (bottom) at final time for $r = 30$ and $\xi \in \{(0.2, -1), (0.2, 0), (0.2, 0.5)\}$.</figcaption></figure>

Estimations of the expectations $\mathbb{E}(E_q(\xi))$ and maxima $\max_{\xi \in \Xi}(E_q(\xi))$ of the relative errors $E_q$ (using a random sample of size 50 in $\Xi$) are plotted on Figure 9 for different values of $r$. For this advection-diffusion problem, we first observe that both order reduction methods provide accurate approximations with low-dimensional spaces. However, for the same dimension of the reduced spaces, the approximation obtained by MTD is more accurate than the approximation obtained by MTI. Indeed, in order to reach a relative error of $10^{-8}$, MTD requires $r = 10$ while MTI requires $r = 30$.
<figure><img src="image_10.png" /><figcaption>Figure 9: Test case 2: Statistical estimations of the expectation and maximum of the relative errors $E_2$ and $E_\infty$ with respect to the reduced dimension $r$, for MTD (left) and MTI (right).</figcaption></figure>

**Error estimates** Figure 10 represents the evolution with time of the exact error $\|e_r(t, \xi)\|_X$ and of the error estimate $\tilde{\Delta}_r(t, \xi)$ for both order reduction methods, for $\xi = (0.2, 0.5)$. Again, we observe the exponential behavior of the error estimates. The errors obtained by MTD are several orders of magnitude lower than the errors obtained by MTI.

<figure><img src="image_11.png" /><figcaption>Figure 10: Test case 2: Evolution with time of the error estimate $\tilde{\Delta}_r$ and true error $\|e_r\|_X$ for MTD (left) and MTI (right), for $r = 30$ and $\xi = (0.2, 0.5)$.</figcaption></figure>

For the first time indices (i.e. $k < 100$), $\tilde{\Delta}_r$ provides a sharper error bound for MTD than for MTI. Also, $\tilde{\Delta}_r$ and the true error have quite similar time evolutions for MTD, and very different time evolutions for MTI. The superiority of MTD over MTI is confirmed on Figure 11, where we observe that the expectations of the effectivity index $\kappa(t, \xi)$ of MTI and MTD are in a ratio of 10. For longer times, the error estimates present similar efficiency for both approaches.
<figure><img src="image_11.png" /><figcaption>Figure 9: Test case 2: Statistical estimations of the expectation and maximum of the relative errors $E_2$ and $E_\infty$ with respect to the reduced dimension $r$, for MTD (left) and MTI (right).</figcaption></figure>

**Error estimates** Figure 10 represents the evolution with time of the exact error $\|e_r(t, \xi)\|_X$ and of the error estimate $\tilde{\Delta}_r(t, \xi)$ for both order reduction methods, for $\xi = (0.2, 0.5)$. Again, we observe the exponential behavior of the error estimates. The errors obtained by MTD are several orders of magnitude lower than the errors obtained by MTI.

<figure><img src="image_12.png" /><figcaption>Figure 10: Test case 2: Evolution with time of the error estimate $\tilde{\Delta}_r$ and true error $\|e_r\|_X$ for MTD (left) and MTI (right), for $r = 30$ and $\xi = (0.2, 0.5)$.</figcaption></figure>

For the first time indices (i.e. $k < 100$), $\tilde{\Delta}_r$ provides a sharper error bound for MTD than for MTI. Also, $\tilde{\Delta}_r$ and the true error have quite similar time evolutions for MTD, and very different time evolutions for MTI. The superiority of MTD over MTI is confirmed on Figure 11, where we observe that the expectations of the effectivity index $\kappa(t, \xi)$ of MTI and MTD are in a ratio of 10. For longer times, the error estimates present similar efficiency for both approaches.
<figure><img src="image_12.png" /><figcaption>Figure 11: Test case 2: Evolution with time of a statistical estimation of the expected effectivity index $E(\kappa(t, \xi))$ for MTD (left) and MTI (right).</figcaption></figure>

**Study of greedy algorithms** Now, we study the behavior of the T-greedy and POD-greedy algorithms used for the construction of reduced spaces for MTD and MTI respectively. For both algorithms, the evolution with the iteration $r$ of the maximum of $\Delta_r^{(0,T)}(\xi)$ over the training set $\Xi_{train}$ is plotted on Figure 12. Note that at iteration $r$, the dimension of the reduced space is actually $r+1$. We observe that the error decreases faster with $r$ for the T-greedy algorithm than for the POD-greedy algorithm. For example, at iteration 20 of the greedy algorithms, the errors for POD-greedy and T-greedy algorithms are respectively of order $10^{-1}$ and 10. When stopping the greedy algorithms at a given precision, the dimension of the obtained reduced space is much smaller for T-greedy than for POD-greedy. We observe on Table 13 that the POD-greedy algorithm can select several times the same point in $\Xi$ (for example, the value $\xi = (-0.3192, 0.9004)$ is selected 6 times during the 30 iterations). This is due to the fact that at each iteration, only the first term of the POD of the selected trajectory is added to the reduced space (here $\ell = 1$). Concerning the T-greedy approach, a point in $\Xi$ can be selected only once. This implies that for a fixed dimension, the POD-greedy approach requires less solutions of the full order dynamical system (here 5 for POD-greedy against 30 for the T-greedy approach when considering a reduced space of dimension 30). To illustrate this point, the number of calls to the full order problem with respect to $\max_{\xi \in \Xi_{train}} \Delta_r^{(0,T)}(\xi)$ is represented on the right curve of Figure 12. The number of calls to the full order problem for POD-greedy algorithm remains smaller, which means that the POD-greedy algorithm uses only the information of snapshots already computed to enrich $X_r$ whereas the T-greedy algorithm requires new snapshots at each iteration. In terms of computational costs, the POD-greedy strategy presents lower offline computational costs since it requires a smaller number of calls to the full order problem, but for a given precision, it leads to larger reduced systems in the online step. For the T-greedy algorithm, we have the opposite conclusions. For a given precision, it leads to lower dimensional reduced spaces, and therefore to lower online costs, but it requires more calls to the full-order problem in the offline step. The MTD approach with POD-greedy remains
<figure><img src="image_13.png" /><figcaption>Figure 11: Test case 2: Evolution with time of a statistical estimation of the expected effectivity index $E(\kappa(t, \boldsymbol{\xi}))$ for MTD (left) and MTI (right).</figcaption></figure>

**Study of greedy algorithms** Now, we study the behavior of the T-greedy and POD-greedy algorithms used for the construction of reduced spaces for MTD and MTI respectively. For both algorithms, the evolution with the iteration $r$ of the maximum of $\Delta_r^{(0,T)}(\boldsymbol{\xi})$ over the training set $\Xi_{train}$ is plotted on Figure 12. Note that at iteration $r$, the dimension of the reduced space is actually $r+1$. We observe that the error decreases faster with $r$ for the T-greedy algorithm than for the POD-greedy algorithm. For example, at iteration 20 of the greedy algorithms, the errors for POD-greedy and T-greedy algorithms are respectively of order $10^{-1}$ and 10. When stopping the greedy algorithms at a given precision, the dimension of the obtained reduced space is much smaller for T-greedy than for POD-greedy. We observe on Table 13 that the POD-greedy algorithm can select several times the same point in $\Xi$ (for example, the value $\boldsymbol{\xi} = (-0.3192, 0.9004)$ is selected 6 times during the 30 iterations). This is due to the fact that at each iteration, only the first term of the POD of the selected trajectory is added to the reduced space (here $\ell = 1$). Concerning the T-greedy approach, a point in $\Xi$ can be selected only once. This implies that for a fixed dimension, the POD-greedy approach requires less solutions of the full order dynamical system (here 5 for POD-greedy against 30 for the T-greedy approach when considering a reduced space of dimension 30). To illustrate this point, the number of calls to the full order problem with respect to $\max_{\boldsymbol{\xi} \in \Xi_{train}} \Delta_r^{(0,T)}(\boldsymbol{\xi})$ is represented on the right curve of Figure 12. The number of calls to the full order problem for POD-greedy algorithm remains smaller, which means that the POD-greedy algorithm uses only the information of snapshots already computed to enrich $X_r$ whereas the T-greedy algorithm requires new snapshots at each iteration. In terms of computational costs, the POD-greedy strategy presents lower offline computational costs since it requires a smaller number of calls to the full order problem, but for a given precision, it leads to larger reduced systems in the online step. For the T-greedy algorithm, we have the opposite conclusions. For a given precision, it leads to lower dimensional reduced spaces, and therefore to lower online costs, but it requires more calls to the full-order problem in the offline step. The MTD approach with POD-greedy remains
a good alternative when treating autonomous linear dynamical systems arising from the discretization of linear parabolic problems with time independent matrices. Furthermore, contrarily to MTI, it does not require to store sequences of time dependent reduced quantities. But when considering non-autonomous dynamical systems (which requires handling time-dependent reduced quantities), the storing costs will be smaller with MTD, since this latter approach provides smaller reduced spaces for a given precision.

<figure><img src="image_13.png" /><figcaption>Figure 12: Test case 2 : Evolution of the maximum of $\tilde{\Delta}_r^{(0,T)}(\xi)$ over $\Xi_{train}$ with the iteration $r$ (left), and number of calls to the full order problem with respect to the maximum of $\tilde{\Delta}_r^{(0,T)}(\xi)$ over $\Xi_{train}$ (right), for T-greedy (blue line) and POD-greedy (red line) algorithms.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13: Test case 2: Selected points in $\Xi$ for the T-greedy (left) and POD-greedy (right) algorithms during 30 iterations. For POD-greedy, the size of a point increases with the number of times it is selected.</figcaption></figure>

**CPU times** As shown in Table 3, CPU-times for both offline and online phases are again similar, when the reduced space dimension is fixed to $r = 30$ for both MTI and MTD. We still observe a small additional cost for the online phase of MTD due to computations with
a good alternative when treating autonomous linear dynamical systems arising from the discretization of linear parabolic problems with time independent matrices. Furthermore, contrarily to MTI, it does not require to store sequences of time dependent reduced quantities. But when considering non-autonomous dynamical systems (which requires handling time-dependent reduced quantities), the storing costs will be smaller with MTD, since this latter approach provides smaller reduced spaces for a given precision.

<figure><img src="image_14.png" /><figcaption>Figure 12: Test case 2 : Evolution of the maximum of $\tilde{\Delta}_r^{(0,T)}(\xi)$ over $\Xi_{train}$ with the iteration $r$ (left), and number of calls to the full order problem with respect to the maximum of $\tilde{\Delta}_r^{(0,T)}(\xi)$ over $\Xi_{train}$ (right), for T-greedy (blue line) and POD-greedy (red line) algorithms.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13: Test case 2: Selected points in $\Xi$ for the T-greedy (left) and POD-greedy (right) algorithms during 30 iterations. For POD-greedy, the size of a point increases with the number of times it is selected.</figcaption></figure>

**CPU times** As shown in Table 3, CPU-times for both offline and online phases are again similar, when the reduced space dimension is fixed to $r = 30$ for both MTI and MTD. We still observe a small additional cost for the online phase of MTD due to computations with
time-dependent reduced quantities, but MTD again provides a better reduced approximation (see Figure 9). Concerning the POD-greedy algorithm, one could imagine to improve the offline cost by storing the solution trajectories computed during the greedy procedure for parameters selected several times, but this would imply additional storage costs.

<table><caption>Table 3: Test case 2: Offline and online CPU-times for both MTD and MTI.</caption><thead><tr><th></th><th>Offline - Basis selection</th><th>Offline - Reduced quantities</th><th>Online</th></tr></thead><tbody><tr><td>T-Greedy</td><td>34907</td><td>77.69</td><td>37.86</td></tr><tr><td>POD-Greedy</td><td>32767</td><td>0.22</td><td>34.6</td></tr></tbody></table>

## 5.3 Test case 3

We study a nonlinear viscous Burger's equation with uncertain parameters. This test case is adapted from [34]. We consider a spatial domain $\Omega = (0, 1)$, a time interval $I = (0, 1)$, homogeneous Dirichlet boundary conditions and an initial condition $u^0 = 0$. We consider $c(u, \xi) = \frac{1}{2}u$, $a(x, \xi) = 0$ and a diffusion coefficient $\mu(\xi) = \xi$, with $\xi \sim U(0.01, 0.06)$ a uniform random variable. The source term is defined by $g(x, t, \xi) = g_1(x, t) + g_2(x, t)$, with

$$g_1(x, t) = 4e^{-(\frac{x-0.2}{0.03})^2} \mathbf{1}_{[0.1,0.3]}(x) \sin(4\pi t), \\ g_2(x, t) = 4 \cdot \mathbf{1}_{[0.6,0.7]}(x) \mathbf{1}_{[0.2,0.4]}(t).$$

The term $g_1$ corresponds to an excitation localized in space and oscillating with time, whereas the term $g_2$ corresponds to a constant excitation over a localized space-time region. We use a finite difference scheme in space (with $n = 300$ nodes, and $d = 298$) with uniform mesh (as in Subsection 5.1) and a semi-implicit Euler scheme in time (with $K = 200$ time steps), yielding the following scheme:

$$(\mathbf{I}_X - \delta t \mathbf{A}(\xi))\mathbf{u}^{k+1}(\xi) = \mathbf{u}^k(\xi) + \delta t(\mathbf{h}(\mathbf{u}^k) + \mathbf{g}^k), \quad (30)$$

where $\mathbf{A} \in \mathbb{R}^{d \times d}$ denotes the discrete second derivative operator, and where the discrete flux $\mathbf{h}(\mathbf{u}^k) \in \mathbb{R}^d$ is defined by $(\mathbf{h}(\mathbf{u}^k))_i = -u_i^k(\mathbf{C}\mathbf{u}^k)_i$ with $\mathbf{C} \in \mathbb{R}^{d \times d}$ the matrix associated with the discrete first derivative operator.

Reduced spaces are constructed with POD-greedy and T-greedy algorithms with a training set of size 60. For the error estimation, we have used a nearest neighbor interpolation of the Lipschitz constant of the flux $\mathbf{h}$. The maximum $M_L(\xi)$ and minimum $m_L(\xi)$ of $\frac{\tilde{L}_X[\mathbf{h}](\mathbf{u}_r(t,\xi),t,\xi)}{L_X[\mathbf{h}](\mathbf{u}_r(t,\xi),t,\xi)}$ over the time interval $I$ have been computed for both MTD and MTI methods (see Table 4). We observe that the values are localized in intervals whose bounds are close to 1, which indicates that the approximations are quite satisfactory. For the efficient evaluation of the nonlinear flux $\mathbf{h}$, an EIM has been used with
time-dependent reduced quantities, but MTD again provides a better reduced approximation (see Figure 9). Concerning the POD-greedy algorithm, one could imagine to improve the offline cost by storing the solution trajectories computed during the greedy procedure for parameters selected several times, but this would imply additional storage costs.

<table><caption>Table 3: Test case 2: Offline and online CPU-times for both MTD and MTI.</caption><thead><tr><th></th><th>Offline - Basis selection</th><th>Offline - Reduced quantities</th><th>Online</th></tr></thead><tbody><tr><th>T-Greedy</th><td>34907</td><td>77.69</td><td>37.86</td></tr><tr><th>POD-Greedy</th><td>32767</td><td>0.22</td><td>34.6</td></tr></tbody></table>

## 5.3 Test case 3

We study a nonlinear viscous Burger's equation with uncertain parameters. This test case is adapted from [34]. We consider a spatial domain $\Omega = (0, 1)$, a time interval $I = (0, 1)$, homogeneous Dirichlet boundary conditions and an initial condition $u^0 = 0$. We consider $c(u, \xi) = \frac{1}{2}u$, $a(x, \xi) = 0$ and a diffusion coefficient $\mu(\xi) = \xi$, with $\xi \sim U(0.01, 0.06)$ a uniform random variable. The source term is defined by $g(x, t, \xi) = g_1(x, t) + g_2(x, t)$, with

$$g_1(x, t) = 4e^{-(\frac{x-0.2}{0.03})^2} \mathbf{1}_{[0.1,0.3]}(x) \sin(4\pi t),$$

$$g_2(x, t) = 4 \cdot \mathbf{1}_{[0.6,0.7]}(x) \mathbf{1}_{[0.2,0.4]}(t).$$

The term $g_1$ corresponds to an excitation localized in space and oscillating with time, whereas the term $g_2$ corresponds to a constant excitation over a localized space-time region. We use a finite difference scheme in space (with $n = 300$ nodes, and $d = 298$) with uniform mesh (as in Subsection 5.1) and a semi-implicit Euler scheme in time (with $K = 200$ time steps), yielding the following scheme:

$$(\mathbf{I}_X - \delta t \mathbf{A}(\xi))\mathbf{u}^{k+1}(\xi) = \mathbf{u}^k(\xi) + \delta t(\mathbf{h}(\mathbf{u}^k) + \mathbf{g}^k), \quad (30)$$

where $\mathbf{A} \in \mathbb{R}^{d \times d}$ denotes the discrete second derivative operator, and where the discrete flux $\mathbf{h}(\mathbf{u}^k) \in \mathbb{R}^d$ is defined by $(\mathbf{h}(\mathbf{u}^k))_i = -u_i^k(\mathbf{C}\mathbf{u}^k)_i$ with $\mathbf{C} \in \mathbb{R}^{d \times d}$ the matrix associated with the discrete first derivative operator.

Reduced spaces are constructed with POD-greedy and T-greedy algorithms with a training set of size 60. For the error estimation, we have used a nearest neighbor interpolation of the Lipschitz constant of the flux $\mathbf{h}$. The maximum $M_L(\xi)$ and minimum $m_L(\xi)$ of $\frac{\tilde{L}_X[\mathbf{h}](\mathbf{u}_r(t, \xi), t, \xi)}{L_X[\mathbf{h}](\mathbf{u}_r(t, \xi), t, \xi)}$ over the time interval $I$ have been computed for both MTD and MTI methods (see Table 4). We observe that the values are localized in intervals whose bounds are close to 1, which indicates that the approximations are quite satisfactory. For the efficient evaluation of the nonlinear flux $\mathbf{h}$, an EIM has been used with
a tolerance $\epsilon = 10^{-10}$ in the stopping criterion (24). For the considered simulations, it corresponds to an average number $m$ of terms in the EIM equal to 103 and 108 for MTD and MTI respectively. As shown in Table 4, the maximum of the relative errors $E_h(t, \xi) = \|h(u_r(t, \xi), t, \xi) - \tilde{h}(u_r(t, \xi), t, \xi)\|_X / \|h(u_r(t, \xi), t, \xi)\|_X$ on the flux is of order $10^{-10}$ as expected.

<table><caption>Table 4: Test case 3: Maximum over the time interval of the relative approximation error of the flux $h$ with EIM, and minimum and maximum values $m_L$ and $M_L$ of the ratio between approximate and exact Lipschitz constants of $h$, for MTD and MTI.</caption><thead><tr><th rowspan="2">ξ</th><th colspan="2">max<sub>t</sub> E<sub>h</sub>(ξ, t)</th><th colspan="2">[m<sub>L</sub>(ξ), M<sub>L</sub>(ξ)]</th></tr><tr><th>MTI</th><th>MTD</th><th>MTI</th><th>MTD</th></tr></thead><tbody><tr><td>0.01</td><td>1.5550 10<sup>-10</sup></td><td>4.2208 10<sup>-10</sup></td><td>[0.9116, 1.1011]</td><td>[1.1055, 1.1297]</td></tr><tr><td>0.035</td><td>3.5803 10<sup>-10</sup></td><td>1.4400 10<sup>-11</sup></td><td>[0.7804, 0.9536]</td><td>[1.0062, 1.0483]</td></tr><tr><td>0.06</td><td>5.0374 10<sup>-10</sup></td><td>1.4201 10<sup>-11</sup></td><td>[0.7567, 1.0009]</td><td>[0.9947, 0.9989]</td></tr></tbody></table>

We first study the behavior of the approximations provided by the MTD and MTI evaluated for the parameter values $\xi \in \{0.01, 0.035, 0.06\}$ and for $r = 15$. The evolution of the solution of the reduced order model obtained with MTD is plotted on Figure 14. We clearly observe different features for different values of the diffusion coefficient. Figure 15 represents the exact solution and the approximations obtained by MTD and MTI at final time and for a given value of the parameter. For $r = 15$, the reduced approximations are in good agreement with the exact solution. Nevertheless, we note that the approximation obtained with MTD is clearly better than with MTI.

<figure><img src="image_15.png" /><figcaption>Figure 14: Test case 3: Evolution over the space-time grid of the approximation $u_r$ computed with MTD for $r = 15$ and for $\xi \in \{0.01, 0.035, 0.06\}$.</figcaption></figure>
a tolerance $\epsilon = 10^{-10}$ in the stopping criterion (24). For the considered simulations, it corresponds to an average number $m$ of terms in the EIM equal to 103 and 108 for MTD and MTI respectively. As shown in Table 4, the maximum of the relative errors $E_h(t, \xi) = \|h(u_r(t, \xi), t, \xi) - \tilde{h}(u_r(t, \xi), t, \xi)\|_X / \|h(u_r(t, \xi), t, \xi)\|_X$ on the flux is of order $10^{-10}$ as expected.

<table><caption>Table 4: Test case 3: Maximum over the time interval of the relative approximation error of the flux $h$ with EIM, and minimum and maximum values $m_L$ and $M_L$ of the ratio between approximate and exact Lipschitz constants of $h$, for MTD and MTI.</caption><thead><tr><th rowspan="2">ξ</th><th colspan="2">max<sub>t</sub> E<sub>h</sub>(ξ, t)</th><th colspan="2">[m<sub>L</sub>(ξ), M<sub>L</sub>(ξ)]</th></tr><tr><th>MTI</th><th>MTD</th><th>MTI</th><th>MTD</th></tr></thead><tbody><tr><td>0.01</td><td>1.5550 10<sup>-10</sup></td><td>4.2208 10<sup>-10</sup></td><td>[0.9116, 1.1011]</td><td>[1.1055, 1.1297]</td></tr><tr><td>0.035</td><td>3.5803 10<sup>-10</sup></td><td>1.4400 10<sup>-11</sup></td><td>[0.7804, 0.9536]</td><td>[1.0062, 1.0483]</td></tr><tr><td>0.06</td><td>5.0374 10<sup>-10</sup></td><td>1.4201 10<sup>-11</sup></td><td>[0.7567, 1.0009]</td><td>[0.9947, 0.9989]</td></tr></tbody></table>

We first study the behavior of the approximations provided by the MTD and MTI evaluated for the parameter values $\xi \in \{0.01, 0.035, 0.06\}$ and for $r = 15$. The evolution of the solution of the reduced order model obtained with MTD is plotted on Figure 14. We clearly observe different features for different values of the diffusion coefficient. Figure 15 represents the exact solution and the approximations obtained by MTD and MTI at final time and for a given value of the parameter. For $r = 15$, the reduced approximations are in good agreement with the exact solution. Nevertheless, we note that the approximation obtained with MTD is clearly better than with MTI.

<figure><img src="image_16.png" /><figcaption>Figure 14: Test case 3: Evolution over the space-time grid of the approximation $u_r$ computed with MTD for $r = 15$ and for $\xi \in \{0.01, 0.035, 0.06\}$.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15: Test case 3: Comparison between the exact solution $u$ (black line) and the approximations $u_r$ obtained with MTD (blue dashed line) and MTI (red dotted line) with $r = 15$, at final time $T = 0.1$ and for $\xi \in \{0.01, 0.035, 0.06\}$.</figcaption></figure>

This superiority of MTD over MTI is confirmed by Figure 16 where we have plotted statistical estimations (with a random sample of size 50) of the expectation and maximum of the relative errors. As we can see, MTD provides with $r = 15$ an approximation with a relative error of $10^{-10}$, whereas MTI only provides an approximation with a relative error of $10^{-2}$ for the same dimension of the reduced space. For $r = 50$, MTI provides an approximation with relative error $10^{-6}$, which is still higher than MTD with $r = 15$.

<figure><img src="image_17.png" /><figcaption>Figure 16: Test case 3: Statistical estimations of the expectations and maxima of the relative errors $E_2$ and $E_\infty$ with respect to the reduced dimension $r$, for MTD (left) and MTI (right).</figcaption></figure>

**Study of error estimates** We now compare the behavior of the error estimate $\tilde{\Delta}_r$ for MTD and MTI. Figure 17 represent the evolution with time of the exact error and of the error estimate for both methods, with and without approximation of the flux and the Lipschitz constants, for $\xi = 0.035$ and $r = 15$. We first remark that the EIM approximation of the nonlinear flux as well as the approximated Lipschitz constants have a small impact on the behavior of the error estimate. Nevertheless we notice that for the first instants for MTD, where the error is smaller than $10^{-10}$, the EIM error is no longer negligible and
<figure><img src="image_17.png" /><figcaption>Figure 15: Test case 3: Comparison between the exact solution $u$ (black line) and the approximations $u_r$ obtained with MTD (blue dashed line) and MTI (red dotted line) with $r = 15$, at final time $T = 0.1$ and for $\xi \in \{0.01, 0.035, 0.06\}$.</figcaption></figure>

This superiority of MTD over MTI is confirmed by Figure 16 where we have plotted statistical estimations (with a random sample of size 50) of the expectation and maximum of the relative errors. As we can see, MTD provides with $r = 15$ an approximation with a relative error of $10^{-10}$, whereas MTI only provides an approximation with a relative error of $10^{-2}$ for the same dimension of the reduced space. For $r = 50$, MTI provides an approximation with relative error $10^{-6}$, which is still higher than MTD with $r = 15$.

<figure><img src="image_18.png" /><figcaption>Figure 16: Test case 3: Statistical estimations of the expectations and maxima of the relative errors $E_2$ and $E_\infty$ with respect to the reduced dimension $r$, for MTD (left) and MTI (right).</figcaption></figure>

**Study of error estimates** We now compare the behavior of the error estimate $\tilde{\Delta}_r$ for MTD and MTI. Figure 17 represent the evolution with time of the exact error and of the error estimate for both methods, with and without approximation of the flux and the Lipschitz constants, for $\xi = 0.035$ and $r = 15$. We first remark that the EIM approximation of the nonlinear flux as well as the approximated Lipschitz constants have a small impact on the behavior of the error estimate. Nevertheless we notice that for the first instants for MTD, where the error is smaller than $10^{-10}$, the EIM error is no longer negligible and
$\|e_r\|_X$ and $\tilde{\Delta}_r$ are slightly impacted. For MTD, a small difference is observed for final times only on $\tilde{\Delta}_r$. This is probably due to the approximation of the Lipchitz constant. These differences disappear when choosing a smaller precision $\varepsilon$ for MTD or when considering larger reduced spaces for MTI. Again, we obtain several orders of magnitude between the errors provided by the MTD and MTI. Here, the expectations of the effectivity index $\kappa(t, \xi)$ (plotted on Figure 18) are of the same order, with or without EIM, which means that the error estimate has the same efficiency for both approaches.

<figure><img src="image_18.png" /><figcaption>Figure 17: Test case 3: Evolution with time of the error estimate $\tilde{\Delta}_r$ (dotted line) and the exact error $\|e_r\|_X$ (solid line), for MTD (left) and MTI (right), for $\xi = 0.035$ and $r = 15$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 18: Test case 3: Evolution with time of a statistical estimation of the expected effectivity index $\mathbb{E}(\kappa(t, \xi))$, for MTD (left) and MTI (right), for $\xi = 0.035$ and $r = 15$.</figcaption></figure>

**Study of greedy algorithms** We finally study the behavior of the greedy algorithms used for the construction of the reduced space $X_r$. Once again, we have plotted the evolution of the maximum of $\tilde{\Delta}_r^{(0,T)}(\xi)$ over $\Xi_{train}$ with the iteration $r$ and the number of calls to the full order problem with respect to the maximum of $\tilde{\Delta}_r^{(0,T)}(\xi)$ on Figure 19.
$\|e_r\|_X$ and $\tilde{\Delta}_r$ are slightly impacted. For MTD, a small difference is observed for final times only on $\tilde{\Delta}_r$. This is probably due to the approximation of the Lipchitz constant. These differences disappear when choosing a smaller precision $\varepsilon$ for MTD or when considering larger reduced spaces for MTI. Again, we obtain several orders of magnitude between the errors provided by the MTD and MTI. Here, the expectations of the effectivity index $\kappa(t, \xi)$ (plotted on Figure 18) are of the same order, with or without EIM, which means that the error estimate has the same efficiency for both approaches.

<figure><img src="image_19.png" /><figcaption>Figure 17: Test case 3: Evolution with time of the error estimate $\tilde{\Delta}_r$ (dotted line) and the exact error $\|e_r\|_X$ (solid line), for MTD (left) and MTI (right), for $\xi = 0.035$ and $r = 15$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 18: Test case 3: Evolution with time of a statistical estimation of the expected effectivity index $\mathbb{E}(\kappa(t, \xi))$, for MTD (left) and MTI (right), for $\xi = 0.035$ and $r = 15$.</figcaption></figure>

**Study of greedy algorithms** We finally study the behavior of the greedy algorithms used for the construction of the reduced space $X_r$. Once again, we have plotted the evolution of the maximum of $\tilde{\Delta}_r^{(0,T)}(\xi)$ over $\Xi_{train}$ with the iteration $r$ and the number of calls to the full order problem with respect to the maximum of $\tilde{\Delta}_r^{(0,T)}(\xi)$ on Figure 19.
We observe the same behavior as for the test case 2. The maximum error during the greedy iterations decreases faster with T-greedy than with POD-greedy (T-greedy reaches a precision of order $10^{-6}$ for the iteration 14, while POD-greedy reaches a precision of $10^{-2}$ for the iteration 49). In terms of computational costs, the number of calls to the full order problem in the offline step is lower for T-greedy algorithm than for POD-greedy (e.g., for reaching a precision of order $10^{-2}$, T-greedy (resp. POD-greedy) algorithm requires 7 (resp. 11) calls to the full order problem). Contrary to the previous linear test case, the use of MTD with T-greedy algorithm clearly reduces the computational costs of both online and offline steps in comparison to a MTI combined with a POD-greedy construction of reduced spaces.

<figure><img src="image_20.png" /><figcaption>Figure 19: Test case 3: Evolution of the maximum of $\tilde{\Delta}_r^{(0,T)}(\xi)$ over $\Xi_{train}$ with the iteration $r$ (left), and number of calls to the full order problem with respect to the maximum of $\tilde{\Delta}_r^{(0,T)}(\xi)$ over $\Xi_{train}$ (right), for T-greedy (blue line) and POD-greedy (red line) algorithms.</figcaption></figure>

**CPU times** The dimension of reduced spaces is fixed to $r = 15$ for both methods. As summarized in the table Table 5, once again for this problem, both offline and online cost are roughly of the same order for both MTI and MTD while MTD gives a finer reduced approximation which is by 8 orders of magnitude (see Figure 16) better than the one provided by MTI.

<table><caption>Table 5: Test case 1: Offline and online CPU-times for both MTD and MTI.</caption><thead><tr><th></th><th>Offline (Basis)</th><th>Offline (Reduction)</th></tr></thead><tbody><tr><td>MTI</td><td>6137.2</td><td>9.22</td></tr><tr><td>MTD</td><td>6253.8</td><td>9.27</td></tr></tbody></table>
We observe the same behavior as for the test case 2. The maximum error during the greedy iterations decreases faster with T-greedy than with POD-greedy (T-greedy reaches a precision of order $10^{-6}$ for the iteration 14, while POD-greedy reaches a precision of $10^{-2}$ for the iteration 49). In terms of computational costs, the number of calls to the full order problem in the offline step is lower for T-greedy algorithm than for POD-greedy (e.g., for reaching a precision of order $10^{-2}$, T-greedy (resp. POD-greedy) algorithm requires 7 (resp. 11) calls to the full order problem). Contrary to the previous linear test case, the use of MTD with T-greedy algorithm clearly reduces the computational costs of both online and offline steps in comparison to a MTI combined with a POD-greedy construction of reduced spaces.

<figure><img src="image_21.png" /><figcaption>Figure 19: Test case 3: Evolution of the maximum of $\tilde{\Delta}_r^{(0,T)}(\xi)$ over $\Xi_{train}$ with the iteration $r$ (left), and number of calls to the full order problem with respect to the maximum of $\tilde{\Delta}_r^{(0,T)}(\xi)$ over $\Xi_{train}$ (right), for T-greedy (blue line) and POD-greedy (red line) algorithms.</figcaption></figure>

**CPU times** The dimension of reduced spaces is fixed to $r = 15$ for both methods. As summarized in the table Table 5, once again for this problem, both offline and online cost are roughly of the same order for both MTI and MTD while MTD gives a finer reduced approximation which is by 8 orders of magnitude (see Figure 16) better than the one provided by MTI.

<table><caption>Table 5: Test case 1: Offline and online CPU-times for both MTD and MTI.</caption><thead><tr><th></th><th>Offline (Basis)</th><th>Offline (Reduction)</th></tr></thead><tbody><tr><th>MTI</th><td>6137.2</td><td>9.22</td></tr><tr><th>MTD</th><td>6253.8</td><td>9.27</td></tr></tbody></table>
# 6 Conclusion

In this paper, we have presented a projection-based model order reduction approach for the solution of parameter-dependent dynamical systems. It relies on a Galerkin projection of the full order dynamical system on time-dependent reduced spaces. This generalizes classical model order reduction methods with time-independent reduced spaces. An a posteriori error estimate using the logarithmic Lipschitz constant associated to the flux has also been proposed and provides an efficient a posteriori error estimate. Using this error estimate, we have derived a greedy algorithm (called T-greedy algorithm) for the construction of a sequence of time-dependent reduced spaces. We have performed several numerical tests on both linear and nonlinear dynamical systems. The order reduction method with time-dependent spaces constructed with a T-greedy algorithm (MTD) has been compared to a method with time-independent spaces constructed with a POD-greedy algorithm (MTI). For the same online costs (time for evaluating the ROM), the error in the predictions of the ROM obtained by MTD may be several orders of magnitude lower than the error in the predictions of the ROM obtained by MTI. In other words, for a given precision of the ROM, the MTD approach yields a significant reduction of online computation costs. Here, the proposed T-greedy algorithm for the selection of parameters relies on an error indicator using a $L^2$-norm in time of an error estimate. This global error estimate could be improved by considering sharper local in time error estimates which could be obtained by improving the local Lipschitz constants used to derive our error bound [34]. Moreover, a goal-oriented variant of this algorithm could be introduced by considering norms taking into account a certain quantity of interest (e.g., the value of the solution at the final time) in order to improve the efficiency of the order reduction method for this approximation of this quantity of interest.

# References

[1] M. BARRAULT, Y. MADAY, N.C. NGUYEN AND A.T. PATERA, *An "empirical interpolation" method: application to efficient reduced-basis discretization of partial differential equations*. Comptes Rendus Mathématique, 339(9):667-672, 2004.

[2] U. BAUR, P. BENNER, B. HAASDONK, C. HIMPE, I. MAIERN AND M. OHLBERGER, *Comparison of methods for parametric model order reduction of instationary problems*, preprint, 2015.

[3] M. BEBENDORF, Y. MADAY AND B. STAMM, *Comparison of some Reduced Representation Approximations*, Quarteroni, Alfio and Rozza, Gianluigi, Reduced Order Methods for Modeling and Computational Reduction, MS&A - Modeling, Simulation and Applications 9:67-100, 2014.
# 6 Conclusion

In this paper, we have presented a projection-based model order reduction approach for the solution of parameter-dependent dynamical systems. It relies on a Galerkin projection of the full order dynamical system on time-dependent reduced spaces. This generalizes classical model order reduction methods with time-independent reduced spaces. An a posteriori error estimate using the logarithmic Lipschitz constant associated to the flux has also been proposed and provides an efficient a posteriori error estimate. Using this error estimate, we have derived a greedy algorithm (called T-greedy algorithm) for the construction of a sequence of time-dependent reduced spaces. We have performed several numerical tests on both linear and nonlinear dynamical systems. The order reduction method with time-dependent spaces constructed with a T-greedy algorithm (MTD) has been compared to a method with time-independent spaces constructed with a POD-greedy algorithm (MTI). For the same online costs (time for evaluating the ROM), the error in the predictions of the ROM obtained by MTD may be several orders of magnitude lower than the error in the predictions of the ROM obtained by MTI. In other words, for a given precision of the ROM, the MTD approach yields a significant reduction of online computation costs. Here, the proposed T-greedy algorithm for the selection of parameters relies on an error indicator using a $L^2$-norm in time of an error estimate. This global error estimate could be improved by considering sharper local in time error estimates which could be obtained by improving the local Lipschitz constants used to derive our error bound [34]. Moreover, a goal-oriented variant of this algorithm could be introduced by considering norms taking into account a certain quantity of interest (e.g., the value of the solution at the final time) in order to improve the efficiency of the order reduction method for this approximation of this quantity of interest.

# References

[1] M. BARRAULT, Y. MADAY, N.C. NGUYEN AND A.T. PATERA, *An "empirical interpolation" method: application to efficient reduced-basis discretization of partial differential equations*. Comptes Rendus Mathématique, 339(9):667-672, 2004.

[2] U. BAUR, P. BENNER, B. HAASDONK, C. HIMPE, I. MAIERN AND M. OHLBERGER, *Comparison of methods for parametric model order reduction of instationary problems*, preprint, 2015.

[3] M. BEBENDORF, Y. MADAY AND B. STAMM, *Comparison of some Reduced Representation Approximations*, Quarteroni, Alfio and Rozza, Gianluigi, Reduced Order Methods for Modeling and Computational Reduction, MS&A - Modeling, Simulation and Applications 9:67-100, 2014.
[4] P. BENNER, S. GUGERCIN AND K. WILLCOX, *A Survey of Model Reduction Methods for Parametric Systems*, preprint, 2013.

[5] P. BINEV, A. COHEN, W. DAHMEN, R. DeVORE, G. PETROVA AND P. WOJ-TASZCZYK, *Convergence Rates for Greedy Algorithms in Reduced Basis Methods*. SIAM Journal on Mathematical Analysis, 43(3):1457–1472, 2011.

[6] A. BUFFA, Y. MADAY, A.T. PATERA, C. PRUD'HOMME AND G. TURINICI, *A priori convergence of the greedy algorithm for the parametrized reduced basis*, ESAIM-Math. Model. Numer. Anal., 46(3):595–603, 2012.

[7] T.T. BUI AND K. WILLCOX AND O. GHATTAS, *Model Reduction for Large-Scale Systems with High-Dimensional Parametric Input Space*, SIAM Journal on Scientific Computing, 6:3270–3288, 2008.

[8] S. CHATURANTABUT AND D.C. SORENSEN, *Nonlinear Model Reduction via Discrete Empirical Interpolation*, SIAM Journal on Scientific Computing, 42(5):2737–2764, 2010.

[9] M. CHENG, T.Y. HOU AND Z. ZHANG AND D.C. SORENSEN, *A dynamically bi-orthogonal method for time-dependent stochastic partial differential equations I: Derivation and algorithms*, Journal of Computational Physics, 242(0):843-868, 2013.

[10] M. CHENG, T.Y. HOU AND Z. ZHANG AND D.C. SORENSEN, *A dynamically bi-orthogonal method for time-dependent stochastic partial differential equations II: Adaptivity and generalizations*, Journal of Computational Physics, 242(0):753–776, 2013.

[11] R. DeVore, G. Petrova and P. Wojtaszczyk *Greedy algorithms for reduced bases in Banach spaces*, Constructive Approximation, 37(3):455–466, 2013.

[12] M. DROHMANN, B. HAASDONK AND M. OHLBERGER, *Reduced Basis Approximation for Nonlinear Parametrized Evolution Equations based on Empirical Operator Interpolation*, SIAM Journal on Scientific Computing, 34(2):A937–A969, 2012.

[13] A posteriori error bounds for the empirical interpolation method, M NU EX A.M. GREPL AND A.T. PATERA, Comptes Rendus Mathématique, 348(9-10):575-579,2010.

[14] J.L. EFTANG, D.J. KNEZEVIC AND A.T. PATERA, *An hp certified reduced basis method for parametrized parabolic partial differential equations*, Mathematical and Computer Modelling of Dynamical Systems, 17(4):395–422, 2011.

[15] S. GLAS, A. MAYERHOFER, K. URBAN, *Two Ways to Treat Time in Reduced Basis Methods* (preprint), 2016
[4] P. BENNER, S. GUGERCIN AND K. WILLCOX, *A Survey of Model Reduction Methods for Parametric Systems*, preprint, 2013.

[5] P. BINEV, A. COHEN, W. DAHMEN, R. DeVORE, G. PETROVA AND P. WOJTASZCZYK, *Convergence Rates for Greedy Algorithms in Reduced Basis Methods*. SIAM Journal on Mathematical Analysis, 43(3):1457–1472, 2011.

[6] A. BUFFA, Y. MADAY, A.T. PATERA, C. PRUD'HOMME AND G. TURINICI, *A priori convergence of the greedy algorithm for the parametrized reduced basis*, ESAIM-Math. Model. Numer. Anal., 46(3):595–603, 2012.

[7] T.T. BUI AND K. WILLCOX AND O. GHATTAS, *Model Reduction for Large-Scale Systems with High-Dimensional Parametric Input Space*, SIAM Journal on Scientific Computing, 6:3270–3288, 2008.

[8] S. CHATURANTABUT AND D.C. SORENSEN, *Nonlinear Model Reduction via Discrete Empirical Interpolation*, SIAM Journal on Scientific Computing, 42(5):2737–2764, 2010.

[9] M. CHENG, T.Y. HOU AND Z. ZHANG AND D.C. SORENSEN, *A dynamically bi-orthogonal method for time-dependent stochastic partial differential equations I: Derivation and algorithms*, Journal of Computational Physics, 242(0):843-868, 2013.

[10] M. CHENG, T.Y. HOU AND Z. ZHANG AND D.C. SORENSEN, *A dynamically bi-orthogonal method for time-dependent stochastic partial differential equations II: Adaptivity and generalizations*, Journal of Computational Physics, 242(0):753–776, 2013.

[11] R. DeVore, G. Petrova and P. WOJTASZCZYK *Greedy algorithms for reduced bases in Banach spaces*, Constructive Approximation, 37(3):455–466, 2013.

[12] M. DROHMANN, B. HAASDONK AND M. OHLBERGER, *Reduced Basis Approximation for Nonlinear Parametrized Evolution Equations based on Empirical Operator Interpolation*, SIAM Journal on Scientific Computing, 34(2):A937–A969, 2012.

[13] A posteriori error bounds for the empirical interpolation method, M NU EX A.M. GREPL AND A.T. PATERA, Comptes Rendus Mathématique, 348(9-10):575-579,2010.

[14] J.L. EFTANG, D.J. KNEZEVIC AND A.T. PATERA, *An hp certified reduced basis method for parametrized parabolic partial differential equations*, Mathematical and Computer Modelling of Dynamical Systems, 17(4):395–422, 2011.

[15] S. GLAS, A. MAYERHOFER, K. URBAN, *Two Ways to Treat Time in Reduced Basis Methods* (preprint), 2016.
[16] M.A. GREPL AND A.T. PATERA, *A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations*, ESAIM: Mathematical Modelling and Numerical Analysis, 39(01):157-181,2005.

[17] B. HAASDONK M. OHLBERGER, *Reduced basis method for finite volume approximations of parametrized linear evolution equations*, ESAIM: Mathematical Modelling and Numerical Analysis, 42(2):277-302, 2008.

[18] B. HAASDONK M. OHLBERGER, *Efficient reduced models and a posteriori error estimation for parametrized dynamical systems by offline/online decomposition*, Mathematical and Computer Modelling of Dynamical Systems, 17(2):145-161,2011

[19] B. HAASDONK, *Convergence Rates of the POD-Greedy method*, M2AN Math. Model. Numer. Anal., 47:859-873, 2015.

[20] B. HAASDONK, *Reduced basis methods for parametrized PDEs - A tutorial introduction for stationary and instationary problems.*, Tech. Rep., 2014, Chapter to appear in P. Benner, A. Cohen, M. Ohlberger and K. Willcox: "Model reduction and approximation for complex systems", Springer.

[21] A. JANON, M. NODET AND C. PRIEUR, *Certified reduced-basis solutions of viscous Burgers equation parametrized by initial and boundary values*, ESAIM: M2AN, 47(2):317-348,2013.

[22] M. KARCHER AND M.A. GREPL, *A Posteriori Error Estimation for Reduced Order Solutions of Parametrized Parabolic Optimal Control Problems*, ESAIM: M2AN, 48(6):1615-1638, 2014.

[23] O. KOCH AND C. LUBICH, *Dynamical Low-Rank Approximation*, SIAM Journal on Matrix Analysis and Applications, 29(2):434-454, 2007.

[24] E. MUSHARBASH, F. NOBILE AND T. ZHOU, *On the dynamically orthogonal approximation of time-dependent random PDEs*, SIAM Journal on Scientific Computing, 2015, 37(2):A776-A810, 2015

[25] A. NOUY, *Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms*. Computer Methods in Applied Mechanics and Engineering, 197:4718-4736, 2008.

[26] A. NOUY, *Proper Generalized Decompositions and Separated Representations for the Numerical Solution of High Dimensional Stochastic Problems*, Archives of Computational Methods in Engineering, 17(4):403-434, 2010.

[27] A. NOUY, *Low-rank tensor methods for model order reduction*, ArXiv e-prints, November 2015.
[16] M.A. GREPL AND A.T. PATERA, *A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations*, ESAIM: Mathematical Modelling and Numerical Analysis, 39(01):157-181, 2005.

[17] B. HAASDONK M. OHLBERGER, *Reduced basis method for finite volume approximations of parametrized linear evolution equations*, ESAIM: Mathematical Modelling and Numerical Analysis, 42(2):277-302, 2008.

[18] B. HAASDONK M. OHLBERGER, *Efficient reduced models and a posteriori error estimation for parametrized dynamical systems by offline/online decomposition*, Mathematical and Computer Modelling of Dynamical Systems, 17(2):145-161, 2011

[19] B. HAASDONK, *Convergence Rates of the POD-Greedy method*, M2AN Math. Model. Numer. Anal., 47:859-873, 2015.

[20] B. HAASDONK, *Reduced basis methods for parametrized PDEs - A tutorial introduction for stationary and instationary problems.*, Tech. Rep., 2014, Chapter to appear in P. Benner, A. Cohen, M. Ohlberger and K. Willcox: "Model reduction and approximation for complex systems", Springer.

[21] A. JANON, M. NODET AND C. PRIEUR, *Certified reduced-basis solutions of viscous Burgers equation parametrized by initial and boundary values*, ESAIM: M2AN, 47(2):317-348, 2013.

[22] M. KARCHER AND M.A. GREPL, *A Posteriori Error Estimation for Reduced Order Solutions of Parametrized Parabolic Optimal Control Problems*, ESAIM: M2AN, 48(6):1615-1638, 2014.

[23] O. KOCH AND C. LUBICH, *Dynamical Low-Rank Approximation*, SIAM Journal on Matrix Analysis and Applications, 29(2):434-454, 2007.

[24] E. MUSHARBASH, F. NOBILE AND T. ZHOU, *On the dynamically orthogonal approximation of time-dependent random PDEs*, SIAM Journal on Scientific Computing, 2015, 37(2):A776-A810, 2015

[25] A. NOUY, *Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms*. Computer Methods in Applied Mechanics and Engineering, 197:4718-4736, 2008.

[26] A. NOUY, *Proper Generalized Decompositions and Separated Representations for the Numerical Solution of High Dimensional Stochastic Problems*, Archives of Computational Methods in Engineering, 17(4):403-434, 2010.

[27] A. NOUY, *Low-rank tensor methods for model order reduction*, ArXiv e-prints, November 2015.
[28] N.C. NGUYEN, G. ROZZA AND A.T. PATERA, *Reduced basis approximation and a posteriori error estimation for the time-dependent viscous Burgers' equation*, Calcolo, 46(3):157-185, 2009.

[29] A. QUARTERONI, R. SACCO AND F. SALERI, *Numerical Mathematics, Texts in Applied Mathematics*, Springer, 2012.

[30] T.P. SAPSIS AND P.F.J. LERMUSIAUX, *Dynamically orthogonal field equations for continuous stochastic dynamical systems*, Physica D: Nonlinear Phenomena, 238(23-24):2347-2360, 2009.

[31] K. STEIH AND K. URBAN, *Space-time Reduced Basis methods for time-periodic parametric partial differential equations*, Ulm University, 2011.

[32] K. URBAN AND A.T. PATERA, *A new error bound for reduced basis approximation of parabolic partial differential equations*, Comptes Rendus Mathématique, 350(3): 203-207, 2012.

[33] S. VOLKWEIN, *Model reduction using proper orthogonal decomposition*, Lecture notes, 2008.

[34] D. WIRTZ, D.C. SORENSEN AND B. HAASDONK, *A-posteriori error estimation for DEIM reduced nonlinear dynamical systems*, SIAM J. Sci. Comput., 36(2):A311-A338, 2014.

[35] M. YANO, A.T. PATERA AND K. URBAN, *A Space-Time Hp-Interpolation-Based Certified Reduced Basis Method for Burgers' Equation*, Models Methods Appl. Sci., 24(9):1903-1935, 2014.
[28] N.C. NGUYEN, G. ROZZA AND A.T. PATERA, *Reduced basis approximation and a posteriori error estimation for the time-dependent viscous Burgers' equation*, Calcolo, 46(3):157-185, 2009.

[29] A. QUARTERONI, R. SACCO AND F. SALERI, *Numerical Mathematics, Texts in Applied Mathematics*, Springer, 2012.

[30] T.P. SAPSIS AND P.F.J. LERMUSIAUX, *Dynamically orthogonal field equations for continuous stochastic dynamical systems*, Physica D: Nonlinear Phenomena, 238(23-24):2347-2360, 2009.

[31] K. STEIH AND K. URBAN, *Space-time Reduced Basis methods for time-periodic parametric partial differential equations*, Ulm University, 2011.

[32] K. URBAN AND A.T. PATERA, *A new error bound for reduced basis approximation of parabolic partial differential equations*, Comptes Rendus Mathématique, 350(3): 203-207, 2012.

[33] S. VOLKWEIN, *Model reduction using proper orthogonal decomposition*, Lecture notes, 2008.

[34] D. WIRTZ, D.C. SORENSEN AND B. HAASDONK, *A-posteriori error estimation for DEIM reduced nonlinear dynamical systems*, SIAM J. Sci. Comput., 36(2):A311-A338, 2014.

[35] M. YANO, A.T. PATERA AND K. URBAN, *A Space-Time Hp-Interpolation-Based Certified Reduced Basis Method for Burgers' Equation*, Models Methods Appl. Sci., 24(9):1903-1935, 2014.