# Generalization of the Neville-Aitken interpolation algorithm on Grassmann manifolds: Applications to reduced order model

Rolando Mosquera, Abdallah El Hamidi, Aziz Hamdouni, Antoine Falaize

► To cite this version:

Rolando Mosquera, Abdallah El Hamidi, Aziz Hamdouni, Antoine Falaize. Generalization of the Neville-Aitken interpolation algorithm on Grassmann manifolds: Applications to reduced order model. International Journal for Numerical Methods in Fluids, 2021, 93 (7), pp.2421-2442. 10.1002/fld.4981. hal-05335379

HAL Id: hal-05335379

https://hal.univ-brest.fr/hal-05335379v1

Submitted on 28 Oct 2025

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
<figure><img src="image_2.png" /></figure>

# Generalization of the Neville-Aitken interpolation algorithm on Grassmann manifolds: Applications to reduced order model

Rolando Mosquera<fnref n="1" />,<fnref n="2" /> ID |

Abdallah El Hamidi<fnref n="2" /> |

Aziz Hamdouni<fnref n="2" /> |

Antoine Falaize<fnref n="2" />

## Correspondence

Rolando Mosquera, IRDL, UMR-6027-CNRS, Université Bretagne Sud, Lorient, France.

Email: rolando.mosquera_meza@univ-lr.fr

## Present address

Rolando Mosquera, Laboratoire des Sciences de l'Ingénieur pour l'Environnement (LaSIE), La Rochelle Université, Pôle Science et Technologie, Avenue Michel Crépeau La Rochelle Cedex 1, 17042, France

## Abstract

An extension of the well-known Neville-Aitken's algorithm for interpolation on the Grassmann manifold $G_m(\mathbb{R}^n)$ in the framework of parametric model order reduction is presented. Interpolation points on $G_m(\mathbb{R}^n)$ are the subspaces spanned by bases obtained by Proper Orthogonal Decomposition of available solutions associated with the chosen parameter sampling. The Neville-Aitken's algorithm is performed recursively via the geodesic barycenter of two points. Three CFD applications are presented: (i) the Von Karman vortex shedding street, (ii) the lid-driven cavity with inflow and (iii) the flow induced by a rotating solid. Numerical results are relevant with respect to accuracy while the asymptotic complexity is comparable to the state of the art.

## KEYWORDS

Grassmann manifold, Neville interpolation algorithm, proper orthogonal decomposition

# 1 | INTRODUCTION

Model order reduction techniques have undergone important developments during the past 20 years. Their use in numerous engineering problems have made it essential tools for quasi-real-time simulation and control of complex problems (see References 1-4). Notable applications are fluid mechanics,⁵⁻⁹ fluid-structure interaction,¹⁰⁻¹³ and structural or material mechanics.¹⁴⁻¹⁶ Basically, it consists in the approximation of a problem with a very large number of degrees of freedom (typically several thousands to millions) by a model with at most one hundred of degrees of freedom.¹⁷⁻¹⁹ The methodology can be decomposed in two main steps:

1. the construction of an appropriate approximation space of the solution,

2. the projection of the governing equations on this approximation space (by Galerkin projection, Petrov-Galerkin projection, residual minimization or any other projection method).

Two main classes of model order reduction methods can be distinguished:

* A posteriori reduction, where the approximation space is constructed from snapshots of the numerical solutions of the full-order problem or from experimental data. Typical examples include the Proper Orthogonal Decomposition (POD).

<footnote n="1">¹IRDL, UMR-6027-CNRS, Université Bretagne Sud, Lorient, France</footnote>

<footnote n="2">²Équipe M2N, LaSIE, UMR-7356-CNRS, La Rochelle Université, La Rochelle, France</footnote>
* A priori reduction, where the two steps above are coupled, that is, the approximation space is determined iteratively in order to minimize the estimated error between the reduced-order solution and the solution of the full-order model. This is typically the Proper Generalized Decomposition (PGD, see References 20,21), the a priori reduction associated with the POD (APR-POD, see References 22,23) or the Reduced Bases method (RB, see References 24,25), though the latter can be seen as mixed with a priori selection of interpolation points and a posteriori computation of the reduced solution.

All these methods can be applied to a given problem to derive an efficient reduced-order model. However, the governing equations usually depend on a parameter of interest (e.g., physical properties, geometry) so that one is interested in the construction of parametric reduced-order models. In this context, three approaches are generally considered in association with the techniques presented above:

1. Construction of parametric reduced-order models by the PGD. With this approach, the reduced solution is searched by a variable splitting method applied to the space, the time and the parameters (see Reference 20 for details), allowing to construct so-called numerical vademecum,⁸,²¹ but raising several difficulties. First, the PGD is computationally intensive for problems in fluid mechanics (see References 26-29). Second, the computation of the solution for parameters not included in the initial discretization of the parameter space requires to interpolate in the approximation space constructed by tensor products of vector subspaces, whose mathematical structure is complex (see Reference 30).

2. Construction of parametric reduced-order models by the RB method. With this approach, several solutions are computed in an offline phase in association with parameters selected by a greedy algorithm. Then, the solution for a given parameter is constructed at the online phase by linear combination of the available solutions. This approach assumes that the solution space is of low complexity as measured by the Kolmogorov widths (see References 31-33). Generally, the nonlinear components are approximated by the Empirical Interpolation Method (EIM introduced in Reference 34) or its discrete version (DEIM, introduced in Reference 35). However, this approach does not appear to be efficient for fluid mechanics problems (see Reference 36).

3. Construction of parametric reduced-order models by the POD. With this approach, several POD bases are computed from the numerical solutions of the governing equations associated with sampling parameters and the POD basis associated with a new parameter is computed from the interpolation into the set of sampling POD bases. The reduced order model (ROM) is then constructed by classical projection of the governing equations on this interpolated basis. This approach is widely used in the literature¹³,³⁷⁻⁴¹ and appears to be efficient for most problems.

The present work is related with the latter approach as follows. We consider the system of equations obtained by a semi-discretization (in space) of a parametrized evolution problem described by a partial differential equation:

$$ \begin{cases} \frac{du_h}{dt} = \mathcal{F}(u_h, \lambda), & \forall t \ge 0, \\ u_h(0, \lambda) = u_h^0, \end{cases} \qquad (1) $$

where $u_h(t, \lambda) \in \mathbb{R}^n$ denotes the solution, $n$ denotes the number of degrees of freedom of the spatial discretization, $\lambda \in \mathbb{R}^p$ denotes the parameter, and $u_h^0$ denotes an initial condition. Suppose that, for a given parameter $\lambda$, we know an orthonormal family $(\phi_1(\lambda), \dots, \phi_m(\lambda))$ of $m$ vectors in $\mathbb{R}^n$. We denote $\Phi(\lambda) = (\phi_1(\lambda), \dots, \phi_m(\lambda))$ the matrix in $\mathbb{R}^{n \times m}$ whose columns are the $\phi_i(\lambda)$, $1 \le i \le m$. Then the solution of (1) is approximated by:

$$ u_h(t, \lambda) \approx \Phi(\lambda) \cdot a(t), \qquad (2) $$

where the vector $a(t) \in \mathbb{R}^m$ denotes the temporal coefficients of the projection of $u_h(t, \lambda)$ on the $m$-dimensional subspace spanned by the vectors $\phi_i(\lambda)$, $1 \le i \le m$. A ROM associated with (1) is then

$$ \begin{cases} \frac{da}{dt} = \hat{\mathcal{F}}(a, \lambda), & \forall t \ge 0, \\ a(0) = \Phi(\lambda)^T \cdot u_h^0, \end{cases} \qquad (3) $$

with

$$ \hat{\mathcal{F}}(a, \lambda) = \Phi(\lambda)^T \cdot \mathcal{F}(\Phi(\lambda) \cdot a, \lambda). \qquad (4) $$
<figure><img src="image_3.png" /><figcaption>FIGURE 1 Naive interpolation that does not preserve the intrinsic structure of $G_1(\mathbb{R}^3)$, that is, the one-dimensional subspaces in $\mathbb{R}^3$. [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

The construction of the matrix $\Phi(\lambda)$ by POD requires a sequence of solutions of the full order (1), as detailed in Section 2 along with the main steps in the construction of the ROM. Now, provided a set of N POD bases $\Phi(\lambda_k)$, $0 \le k \le N$ associated with sampling parameters $\lambda_k$, $0 \le k \le N$, the construction of the ROM (3) for a new parameter $\lambda^* \ne \lambda_k$, $0 \le k \le N$, without the costly computation of the snapshots for (1) is as follows. First, an approximation $\Phi^*$ to the basis $\Phi(\lambda^*)$ is built by interpolating into a set of known bases $\Phi(\lambda_k)$, $0 \le k \le N$, and second, this new basis is used in equation (4) to construct the reduced vector field $\hat{\mathcal{F}}^*(a, \lambda^*) = \Phi^{*1} \cdot \mathcal{F}(\Phi^* \cdot a, \lambda^*)$ so that the ROM for the new parameter is given by

$$ \left\{ \begin{array}{l} \frac{da}{dt} = \hat{\mathcal{F}}^*(a, \lambda^*), \quad \forall t \ge 0, \\ a(0) = \Phi^{*1} \cdot u_h^0. \end{array} \right. \qquad (5) $$

A natural question arising in this context is how to interpolate the POD bases. In a previous work,³⁷ the authors pointed out that the appropriate objects to interpolate are not the POD bases but the underlying spanned subspaces. The set of all m-dimensional vector spaces in $\mathbb{R}^n$ is a compact Riemannian manifold, called the Grassmann manifold and denoted $G_m(\mathbb{R}^n)$. A naive interpolation is not compatible with the intrinsic structure of $G_m(\mathbb{R}^n)$, as illustrated in Figure 1 for one-dimensional subspaces in $\mathbb{R}^3$. For a new parameter value $\lambda^*$, the interpolation procedure introduced in Reference 37 is as follows. First, we choose a reference point $\bar{\Phi}_{\text{ref}} := \text{span}(\varphi_1(\lambda_{\text{ref}}), \dots, \varphi_m(\lambda_{\text{ref}}))$ on $G_m(\mathbb{R}^n)$ and denote the tangent space to $G_m(\mathbb{R}^n)$ at $\bar{\Phi}_{\text{ref}}$ by $T_{\bar{\Phi}_{\text{ref}}} G_m(\mathbb{R}^n)$. Second, all the points $(\bar{\Phi}_k)_{0 \le k \le N}$ are mapped on the tangent space $T_{\bar{\Phi}_{\text{ref}}} G_m(\mathbb{R}^n)$, where $\bar{\Phi}_k := \text{span}(\varphi_1(\lambda_k), \dots, \varphi_m(\lambda_k))$, for $0 \le k \le N$. This can be achieved thanks to the logarithm mapping at the reference point $\bar{\Phi}_{\text{ref}}$, denoted by $\log_{\bar{\Phi}_{\text{ref}}} : G_m(\mathbb{R}^n) \to T_{\bar{\Phi}_{\text{ref}}} G_m(\mathbb{R}^n)$, well defined in a neighborhood of the point $\bar{\Phi}_{\text{ref}}$. Third, a classical vector interpolation is performed in the vector space $T_{\bar{\Phi}_{\text{ref}}} G_m(\mathbb{R}^n)$ for the target parameter $\lambda^*$. Finally, one returns to the manifold $G_m(\mathbb{R}^n)$ via the reciprocal exponential mapping from the tangent space $T_{\bar{\Phi}_{\text{ref}}} G_m(\mathbb{R}^n)$ to get the desired approximation of $\bar{\Phi}^* = \text{span}(\varphi_1(\lambda^*), \dots, \varphi_m(\lambda^*))$ on $G_m(\mathbb{R}^n)$. This method has been successfully applied to several problems from structural mechanics to aeroelasticity in References 37-39. However, the interpolation results can be sensitive to the choice of the reference point.⁴¹ The strategy that consists in the use of a reference point to construct an approximation to a set of points on a Riemannian manifold has been considered in other works.⁴²,⁴³ In a recent work,⁴¹ the authors developed the Grassmann Inverse Distance Weighting (G-IDW) interpolation algorithm as an extension of the well-known Inverse Distance Weighting (IDW) method to Grassmann manifolds. It consists in the minimization of a quadratic function in the geodesic distance to interpolation points, with appropriate weights. This G-IDW method does not require a reference point and yields relevant results. However, it relies on an iterative procedure. The objective of the work presented in this paper is to develop an explicit algorithm to interpolate over the Grassmann manifold without relying on a reference vector space tangent to the manifold at a selected point.

In this work, we propose an extension of the well-known Neville-Aitken's algorithm⁴⁴,⁴⁵ to construct an explicit interpolation algorithm on the Grassmann manifold without relying on a reference vector space. In vector spaces, the Neville-Aitken algorithm computes the Lagrange interpolation polynomial recursively from the interpolation of two points. Replacing straight lines by geodesics, we extend this algorithm to Grassmann manifolds. The resulting algorithm has several interesting features. First, it does not rely on the choice of a reference vector space so that the set of sampling points is not limited to a neighborhood of a particular reference point, which is a necessary condition to guarantee the existence of the inverse of the exponential map in the method proposed in reference.³⁷ Second, it is explicit (i.e., noniterative) which ensures that an interpolated basis can be computed in finite time, contrarily to the method proposed in Reference 41.

This paper is organized as follows. Reminders on the POD-Galerkin method are given in Section 2. In Section 3, the Grassmann manifold and the underlying tools used in the paper are briefly presented. In Section 4, the Neville-Aitken's
algorithm in the framework of vector spaces is recalled and its extension to Grassmann manifolds is introduced. The performances of the proposed method are illustrated in Section 5 through three independent CFD applications, namely: the Von Karman vortex shedding street, the lid-driven cavity with inflow and the flow induced by a rotating solid. Our results are compared to those obtained by the exact POD method and to the interpolation procedure developed in Reference 37. A conclusion is presented in Section 6, where the benefits and the limitations of our work are highlighted and interesting perspectives are proposed.

## 2 | REMINDERS ON THE POD-GALERKIN METHOD FOR THE CONSTRUCTION OF ROMS

The matrix $\Phi(\lambda)$ in the framework of POD is determined from a set of $n'$ snapshots of the approximate solution to the governing equations (1) at chosen instants $t_i$, $1 \le i \le n'$ as follows (see also²⁵ for a detailed presentation). Denote $S(\lambda) = (u_h(t_1, \lambda), \dots, u_h(t_{n'}, \lambda)) \in \mathbb{R}^{n' \times n'}$ the matrix whose columns are the selected snapshots and consider the thin Singular Value Decomposition (SVD) of $S(\lambda)$:

$$S(\lambda) = \tilde{\Phi}(\lambda) \cdot \Sigma(\lambda) \cdot V^{\top}(\lambda), \qquad (6)$$

where $\tilde{\Phi}(\lambda) \in \mathbb{R}^{n' \times n'}$ is the matrix whose columns are the left-singular orthogonal unit vectors of $S(\lambda)$, $\Sigma(\lambda) \in \mathbb{R}^{n' \times n'}$ is the diagonal matrix whose entries are the singular values of $S(\lambda)$ with $\Sigma(\lambda) = \text{diag}(\sigma_1, \dots, \sigma_{n'})$ and $\sigma_1 \ge \dots \ge \sigma_{n'}$, and $V(\lambda) \in \mathbb{R}^{n' \times n'}$ is the matrix whose columns are the right-singular orthogonal unit vectors of $S(\lambda)$. Denoting $\Phi(\lambda) = (\varphi_1(\lambda), \dots, \varphi_m(\lambda)) \in \mathbb{R}^{n \times m}$ the truncation of $\tilde{\Phi}(\lambda)$ to the first $m$ columns, we have the classical approximation result:

$$\|S - \Phi \cdot \Phi^{\top} \cdot S\|_{F} = \sum_{i=m+1}^{n'} \sigma_{i}^{2}, \qquad (7)$$

where $\| \cdot \|_{F}$ denotes the Frobenius norm. The truncation index $m$ can be chosen to achieve a given precision $\epsilon \ll 1$ so that:

$$\frac{\sum_{i=m+1}^{n'} \sigma_i^2}{\sum_{i=1}^{n'} \sigma_i^2} \le \epsilon. \qquad (8)$$

For a given parameter $\lambda$, the construction of the ROM by the POD-Galerkin method is done in two steps.

Step 1: Construction of the POD basis. Given the matrix of snapshots $S(\lambda) = (u_h(t_1, \lambda), \dots, u_h(t_{n'}, \lambda)) \in \mathbb{R}^{n \times n'}$ of numerical solutions so the governing equations (1), the POD basis $\Phi(\lambda)$ is computed to a given precision $\epsilon \ll 1$ by the truncated SVD in Equation (6) where the truncation index can be selected according to equation (8).

Step 2: Construction of the ROM. Defining the reduced order solution as $\hat{u}(t, \lambda) = \Phi(\lambda) a(t)$ with $a(t) \in \mathbb{R}^m$, the reduced order vector field is computed by

$$\hat{\mathcal{F}}(a, \lambda) = \Phi(\lambda)^{\top} \cdot \mathcal{F}(\Phi(\lambda) \cdot a, \lambda),$$

leading to the following ROM

$$\left\{ \begin{array}{ll} \frac{da}{dt} = \hat{\mathcal{F}}(a, \lambda), & \forall t \ge 0, \\ a(0) = \Phi(\lambda)^{\top} \cdot u_h^0. \end{array} \right.$$

Now, given a new parameter $\lambda^*$ for which the snapshots $S(\lambda^*) = (u_h(t_1, \lambda^*), \dots, u_h(t_{n'}, \lambda^*)) \in \mathbb{R}^{n \times n'}$ are not known, the construction of the POD basis in step 1 above is replaced by an interpolation in the set of known reduced basis associated with the sampling parameters $\lambda_1, \dots, \lambda_N$ to provide the interpolated basis $\Phi^*$. Following References 37,38, this interpolation is carried out on the Grassmann manifold. The next section is devoted to a recall of the mathematical tools necessary to construct such Grassmann interpolation algorithms.
## 3 | THE GRASSMANN MANIFOLD

Let $m, n$ be two positive integers with $m \le n$. The set of all $\mathbb{R}$-subspaces of $\mathbb{R}^n$ of dimension $m$ is a differentiable manifold of dimension $m \times (n-m)$, called the Grassmann manifold and usually denoted by $G_m(\mathbb{R}^n)$. Let $\mathbb{R}_*^{n \times m}$ be the set of all matrices of size $n \times m$, whose column vectors are linearly independent; one defines on this set the following equivalence relation

$$X \sim Y, \text{ if and only if } \exists A \in \mathrm{GL}_m(\mathbb{R}) \text{ such that } X = YA,$$

where $\mathrm{GL}_m(\mathbb{R})$ denotes the linear group of $\mathbb{R}^m$, that is, the set of $m \times m$ invertible matrices. The Grassmann manifold can be realized as the set of such equivalence classes. More precisely,

$$G_m(\mathbb{R}^n) = \{ \bar{X} : X \in \mathbb{R}_*^{n \times m} \},$$

where

$$\bar{X} = \{ X A : A \in \mathrm{GL}_m(\mathbb{R}) \}.$$

The Grassmann manifold is a compact topological space and, thanks to the canonical projection denoted by

$$\pi : \mathbb{R}_*^{n \times m} \to G_m(\mathbb{R}^n),$$

it can be endowed by a differentiable manifold structure so that $\pi$ is a submersion (see Reference 46).

### 3.1 | Geodesic and exponential map

The set $\mathbb{R}_*^{n \times m}$ can be endowed with a metric $g$ that is right-invariant by $\mathrm{GL}_m(\mathbb{R})$, that is, such that

$$g_{XA}(Y_1A, Y_2A) = g_X(Y_1, Y_2).$$

This metric $g$ is defined by

$$g_X(Y_1, Y_2) = \mathrm{tr} \left[ (X^\top X)^{-1} Y_1^\top Y_2 \right], \qquad (9)$$

for every $Y_1, Y_2 \in \mathbb{R}^{n \times m}$, $X \in \mathbb{R}_*^{n \times m}$ and $A \in \mathrm{GL}_m(\mathbb{R})$. It induces a Riemannian metric $\bar{g}$ on $G_m(\mathbb{R}^n)$ through the submer-
sion $\pi$. Moreover, one constructs a Riemannian connection on $G_m(\mathbb{R}^n)$, which allows the definition of geodesics and
the geodesic exponential mapping on $G_m(\mathbb{R}^n)$ as detailed in References 47,48. In what follows, the distance on $G_m(\mathbb{R}^n)$
induced by the metric $\bar{g}$ will be denoted by $d$, and the ball centered at $y$ of radius $r > 0$ in $G_m(\mathbb{R}^n)$ will be denoted
by $B_d(y, r)$.

For every $X \in \mathbb{R}_*^{n \times m}$, we introduce $H_X = [T_X \pi^{-1}(\bar{X})]^{\perp}$, the orthogonal space to the tangent space at $X$ of the fiber $\pi^{-1}(\bar{X})$, where the orthogonality is taken with respect to the scalar product induced by the metric $g$. Then the mapping

$$d_X \pi : H_X \to T_{\bar{X}} G_m(\mathbb{R}^n),$$

is an isomorphism. Now, given a vector $v \in T_{\bar{X}} G_m(\mathbb{R}^n)$, if we perform the thin SVD of the matrix $(d_X \pi)^{-1}(v) (X^\top X)^{-1/2}$ as

$$(d_X \pi)^{-1}(v) (X^\top X)^{-1/2} = U \Sigma V^\top,$$

the geodesic of initial point $\gamma(0) = \bar{X}$ and initial velocity $\gamma'(0) = v$ can be written in the following form

$$\gamma(t) = \overline{X(X^\top X)^{-1/2} V \cos(t\Sigma) + U \sin(t\Sigma)}.$$
Thanks to the compactness of the Grassmann manifold, the geodesics are defined on all $\mathbb{R}$ by the Hopf–Rinow theorem. The geodesic exponential mapping is defined by $\exp_{\bar{X}}(v) = \gamma(1)$ and has a matrix representation given by:

$$X(X^\top X)^{-1/2} V \cos(\Sigma) + U \sin(\Sigma).$$

## 3.2 | Injectivity radius and logarithm map

Thereafter we need the inverse of the exponential application. First, we can notice that $d_{0\bar{X}}\exp_{\bar{X}} = \text{Id}$, so that the theorem of the inverse function ensures that $\exp_{\bar{X}}$ is a local diffeomorphism at $0_{\bar{X}}$ (the null vector of $T_{\bar{X}}G_m(\mathbb{R}^n)$). To study more precisely the domain where the exponential map admits an inverse, we introduce the notion of radius of injectivity in the following definition:

**Definition 1** (Injectivity radius). Let $\bar{X} \in G_m(\mathbb{R}^n)$. The injectivity radius of $(G_m(\mathbb{R}^n), \bar{g})$ at $\bar{X}$ is defined by

$$r_I(\bar{X}) = \sup \{ \rho > 0 : \exp_{\bar{X}} : B(0_{\bar{X}}, \rho) \to \exp_{\bar{X}}(B(0_{\bar{X}}, \rho)) \text{ is a diffeomorphism} \},$$

and the injectivity radius of $G_m(\mathbb{R}^n)$ is defined by

$$r_I = \inf\{r_I(\bar{X}) : \bar{X} \in G_m(\mathbb{R}^n)\}.$$

For example, the injectivity radius of the unit sphere $S^n$ is $\pi$ and that of the Euclidean space is infinity. The value of the injectivity radius of the Grassmann manifolds $G_m(\mathbb{R}^n)$ is given by the following theorem 1. which proof can be found in Reference 49.

**Theorem 1.** Let $n, m$ be two positive integers such that $\min\{m, n - m\} \ge 2$. Then the injectivity radius of $G_m(\mathbb{R}^n)$ is given by $r_I = \pi/2$.

This theorem gives the radius of the ball where the reciprocal mapping of the geodesic exponential $\exp_{\bar{X}}^{-1}$ (or $\log_{\bar{X}}$) is defined. So, for $\bar{Y} \in G_m(\mathbb{R}^n)$ with $d(\bar{X}, \bar{Y}) < \pi/2$, the reciprocal mapping $\exp_{\bar{X}}^{-1}(\bar{Y})$ is represented by the matrix

$$U_{\alpha}\tan^{-1}(\Sigma_{\alpha}) V_{\alpha}^{\top} (X^{\top}X)^{1/2},$$

with the thin SVD

$$U_{\alpha}\Sigma_{\alpha}V_{\alpha}^{\top} = [Y(X^{\top}Y)^{-1}(X^{\top}X) - X](X^{\top}X)^{-1/2}.$$

We can deduce that the geodesic between $\bar{X}$ and $\bar{Y}$ in $G_m(\mathbb{R}^n)$ is given by

$$\gamma_{\bar{X}, \bar{Y}}(t) = X (X^{\top} X)^{-1/2} V_{\alpha} \cos(t \tan^{-1}(\Sigma_{\alpha})) + U_{\alpha} \sin(t \tan^{-1}(\Sigma_{\alpha})). \quad (10)$$

After this reminders, we introduce in the next section an extension to the Grassmann manifold of the well-known Neville–Aitken algorithm.45,50

# 4 | NEVILLE-AITKEN’S ALGORITHM FOR LAGRANGE INTERPOLATION

In this section, we first recall the classical Neville–Aitken’s algorithm for vector-valued functions, and second introduce our extension to the Grassmann manifolds.

## 4.1 | Recall of the classical Neville–Aitken’s algorithm

Consider a function $f: \Lambda \subset \mathbb{R} \to V$, where $V$ is a vector space, and a subset $\{\lambda_0, \lambda_1, \dots\} \subset \Lambda$, with $f(\lambda_k) = y_k \in V$, for every $\lambda_k$ and $\lambda_i \ne \lambda_j$ for all $i \ne j$. Given an integer $N \ge 1$, the Lagrange interpolation polynomial of $f$ at
{ λ₀, λ₁, ..., λₙ } will be denoted by L[λ₀, ..., λₙ; y₀, ..., yₙ]. It is the unique polynomial of degree less or equal to N satisfying:

$$L[\lambda_0, \dots, \lambda_N; y_0, \dots, y_N](\lambda_k) = y_k, \quad \text{for every } k \in \{0, 1, \dots, N\}. \tag{11}$$

Let us consider the affine functions α[i,j] defined on Λ by

$$\alpha_{[i,j]}(\lambda) = \frac{\lambda - \lambda_i}{\lambda_j - \lambda_i}, \tag{12}$$

for all i,j∈{0, 1, ...} with i≠j. Then, the Neville-Aitken's algorithm can be summarized by the following induction formula:

$$L[\lambda_0, \dots, \lambda_N; y_0, \dots, y_N](\lambda) = (1 - \alpha_{[0,N]}(\lambda)) L[\lambda_0, \dots, \lambda_{N-1}; y_0, \dots, y_{N-1}](\lambda) \\ + \alpha_{[0,N]}(\lambda) L[\lambda_1, \dots, \lambda_N; y_1, \dots, y_N](\lambda), \tag{13}$$

with the initialization

$$L[\lambda_k; y_k](\lambda) = y_k \quad \text{for every } k \in \{0, 1, \dots, N\}. \tag{14}$$

That is, for a fixed λ ∈ Λ = [min₀≤k≤N λₖ, max₀≤k≤N λₖ], the quantity L[λ₀, ..., λₙ; y₀, ..., yₙ](λ) is the barycenter of L[λ₀, ..., λₙ-₁; y₀, ..., yₙ-₁](λ) and L[λ₁, ..., λₙ; y₁, ..., yₙ](λ), weighted by (1 - α[0,N](λ)) and α[0,N](λ), respectively. Therefore, Neville-Aitken's algorithm reduces the Lagrange interpolation to the recursive computations of barycenter of two points at each step, with appropriate weights. Thus, for a fixed value of the parameter λ, we get

$$\begin{align*} L[\lambda_0, \lambda_1; y_0, y_1](\lambda) &= (1 - \alpha_{[0,1]}(\lambda)) y_0 + \alpha_{[0,1]}(\lambda) y_1, \\ L[\lambda_1, \lambda_2; y_1, y_2](\lambda) &= (1 - \alpha_{[1,2]}(\lambda)) y_1 + \alpha_{[1,2]}(\lambda) y_2, \\ L[\lambda_0, \lambda_1, \lambda_2; y_0, y_1, y_2](\lambda) &= (1 - \alpha_{[0,2]}(\lambda)) L[\lambda_0, \lambda_1; y_0, y_1](\lambda) + \alpha_{[0,2]}(\lambda) L[\lambda_1, \lambda_2; y_1, y_2](\lambda), \end{align*}$$

with deg(L[λ₀, λ₁; y₀, y₁]) ≤ 1, deg(L[λ₁, λ₂; y₁, y₂]) ≤ 1, deg(L[λ₀, λ₁, λ₂; y₀, y₁, y₂]) ≤ 2, and so on to obtain the Lagrange interpolation polynomial L[λ₀, ..., λₙ; y₀, ..., yₙ], whose degree is less or equal to N.

## 4.2 Extension of the Neville-Aitken's algorithm to Grassmann manifolds

As pointed before, if we are able to define the notion of barycenter of two points on the Grassmann manifold, then we will be able to carry out interpolation on it, via Neville-Aitken's algorithm. More precisely, consider a function f : Λ ⊂ ℝ → Gₘ(ℝⁿ) such that yᵢ := f(λᵢ) are given, for i ∈ {0, 1, ...}. We will describe one way to handle the interpolation of f at { λ₀, λ₁, ..., λₙ}, for any N ≥ 1. At first, consider two (distinct) values λᵢ and λⱼ whose images yᵢ and yⱼ belong to the Grassmann manifold Gₘ(ℝⁿ) and satisfy the condition d(yᵢ, yⱼ) ≤ π/2. Thanks to Theorem 1 there is a unique geodesic γ[yᵢ, yⱼ] : [0, 1] → Gₘ(ℝⁿ) satisfying γ[yᵢ, yⱼ](0) = yᵢ and γ[yᵢ, yⱼ](1) = yⱼ. Then, the map defined on [λᵢ, λⱼ] by

$$\Gamma[\lambda_i, \lambda_j; y_i, y_j] := \gamma[y_i, y_j] \circ \alpha_{[i,j]}, \tag{15}$$

is the unique geodesic on Gₘ(ℝⁿ) satisfying Γ[λᵢ, λⱼ; yᵢ, yⱼ](λᵢ) = yᵢ and Γ[λᵢ, λⱼ; yᵢ, yⱼ](λⱼ) = yⱼ, where α[i,j] is defined in Equation (12). Indeed, geodesics are invariant by affine transformations and we know from Equation (10) a matrix representation of the geodesic connecting two points on the manifold. This matrix representation will be used in the Algorithm 1. In practice to verify that two points X, Y ∈ Gₘ(ℝⁿ) verify d(X, Y) ≤ π/2, we calculate XᵀY = UΣVᵀ and so we have

$$d(\bar{X}, \bar{Y}) = \left( \sum_{i=1}^{m} \theta_i^2 \right)^{1/2}, \tag{16}$$

where Σ = diag(σ₁, ..., σₘ) and θᵢ = arccos(σᵢ) for all 1 ≤ i ≤ m.
**Definition 2.** For a given $\lambda \in [\lambda_i, \lambda_j]$, the point $\Gamma[\lambda_i, \lambda_j; y_i, y_j](\lambda) \in G_m(\mathbb{R}^n)$ is called the barycenter (in the sense of geodesics) on $G_m(\mathbb{R}^n)$ of $y_i$ and $y_j$, weighted by $(1 - \alpha_{[i,j]}(\lambda))$ and $\alpha_{[i,j]}(\lambda)$, respectively.

The barycenter $\Gamma[\lambda_i, \lambda_j; y_i, y_j](\lambda)$ provides an interpolation on $G_m(\mathbb{R}^n)$ of the function $f$ at $\lambda_i$ and $\lambda_j$. Moreover, it coincides with the barycenter in the sense of Karcher.⁵¹ More precisely, we have the following result.

**Proposition 1.** Let $\lambda_i < \lambda_j$ be two parameter values in $\{\lambda_0, \lambda_1, \dots\}$ and $y_i, y_j$ be two points on the Grassmann manifold $G_m(\mathbb{R}^n)$ associated with $\lambda_i$ and $\lambda_j$, respectively. Then

$$ \Gamma[\lambda_i, \lambda_j; y_i, y_j](\lambda) = \arg\min_{y \in G_m(\mathbb{R}^n)} F(y), $$

where

$$ F(y) = \frac{1}{2} \left[ (1 - \alpha_{[i,j]}(\lambda)) \, d^2(y, y_i) + \alpha_{[i,j]}(\lambda) \, d^2(y, y_j) \right], \quad (16) $$

and $d$ denotes the geodesic distance on $G_m(\mathbb{R}^n)$.

*Proof.* From the theorem 4.1 in Reference 41, it is sufficient to show that $\nabla F(\Gamma[\lambda_i, \lambda_j; y_i, y_j](\lambda)) = 0$, or equivalently:

$$ (1 - \alpha_{[i,j]}(\lambda)) \exp_{\Gamma[\lambda_i, \lambda_j; y_i, y_j](\lambda)}^{-1}(y_i) + \alpha_{[i,j]}(\lambda) \exp_{\Gamma[\lambda_i, \lambda_j; y_i, y_j](\lambda)}^{-1}(y_j) = 0. $$

In order to alleviate the notations, we denote $\Gamma[\lambda_i, \lambda_j; y_i, y_j](\lambda) = \hat{y}$ in the sequel. Then, the two geodesic paths $\gamma_1$ and $\gamma_2$ mapping the interval $[0, 1]$ to $G_m(\mathbb{R}^n)$ are

$$ \gamma_1(t) = \Gamma[\lambda_i, \lambda_j; y_i, y_j](\lambda + t(\lambda_i - \lambda)), $$

$$ \gamma_2(t) = \Gamma[\lambda_i, \lambda_j; y_i, y_j](\lambda + t(\lambda_j - \lambda)), $$

from $\hat{y}$ to $y_i$ and from $\hat{y}$ to $y_j$, respectively. Therefore, we get

$$ \gamma'_1(0) = \exp_{\hat{y}}^{-1}(y_i) = (\lambda_i - \lambda)\Gamma[\lambda_i, \lambda_j; y_i, y_j](\lambda), \quad (17) $$

$$ \gamma'_2(0) = \exp_{\hat{y}}^{-1}(y_j) = (\lambda_j - \lambda)\Gamma[\lambda_i, \lambda_j; y_i, y_j](\lambda). \quad (18) $$

Now, combining (17) and (18), we obtain

$$ (1 - \alpha_{[i,j]}(\lambda)) \exp_{\hat{y}}^{-1}(y_i) + \alpha_{[i,j]}(\lambda) \exp_{\hat{y}}^{-1}(y_j) = 0, $$

which ends the proof.

In the same way, we can extend the analogue of (13) to determine the interpolation $\Gamma[\lambda_i, \lambda_j, \lambda_k; y_i, y_j, y_k]$ of $f$ at $\{\lambda_i, \lambda_j, \lambda_k\}$. Indeed, we compute at first $\Gamma[\lambda_i, \lambda_j; y_i, y_j](\lambda)$ and $\Gamma[\lambda_j, \lambda_k; y_j, y_k](\lambda)$ by (14) and then

$$ \Gamma[\lambda_i, \lambda_j, \lambda_k; y_i, y_j, y_k](\lambda) = \Gamma[\lambda_i, \lambda_k; \Gamma[\lambda_i, \lambda_j; y_i, y_j](\lambda), \Gamma[\lambda_j, \lambda_k; y_j, y_k](\lambda)](\lambda), $$

and so on to obtain the interpolation $\Gamma[\lambda_0, \dots, \lambda_N; y_0, \dots, y_N]$ of $f$ at $\{\lambda_0, \dots, \lambda_N\}$, that is:

$$ \Gamma[\lambda_0, \dots, \lambda_N; y_0, \dots, y_N](\lambda) = \Gamma[\lambda_0, \lambda_N; \Gamma[\lambda_0, \dots, \lambda_{N-1}; y_0, \dots, y_{N-1}](\lambda), \Gamma[\lambda_1, \dots, \lambda_N; y_1, \dots, y_N](\lambda)](\lambda). $$

The resulting Neville-Aitken algorithm on Grassmann manifolds is given in Algorithm 1. Once the interpolated basis $\Phi^I(\lambda^*)$ has been computed, the POD-ROM associated with the parameter $\lambda^*$ is built as in the step 2 described in Section 2.

**Remark 1** (Computational complexity). The asymptotic complexity associated with the Grassmannian Neville-Aitken interpolation method (Algorithm 1) is estimated as follows. The calculation of a geodesic connecting two points
**Algorithm 1.** Neville–Aitken’s algorithm on Grassmann manifolds

**Data:**

i. List of $(N + 1)$ parameter values: $\Lambda = \{\lambda_0, \lambda_1, \dots, \lambda_N\} \subset \mathbb{R}$.

ii. List of $(N + 1)$ bases: $\{\Phi(\lambda_0), \Phi(\lambda_1), \dots, \Phi(\lambda_N)\} \subset \mathbb{R}^{n \times m}$.

iii. The new parameter value: $\lambda^* \in \left[ \min_{0 \le k \le N} \lambda_k, \max_{0 \le k \le N} \lambda_k \right]$.

**Result:** The interpolated basis $\Phi^I(\lambda^*)$

for $i = 0$ to $N$ do
$Y(i, 0) \leftarrow \Phi(\lambda_i)$

end for

for $j = 0$ to $N - 1$ do

for $i = 0$ to $N - j - 1$ do
$Y(i, j + 1) \leftarrow \Gamma[\lambda_i, \dots, \lambda_{i+j+1}; y_i, \dots, y_{i+j+1}](\lambda^*)$

end for

end for

$\Phi^I(\lambda^*) \leftarrow Y(0, N)$

$\bar{\Phi}_1, \bar{\Phi}_2$ of the Grassmann manifold $G_m(\mathbb{R}^n)$ involves the calculation of a thin SVD which asymptotic cost is $\mathcal{O}(n m^2)$. So, given $N + 1$ points $\bar{\Phi}_0, \dots, \bar{\Phi}_N$ of the Grassmann manifold associated with the parameters $\{\lambda_i\}_{0 \le i \le N} \subset \mathbb{R}$, the total asymptotic complexity to conduct the Grassmannian Neville–Aitken’s interpolation for a new parameter is $\frac{N(N+1)}{2} \times \mathcal{O}(n m^2)$. We can notice that this total complexity depends linearly on the size $n$ of the high-fidelity model.

*Remark 2* (Multiparametric interpolation). The Grassmannian Neville–Aitken interpolation method (Algorithm 1) have been derived for one dimensional parameters $\lambda \in \mathbb{R}$. It can be extended to tackle multiparametric interpolation with $\lambda \in \mathbb{R}^p$ by interpolating successively in each direction of the parameter space as described in Algorithm 2. However, the resulting interpolated basis is sensitive to the ordering of the scalar components in vector parameter $\lambda$.

5 | APPLICATIONS

The proposed Neville–Aitken’s algorithm on $G_m(\mathbb{R}^n)$ is applied in this section to the interpolation of POD bases for the construction of ROMs associated with three fluid dynamics applications, namely (i) the Von Karman vortex shedding street, (ii) the lid-driven cavity with inflow, and (iii) the flow induced by a rotating solid.

1. We derive a High-Dimensional Model (HDM) by a standard finite-element discretization of the governing equations.

2. We build $(N + 1)$ POD bases by applying the method of snapshots⁵² to the solution obtained from the simulation of the HDM for a set of $(N + 1)$ sampling parameters $\Lambda_N = \{\lambda_0, \dots, \lambda_N\}$.

3. We interpolate the set of POD bases at a given new parameter $\lambda^* \in [\min_{0 \le k \le N} \lambda_k, \max_{0 \le k \le N} \lambda_k]$.

4. We build a ROM (POD-ROM) by a Galerkin projection of the governing equations on the POD bases obtained from interpolation as described in Section 2.

5. The results are compared with the POD-ROM built from the POD basis computed directly from the snapshots of the HDM for $\lambda^*$.

Details on the POD-ROM can be found in References 53-55. All the numerical tests have been performed using the Python/C++ finite element library DOLFIN⁵⁶ on a computer equipped with a processor<fnref n="1" /> Intel Xeon E5-2620 v4 and 64Go of RAM.

<footnote n="1">¹two sockets, eight cores for each socket, two threads for each core, cadenced at 2.10 GHz with a cache of 20 MB.</footnote>
**Algorithm 2.** Multiparametric interpolation over the Grassmann manifolds

**Data:**

i. Set of $p$ lists of $(N_i + 1)_{1 \le i \le p}$ parameter values: $\Lambda_i = \{\hat{\lambda}_{i,0}, \dots, \hat{\lambda}_{i,N_i}\} \subset \mathbb{R}$.

ii. List of $N = \prod_{i=1}^p (N_i + 1)$ bases $\{\Phi(\lambda_1), \dots, \Phi(\lambda_N)\} \subset \mathbb{R}_*^{n \times m}$ associated with the parameters $\lambda_\kappa = (\hat{\lambda}_{1,k_1}, \dots, \hat{\lambda}_{p,k_p})$, with multi-index $\kappa = (k_1, \dots, k_p)$, $0 \le k_i \le N_i$ and $0 \le |\kappa| \le \sum_{i=1}^p N_i$.

iii. The new parameter value: $\lambda^* \in \prod_{i=1}^p \left[ \min_{0 \le k_i \le N_i} \hat{\lambda}_{i,k_i}, \max_{0 \le k_i \le N_i} \hat{\lambda}_{i,k_i} \right]$.

**Result:** The interpolated basis $\Phi^I(\lambda^*)$

for $0 \le |\kappa| \le N$ do

$$ \hat{\Phi}_0(k) \leftarrow \Phi(\lambda_\kappa) $$

end for

for $s = 1$ to $p$ do

for $0 \le |\hat{k}_s| \le \prod_{i=s+1}^p (N_i + 1)$ with $\hat{k}_s = (k_{s+1}, \dots, k_p)$ do

$\hat{\Phi}_s(\hat{k}_s) \leftarrow$ One-dimensional interpolation of $(\hat{\Phi}_{s-1}(0, k_{s+1}, \dots, k_p), \dots, \hat{\Phi}_{s-1}(N_s, k_{s+1}, \dots, k_p))$ w.r.t the $s$-th component of $\lambda^*$

end for

end for

$$ \Phi^I(\lambda^*) \leftarrow \hat{\Phi}_p(\hat{k}_p) $$

## 5.1 | Von Karman vortex shedding street

Here, we consider the planar flow of an incompressible Newtonian fluid around a solid cylinder in a channel. The configuration is shown in Figure 2. The governing incompressible Navier-Stokes equations are described thereafter. The physical domain is $\Omega \subset \mathbb{R}^2$ and the time horizon is $T \subset \mathbb{R}_+$. We search the velocity $v$ and the pressure $p$ in adequate functional spaces$^{57}$ such that

$$ \left\{ \begin{array}{ll} \rho \left( \frac{\partial v}{\partial t} + v \cdot \nabla v \right) = \mu \Delta v - \nabla p + f, & \forall \mathbf{x} \in \Omega, \quad \forall t \in T, \\ \nabla \cdot v = 0, & \forall \mathbf{x} \in \Omega, \quad \forall t \in T, \\ v = v_0, & \forall \mathbf{x} \in \Omega, \quad t = 0, \\ v = v_D, & \forall \mathbf{x} \in \Gamma_D, \quad \forall t \in T, \\ \sigma \cdot n = 0, & \forall \mathbf{x} \in \Gamma_N = \Gamma \setminus \Gamma_D, \quad \forall t \in T, \end{array} \right. \qquad (19) $$

where $f$ denotes body forces density, $\mu$ is the dynamic viscosity, $v_0$ denotes the initial condition, $v_D$ denotes the constant Dirichlet boundary condition for the velocity imposed on $\Gamma_D$, $\sigma = 2\mu D(v) - pI_2$ denotes the Cauchy stress tensor with $D(v) = \frac{1}{2}(\nabla v + \nabla v^T)$, $I_2 \in \mathbb{R}^{2 \times 2}$ the identity tensor, and $\Gamma_N$ denotes the boundary on which outflow boundary condition is imposed.

The dimensionless parameter in this application is the Reynolds number $\lambda = \frac{\rho v_\infty D}{\mu}$ where $\rho = 1 \text{ kg.m}^{-2}$ is the fluid density, $v_\infty = 1 \text{ m.s}^{-1}$ is the magnitude of the velocity at inflow and $D = 1 \text{ m}$ is the cylinder diameter.

The HDM is derived by a standard variational formulation of (19) over the Taylor-Hood finite-elements space (P2 for velocity and P1 for the pressure). The mesh includes 7893 nodes. For each sampling parameter $\lambda \in \Lambda_N$, the simulation is performed by a standard predictor/corrector scheme (see e.g., Reference 58, section 7.3) and a set of $n' = 200$ regularly spaced snapshots of both the velocity $(v(x, t_i, \lambda))_{1 \le i \le n'}$ and the pressure $(p(x, t_i, \lambda))_{1 \le i \le n'}$ are selected over a period $T = 20 \text{ s}$ which does not include the transient period. To ensure homogeneous boundary conditions for the POD modes, the velocity (respectively the pressure) is split into a mean field $\bar{v}(x, \lambda) = \frac{1}{n'} \sum_{i=1}^{n'} v(x, t_i, \lambda)$ (respectively $\bar{p}(x, \lambda) = \frac{1}{n'} \sum_{i=1}^{n'} p(x, t_i, \lambda)$) and a fluctuating field:

$$ \bar{v}(x, t_i, \lambda) = v(x, t_i, \lambda) - \bar{v}(x, \lambda), \qquad (20) $$

$$ \tilde{p}(x, t_i, \lambda) = p(x, t_i, \lambda) - \bar{p}(x, \lambda). \qquad (21) $$
<figure><img src="image_4.png" /><figcaption>FIGURE 2 Description of the configuration for the application of Section 5.1 (Von Karman vortex shedding street). The horizontal (respectively vertical) dimension of the channel is 40 m (respectively 10 m). The solid cylinder is centered at x = (10, 5) and has a diameter D = 1 m. The Dirichlet boundary conditions are inflow v = (1, 0) on Γ<sub>D,1</sub>, slip condition v · n = 0 on top and bottom boundaries Γ<sub>D,2</sub>, and no-slip condition v = (0, 0) on the disk boundary Γ<sub>D,0</sub>. The Neumann boundary condition is outflow σ · n = 0 on Γ<sub>N</sub>, with σ the Cauchy stress tensor and n the outward normal unit vector on the boundary [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

Then, we compute the POD basis $\Phi(x, \lambda) = (\varphi_i(x, \lambda))_{1\le i \le n'}$ (respectively $\Psi(x, \lambda) = (\psi_i(x, \lambda))_{1\le i \le n'}$) for the fluctuating velocity (respectively pressure) associated with each parameter by the method known as snapshots-POD.⁵² The POD bases are truncated to $m$ modes for every parameters $\{\lambda_0, \dots, \lambda_N\}$ ($m=10$ in this section) so that the instantaneous fields are approximated as

$$ \hat{v}(x, t, \lambda) = \bar{v}(x, \lambda) + \sum_{i=1}^{m} a_i(t, \lambda) \varphi_i(x, \lambda), \qquad (22) $$

$$ \hat{p}(x, t, \lambda) = \bar{p}(x, \lambda) + \sum_{i=1}^{m} b_i(t, \lambda) \psi_i(x, \lambda). \qquad (23) $$

Then we construct so-called POD-ROMs by a standard Galerkin projection of the momentum equation (19) onto (i) the POD basis for the velocity and (ii) the gradient of the POD basis for the pressure (see e.g., Reference 59). This yields dynamical systems for the temporal coefficients $(a_i)_{1\le i \le m}$ and $(b_i)_{i \le i \le m}$ in (22) and (23).

First, we randomly select 100 cases to give an overview of the performances of the method and comparison with the method proposed in Reference 37 (Section 5.1.1). Second, we illustrate the sensitivity of the proposed method to the choice for the sampling parameters by a numerical example (Section 5.1.2). Third, we provide a numerical analysis of the order of the method (Section 5.1.3).

The criteria used in the sequel to evaluate the quality of the interpolation are:

Grassmann distance (15) between the subspace spanned by the interpolated POD basis $\Phi^I(\lambda^*)$ and the subspace spanned by the target POD basis $\Phi(\lambda^*)$ computed from the snapshots of the HDM;
Projection error for the fluctuating velocity defined for a given POD basis $\Phi$ as follows:

$$ \varepsilon_{\Phi}^{\text{proj}} \triangleq \frac{\| \tilde{\mathbf{V}}_{\lambda^*} - P_{\Phi} \tilde{\mathbf{V}}_{\lambda^*} \|_{F}^{2}}{\| \tilde{\mathbf{V}}_{\lambda^*} \|_{F}^{2}}, \qquad (24) $$

where the columns of the matrix $\tilde{\mathbf{V}}_{\lambda}$ are the snapshots of the discrete fluctuating velocity field obtained by the simulation of the HDM for the parameter $\lambda$, $P_{\Phi} \tilde{\mathbf{V}} = \sum_{i=1}^{m} (\varphi_i | \tilde{\mathbf{V}})_{L^2(\Omega)}$ $\varphi_i$ denotes the matrix whose columns are the columns of $\tilde{\mathbf{V}}$ projected onto the subspace spanned by the POD basis $\Phi$, and $\| \cdot \|_F$ denotes the Frobenius norm;
for the velocity field reconstructed from the temporal coefficients obtained by the ROM

ROM error

$$ \varepsilon_{\Phi}^{\text{rom}} \triangleq \frac{\| \tilde{\mathbf{V}}_{\lambda^*} - \hat{\mathbf{V}}_{\Phi} \|_{F}^{2}}{\| \tilde{\mathbf{V}}_{\lambda^*} \|_{F}^{2}}, \qquad (25) $$
<table><caption>**TABLE 1** Average over the 100 test cases detailed in Table A1 of (i) the Grassmann distance $d(\bar{\Phi}(\lambda^*), \bar{\Phi}^1(\lambda^*))$ defined in (15), (ii) the projection error $\epsilon_{\Phi^i(\lambda^*)}^{\text{proj}}$ defined in (24) and (iii) the ROM error $\epsilon_{\Phi^i(\lambda^*)}^{\text{rom}}$ defined in (25) with $\Phi^i(\lambda^*)$ the interpolated basis and $\Phi(\lambda^*)$ the reference basis computed directly from the HDM snapshots. The correspondence of the labels for the interpolation methods is given in Section 5.1.1</caption><thead><tr><th></th><th>GN</th><th>TGL</th><th>TGL-R2</th><th>TGL-R4</th></tr></thead><tbody><tr><td>Grassmann distance</td><td>1.730e+00</td><td>2.070e+00</td><td>2.035e+00</td><td>2.121e+00</td></tr><tr><td>Projection error</td><td>1.472e-01</td><td>2.497e-01</td><td>2.427e-01</td><td>2.734e-01</td></tr><tr><td>ROM error</td><td>2.096e-01</td><td>3.237e-01</td><td>3.129e-01</td><td>3.357e-01</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**FIGURE 3** Projection error (24) for the randomly selected cases (Von Karman vortex shedding street, Section 5.1.1) [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

where the matrix $\hat{V}_{\Phi}$ is the reconstruction of the fluctuating velocity $\hat{v} - \bar{v}$ in (22) for the coefficients $(a_i)_{1 \le i \le m}$ obtained by the simulation of the POD-ROM associated with the basis $\Phi$.

### 5.1.1 | Randomly selected cases

We construct 100 cases by randomly selecting $(N + 1) = 5$ sampling parameters $\Lambda_N = \{\lambda_0, \dots, \lambda_N\}$, $\lambda_0 < \dots < \lambda_N$, and a target parameter $\lambda^* \in (\lambda_0, \lambda_N)$ in a set of available parameters $\Lambda_{\text{available}} = \{70, 75, \dots, 95, 100, 110, \dots, 490, 500\}$. These test cases are given in Table A1. We use four methods to interpolate the sets $(\Phi(x, \lambda))_{\lambda \in \Lambda_N}$ and $(\Psi(x, \lambda))_{\lambda \in \Lambda_N}$ at the parameter $\lambda^*$:

* GN the Neville-Aitken's algorithm on $G_m(\mathbb{R}^n)$ proposed in Section 4.2;

* TGL the Lagrange interpolator on $T_{\bar{\Phi}(\lambda_{\text{ref}})} G_m(\mathbb{R}^n)$ proposed in reference³⁷ with the reference point chosen as the point in the sample bases associated with the parameter closest to the target parameter: $\lambda_{\text{ref}} = \lambda_i$ with $i = \operatorname{argmin}_{j \in \{0, \dots, N\}} |\lambda^* - \lambda_j|$;

* TGL-R2 same method as TGL with the reference point chosen as the middle point in the sample set $\lambda_{\text{ref}} = \lambda_2$;

* TGL-R4 same method as TGL with the reference point chosen as the last point in the sample set $\lambda_{\text{ref}} = \lambda_4$.

The averages over the 100 cases for the criteria (Grassmann distance, projection error and ROM error) are given in Table 1, where we see that the proposed Neville-Aitken's interpolation over $G_m(\mathbb{R}^n)$ yields globally better results compared with interpolation over $T_{\bar{\Phi}(\lambda_{\text{ref}})} G_m(\mathbb{R}^n)$ for all reference points $\Phi(\lambda_{\text{ref}})$. The projection errors for every 100 cases for the methods **GN** and **TGL** are shown in Figure 3, and the detailed values for this criterion are given in Table 2 for some selected cases. It should be remarked in Table 2 that the natural choice for $\lambda_{\text{ref}}$ in the method **TGL** as the parameter in the sampling set $\Lambda$ closest to the target $\lambda^*$ does not systematically yield the better results, so that the selection for $\lambda_{\text{ref}}$ is not trivial.
<table><caption>**TABLE 2** Projection error $\epsilon_{\Phi^I(\lambda^*)}^{\text{proj}}$ defined in (24) with $\Phi^I(\lambda^*)$ the interpolated basis in some selected cases for the application of Section 5.1.1 (Von Karman vortex shedding street). The cases are detailed in Table A1. The correspondence of the labels for the interpolation methods is given in Section 5.1.1</caption><thead><tr><th>Case</th><th>GN</th><th>TGL</th><th>TGL-R2</th><th>TGL-R4</th></tr></thead><tbody><tr><td>10</td><td>1.35e-01</td><td>1.89e-01</td><td>1.94e-01</td><td>2.03e-01</td></tr><tr><td>20</td><td>2.71e-01</td><td>3.24e-01</td><td>3.25e-01</td><td>4.52e-01</td></tr><tr><td>30</td><td>1.89e-02</td><td>2.79e-02</td><td>2.79e-02</td><td>3.66e-02</td></tr><tr><td>40</td><td>2.68e-03</td><td>4.60e-03</td><td>4.07e-03</td><td>4.44e-03</td></tr><tr><td>50</td><td>4.58e-02</td><td>2.59e-01</td><td>1.95e-01</td><td>2.58e-01</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**FIGURE 4** Drag effort $F_x$ and lift effort $F_y$ from (26) on the disk boundary for the Von Karman vortex shedding street of Section 5.1.1 for the randomly selected case 10 detailed in Table A1. **HDM**: efforts computed directly from the snapshots of the high-dimensional model; **GN** and **TGL**: efforts computed from the snapshots reconstructed by (22) and (23) for the POD basis interpolated by the proposed Neville–Aitken interpolation over $G_m(\mathbb{R}^n)$ and by a Lagrange interpolation over the tangent space $T_{\bar{\Phi}_{(\lambda_{ref})}} G_m(\mathbb{R}^n)$ (see Section 5.1.1) [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

Finally, we compare the aerodynamic efforts imposed by the fluid on the cylinder, defined as follows:

$$ \left( \frac{F_x}{F_y} \right) = \int_{\Gamma_{D,0}} \sigma(x, t) \cdot n(x) \, d\ell, \qquad (26) $$

where $\Gamma_{D,0}$ denotes the disk boundary, $\sigma$ is the fluid stress tensor, $n$ is the outward unit vector normal to $\Gamma_{D,0}$, and $d\ell$ is a line element. Results are shown in Figure 4 for the case 10 detailed in Table A1, where we see that the efforts obtained by the POD-ROM built from bases interpolated by the proposed Grassmannian Neville–Aitken interpolator are closer to the efforts computed directly from the HDM snapshots compared with the tangent interpolation method **TGL**.

### 5.1.2 | Numerical tests for the influence of the sampling set

In this section, we heuristically investigate the influence of the selection of the sampling set on the interpolation results from numerical tests. To do so, we consider a reference sampling set with $N + 1 = 3$ parameters $\Lambda_{\text{ref}} = \{\lambda_0, \lambda_1, \lambda_2\}$ with associated POD basis $\{\Phi(\lambda_i)\}_{0 \le i \le 2}$. Bases interpolated in this sampling set at the target parameter $\lambda^*$ are denoted $\Phi_{\text{ref}}^I(\lambda^*)$. We also define modified sampling sets $\Lambda_{\epsilon} = \{\lambda_0, \lambda_1 + \epsilon, \lambda_2\}$ obtained by moving the midpoint from $\lambda_1$ to $\lambda_1 + \epsilon$. Bases interpolated in the latter sampling set at the target parameter $\lambda^*$ are denoted $\Phi_{\epsilon}^I(\lambda^*)$. We focus on the following quantity of interest inspired from the Lebesgue constant defined for interpolation of vector spaces but not defined for interpolation on manifolds:

$$ L_{\epsilon}(\lambda^*) = \frac{d(\bar{\Phi}_{\text{ref}}^I(\lambda^*), \bar{\Phi}_{\epsilon}^I(\lambda^*))}{d(\bar{\Phi}(\lambda_1), \bar{\Phi}(\lambda_1 + \epsilon))}. \qquad (27) $$
<figure><img src="image_7.png" /><figcaption>FIGURE 5 Numerical investigation of the influence of the sampling set on interpolation results, with the quantity $L_{\epsilon}(\lambda^{*})$ defined in (27) (see details are given in Section 5.1.2) [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 6 Numerical analysis of the order of interpolation methods over $G_m(\mathbb{R}^n)$ (see Section 5.1.3). **GN**: Neville–Aitken algorithm on $G_m(\mathbb{R}^n)$; **TGL**: Lagrange interpolation over $T_{\bar{\Phi}(\lambda_{\text{ref}})}G_m(\mathbb{R}^n)$, with $\lambda_{\text{ref}}$ the parameter closest to the target $\lambda^*$ in the sampling set. The Grassmann distance between the target and the interpolated bases is shown relatively to the maximum distance obtained for the method **GN**, and $h_{\text{max}} = 10^{-4}$ [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

For the numerical tests, we choose the reference parameters $\Lambda_{\text{ref}} = \{250, 300, 350\}$ and the modified parameters $\Lambda_{-5}$, $\Lambda_{-1}$, $\Lambda_{+1}$, $\Lambda_{+5}$, and $L_{\epsilon}(\lambda^*)$ is computed for $\lambda^* \in \{250, 260, \dots, 290, 295, 299, 300, 301, 305, 310, \dots, 350\}$. We see in Figure 5 that the variations of $L_{\epsilon}(\lambda^*)$ for every $\epsilon$ is small and bounded so that, in this example, the proposed method is not highly sensitive to the sampling set. Note that this test may not be the most pertinent to study a kind of generalization of the notion of Lebesgue constant and that it is provided only for illustration purpose.

### 5.1.3 | Numerical tests for the method order

In this section, we examine numerically the order of the proposed method **GN** and compare with the method **TGL** (see labels definitions in Section 5.1.1). For this test, we use $N + 1 = 2$ points $\Lambda_h = \{\lambda^* - h, \lambda^* + h\}$ with the target $\lambda^* = 150$, and we compute the Grassmann distance $d(\bar{\Phi}(\lambda^*), \bar{\Phi}_h^1(\lambda^*))$ for $h \in \{10^{-4}, 10^{-5}, 10^{-6}, 10^{-7}\}$. The results are given in Figure 6, where we see that the order of the proposed method with respect to the Grassmann distance is $N + 1$ in this case, which coincides with the order of the method proposed in Reference 37.

## 5.2 | Multiparametric interpolation: the lid-driven cavity with inflow

Here, we consider the planar flow of an incompressible Newtonian fluid in a lid-driven cavity with an inflow and an outflow (see description in Figure 7). The governing equations are similar to the previous case of Section 5.1. We define
<figure><img src="image_9.png" /><figcaption>FIGURE 7 Description of the configuration for the application of Section 5.2 (lid-driven cavity with inflow). The dimensions of the cavity are 1 m × 1 m and the inflow/outflow extends by 0.1 m. The Dirichlet boundary condition for the velocity is $v|_{\Gamma_{D,1}} = (1, 0)$ and $v|_{\Gamma_{D,2}} = (v_{\text{inflow}}, 0)$ with $v_{\text{inflow}} \in (0, 1)$ the second parameter. The outflow boundary condition is $\sigma \cdot n = 0$ on $\Gamma_N$ (with $\sigma$ the Cauchy stress tensor and $n$ the outward normal unit vector on the boundary). Homogeneous Dirichlet boundary conditions are imposed on the remaining boundary $\Gamma_{D,0} = \partial\Omega \setminus (\Gamma_{D,1} \cup \Gamma_{D,2} \cup \Gamma_N)$. The snapshot in (B) represents the magnitude of the velocity at time $t = 5$ s for the parameter $\lambda = (10^{-3}, 0.6)$ as defined in (28). (A) Schematic; (B) Snapshot [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

<table><caption>TABLE 3 Description of the cases for the application of Section 5.2 (lid-driven cavity with inflow). The vector parameter $\lambda$ is defined in (28). The set of sampling parameters for which the simulation of the high-dimensional model is performed is $\Lambda_{\mu}$ for the first component of $\lambda$ and $\Lambda_{v_{\text{inflow}}}$ for the second component. The interpolation is performed for the new parameter $\lambda^*$</caption><thead><tr><td>Case</td><td>Sampling Λ<sub>μ</sub></td><td>Sampling Λ<sub>v<sub>inflow</sub></sub></td><td>Target λ*</td></tr></thead><tbody><tr><td>1</td><td>(0.001, 0.003, 0.004)</td><td>(0.1, 0.5, 0.6)</td><td>(0.002, 0.3)</td></tr><tr><td>2</td><td>(0.001, 0.002, 0.004)</td><td>(0.1, 0.5, 0.6)</td><td>(0.003, 0.3)</td></tr><tr><td>3</td><td>(0.001, 0.003, 0.004)</td><td>(0.1, 0.5, 0.6)</td><td>(0.002, 0.4)</td></tr></tbody></table>

two Reynolds numbers associated with this configuration. The first is associated with the lid velocity $R_1 = \frac{\rho v_{\text{lid}} \ell}{\mu}$ where $\rho = 1\text{kg.m}^{-3}$ is the fluid density, $v_{\text{lid}} = 1\text{m.s}^{-1}$ is the magnitude of the lid velocity $v|_{\Gamma_{D,1}} = (1, 0)$, $\ell = 1\text{m}$ is the length of the cavity's sides and $\mu$ is the dynamic viscosity so that $R_1 = \mu^{-1}$. The second is the Reynolds number associated with the velocity at inflow: $R_2 = \frac{\rho v_{\text{inflow}} \ell}{\mu}$ with $v_{\text{inflow}}$ the magnitude of inflow velocity $v|_{\Gamma_{D,2}} = (v_{\text{inflow}}, 0)$. Thus, the dimensionless vector parameter selected for this application is

$$ \lambda = \left( \frac{1}{R_1}, \frac{R_2}{R_1} \right) = (\mu, v_{\text{inflow}}). \qquad (28) $$

The HDM is derived exactly as in the previous application (Section 5.1). The mesh includes here 13,656 nodes.

The test cases are given in Table 3. The test procedure is identical to that in previous Section 5.1, with the difference that the snapshots are taken from the initial condition over a transient period of 10 s. We compute sets of POD bases for the fluctuating velocity and pressure, and we derive interpolated bases with the proposed method and that proposed in Reference 37 as the one-dimensional interpolation method in step 5 of Algorithm 2. Every bases are truncated to $m = 20$ modes. The average of the criteria (Grassmann distance, projection error and ROM error) are given in Table 4 where the results obtained by the proposed interpolation over $G_m(\mathbb{R}^n)$ are globally slightly better than the results obtained by interpolation over $T_{\bar{\Phi}(\lambda_{\text{ref}})} G_m(\mathbb{R}^n)$. The details for the projection error and the ROM error are given in Tables 5 and 6, respectively. The simulation of these reduce-order dynamical systems yields the temporal coefficients $(a_i(t))_{1 \le i \le m}$ and $(b_i(t))_{1 \le i \le m}$ in (22) and (23). The coefficients for the velocity $(a_i(t))_{1 \le i \le 4}$ are shown in Figure 8. We see the coefficients obtained by the proposed method GN are closer to the coefficients obtained by projection of the HDM solution when compared with those obtained by the method TGL.
<table><thead><tr><th></th><th>GN</th><th>TGL</th></tr></thead><tbody><tr><td>Grassmann distance</td><td>1.149e+00</td><td>1.378e+00</td></tr><tr><td>Projection error</td><td>5.195e-02</td><td>7.629e-02</td></tr><tr><td>ROM error</td><td>6.718e-02</td><td>1.002e-01</td></tr></tbody></table>

<table><thead><tr><th>Case</th><th>GN</th><th>TGL</th></tr></thead><tbody><tr><td>1</td><td>5.605e-02</td><td>1.146e-01</td></tr><tr><td>2</td><td>8.071e-02</td><td>8.950e-02</td></tr><tr><td>3</td><td>1.908e-02</td><td>2.468e-02</td></tr></tbody></table>

<table><caption>TABLE 4 Average over the three test cases detailed in Table 3 of (i) the Grassmann distance $d(\bar{\Phi}(\lambda^*), \bar{\Phi}^I(\lambda^*))$ defined in (15), (ii) the projection error $\epsilon_{\Phi^I(\lambda^*)}^{\text{proj}}$ defined in (24) and (iii) the ROM error $\epsilon_{\Phi^I(\lambda^*)}^{\text{rom}}$ defined in (25) with $\Phi^I(\lambda^*)$ the interpolated basis and $\Phi(\lambda^*)$ the reference basis computed directly from the HDM snapshots. The correspondence of the labels for the interpolation methods is given in Section 5.1.1</caption><thead><tr><th>Case</th><th>GN</th><th>TGL</th></tr></thead><tbody><tr><td>1</td><td>6.992e-02</td><td>1.497e-01</td></tr><tr><td>2</td><td>1.106e-01</td><td>1.236e-01</td></tr><tr><td>3</td><td>2.090e-02</td><td>2.775e-02</td></tr></tbody></table>

<caption>TABLE 5 Details of the projection errors for the application of Section 5.2 for the test cases detailed in Table 3. The correspondence of the labels for the interpolation methods is given in Section 5.1.1</caption>

<figure><img src="image_10.png" /><figcaption>FIGURE 8 Temporal coefficients $(a_i(t))_{1 \le i \le 4}$ in (22) returned by the dynamical system build from the proper orthogonal decomposition basis obtained directly from the snapshots (HDM) and from the two interpolation methods (GN and TGL) for the application of Section 5.2 (lid-driven cavity with inflow) in the case 1 (see Table 3) [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

## 5.3 Flows induced by a rotating body

Here, we consider a two dimensional circular spatial domain $\Omega = \Omega_S \cup \Omega_F$ filled with a rotating ellipsoidal body $\Omega_S$ immersed in an incompressible Newtonian fluid $\Omega_F$ (see description in Figure 9 and details in Reference 13). The parameter is the Reynolds number defined as

$$ \lambda = \frac{\rho v_{\infty} \ell}{\mu}, \qquad (29) $$
**FIGURE 9** Description of the configuration for the application of Section 5.3 (flow induced by a rotating body). The domain radius is 1 m and the ellipse principal radius is $R=0.8$ m with an aspect ratio of 0.2. The boundary condition is outflow $\sigma \cdot n = 0$ on $\Gamma_N$ (with $\sigma$ the Cauchy stress tensor and $n$ the outward normal unit vector on the boundary). (A) Schematic; (B) Mesh [Colour figure can be viewed at wileyonlinelibrary.com]

<figure><img src="image_11.png" /><figcaption>**TABLE 7** Description of the cases for the application of Section 5.3 (flow induced by a rotating body). The parameter is the Reynolds number as defined in (29). The set of sampling parameters for which the simulation of the full order model is performed is $\Lambda_N$. The interpolation is performed for the new parameter $\lambda^*$</figcaption></figure>

<table><caption>**TABLE 8** Results for the projection error as defined in (24) for the application of Section 5.3 (flow induced by a rotating body). **HDM**: error associated with the proper orthogonal decomposition (POD) basis computed directly from the High-Dimensional Model (HDM) snapshots; **GN** and **TGL**: error associated with the POD basis interpolated by the proposed Neville–Aitken interpolation over $G_m(\mathbb{R}^n)$ and by a Lagrange interpolation over $T_{\Phi(\lambda_{ref})} G_m(\mathbb{R}^n)$ (see labels defined in Section 5.1.1)</caption><thead><tr><td>Case</td><td>Sampling Λ<sub>N</sub></td><td>Target λ*</td></tr></thead><tbody><tr><td>1</td><td>(500, 1000, 1500, 2500, 3000)</td><td>2800</td></tr><tr><td>2</td><td>(500, 1000, 1500, 1800, 2500, 2800, 3000)</td><td>800</td></tr><tr><td>3</td><td>(500, 1500, 3000)</td><td>2800</td></tr></tbody></table>

<table><thead><tr><td>case</td><td>HDM</td><td>GN</td><td>TGL</td></tr></thead><tbody><tr><td>1</td><td>6.493e-04</td><td>1.874e-03</td><td>1.880e-03</td></tr><tr><td>2</td><td>2.747e-04</td><td>2.457e-03</td><td>2.457e-03</td></tr><tr><td>3</td><td>6.493e-04</td><td>2.141e-03</td><td>2.143e-03</td></tr></tbody></table>

with the density $\rho = 1 \text{ kg.m}^{-3}$, the dynamic viscosity $\mu = 0.01 \text{ kg.m}^{-1}\text{s}^{-1}$, the velocity at the ellipse tips $v_{\infty} = R \dot{\theta} (\text{m.s}^{-1})$ with $\dot{\theta}$ the angular velocity ($\text{s}^{-1}$), and $\ell = 2R$ (m) the ellipse principal diameter. The governing equations are derived by extending the incompressible Navier–Stokes equations to the solid domain by the *fictitious domains* method introduced in References 60,61 and developed in Reference 62.

Then, the momentum equation and the continuity equation are solved together by a monolithic formulation on the mixed finite element space known as the *mini space* (i.e., linear vector Lagrange element enriched with cubic vector bubble element for velocity and piecewise linear element for pressure, see Reference 63 for details). The finite-elements mesh includes 52,669 nodes and is not conforming with the body’s boundaries (see Figure 9(B)). The test cases (i.e., the set of sampling parameters and the target parameter) are given in Table 7. The test procedure is identical to that described in Section 5.1, except for the construction of the POD-ROM which has to comply with the multiphase description of the fluid-structure interaction (the interested reader is referred to Reference 13). As in previous Section 5.2, we are concerned here with the *transient period* that coincide with an half rotation of the ellipse. We compute sets of POD bases for the fluctuating velocity, and we derive interpolated bases by the proposed method and that proposed in Reference 37. Every bases are truncated here to $m = 40$ modes. The projection error for the fluctuating velocity defined in (24) are given in Table 8. The error associated with the reduced-order dynamical systems (in Table 9). For this configuration, the proposed method gives results similar to that proposed in Reference 37. Additionally, the temporal coefficients $(a_i(t))_{i=1}^4$ for the fluctuating velocity are shown in Figure 10.

6 | CONCLUSION

In this work, we exploited the recursive character of the well-known Neville–Aitken’s interpolation to perform interpolation over the Grassmann manifold using the geodesic barycenter of two points. This method has been used to interpolate
<table><caption>TABLE 9 Reduced order model error defined in (25) for the application of Section 5.3 (flow induced by a rotating body). **HDM**: error associated with the proper orthogonal decomposition (POD) basis computed directly from the high-dimensional model (HDM) snapshots; **GN** and **TGL**: error associated with the POD basis interpolated by the proposed Neville–Aitken interpolation over $G_m(\mathbb{R}^n)$ and by a Lagrange interpolation over $T_{\Phi(\lambda_{ref})} G_m(\mathbb{R}^n)$ (see labels defined in Section 5.1.1)</caption><thead><tr><th>case</th><th>HDM</th><th>GN</th><th>TGL</th></tr></thead><tbody><tr><td>1</td><td>2.964e-02</td><td>2.108e-02</td><td>2.189e-02</td></tr><tr><td>2</td><td>3.029e-03</td><td>1.817e-02</td><td>1.818e-02</td></tr><tr><td>3</td><td>2.964e-02</td><td>1.474e-02</td><td>1.475e-02</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>FIGURE 10 Temporal coefficients $(a_i(t))_{i=1}^4$ in (22) returned by the dynamical system build from the reference proper orthogonal decomposition basis obtained directly from the snapshots (**HDM**) and from the two interpolation methods (**GN** and **TGL**) for the application of Section 5.3 (flow induced by a rotating body) in the case 2 (see Table 7) [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

POD basis for the construction of POD-ROMs associated with three CFD applications. The benefits of the proposed method are (i) that it does not depend on the nontrivial choice of a reference point on $G_m(\mathbb{R}^n)$ (as shown in Section 5.1.1), and (ii) that it is globally more accurate than the method proposed in Reference 37, with however comparable results in some cases. The asymptotic complexity is linear with respect to the dimension of the underlying vector space, and thus proves comparable to that of the method proposed in Reference 37 for a reasonable number of interpolation points (which can always be achieved by some sort of clustering in the parameter space). Numerical tests show that the proposed method is of the same order than the method in Reference 37. Moreover, it is direct and does not require the resolution of a fixed point problem associated with a minimization process as in Reference 41 in the framework of the IDW method on Grassmann manifolds.

Multiparametric interpolation can be performed with the proposed method, provided a structured sampling of the multi-dimensional parameter space. However, the asymptotic complexity of the interpolation grows with the number of dimensions; this is a limitation of our approach and a direct perspective to this work is to develop sparsity-based method for such cases.

A more general perspective is the development of a general theoretical framework for the interpolation of functions with values in Riemannian manifolds, in order to handle clearly some notions such as stability and interpolation error which are still heuristic.

## ACKNOWLEDGEMENTS

This work was partially supported by the French funding project Projet CPER 2015-2020.

## DATA AVAILABILITY STATEMENT

Data sharing is not applicable to this article as no new data were created or analyzed in this study.
ORCID

Rolando Mosquera ID https://orcid.org/0000-0002-5890-9829

REFERENCES

1. Schilders WH, Van der Vorst HA, Rommes J. *Model Order Reduction: Theory, Research Aspects and Applications*. Vol 13. Berlin, Germany: Springer; 2008.

2. Quarteroni A, Rozza G. *Reduced Order Methods for Modeling and Computational Reduction*. Vol 9. Berlin, Germany: Springer; 2014.

3. Guibert R, Mcleod K, Caiazzo A, et al. Group-wise construction of reduced models for understanding and characterization of pulmonary blood flows from medical images. *Med Image Anal.* 2014;18(1):63-82.

4. Benner P, Gugercin S, Willcox K. A survey of projection-based model reduction methods for parametric dynamical systems. *SIAM Rev.* 2015;57(4):483-531.

5. Benner P. *Model Order Reduction, Volume 3: Applications*. Berlin, Germany: Walter de Gruyter GmbH & Co KG; 2021.

6. Chinesta F, Huerta A, Rozza G, Willcox K. Model order reduction. *Encyclopedia of Computational Mechanics*. New York, NY: Wiley; 2016.

7. Lassila T, Manzoni A, Quarteroni A, Rozza G. Model order reduction in fluid dynamics: challenges and perspectives. *Reduced Order Methods for modeling and computational reduction*. Berlin, Germany: Springer; 2014:235-273.

8. Chinesta F, Keunings R, Leygue A. *The Proper Generalized Decomposition for Advanced Numerical Simulations: A Primer*. Berlin, Germany: Springer Science & Business Media; 2013.

9. Lestandi L, Bhaumik S, Sengupta TK, Avatar GKC, Azaíez MPOD. applied to numerical study of unsteady flow inside lid-driven cavity. *J Math Study.* 2018;51(2):150-176.

10. Dowell EH, Hall KC. Modeling of fluid-structure interaction. *Annu Rev Fluid Mech.* 2001;33(1):445-490.

11. Lieu T, Farhat C, Lesoinne M. Reduced-order fluid/structure modeling of a complete aircraft configuration. *Comput Methods Appl Mech Eng.* 2006;195(41-43):5730-5742.

12. Liberge E, Hamdouni A. Reduced order modelling method via proper orthogonal decomposition (POD) for flow around an oscillating cylinder. *J Fluids Struct.* 2010;26(2):292-311.

13. Falaize A, Liberge E, Hamdouni A. POD-based reduced order model for flows induced by rigid bodies in forced rotation. *J Fluids Struct.* 2019;91:102593.

14. Niroomandi S, Alfaro I, Cueto E, Chinesta F. Model order reduction for hyperelastic materials. *Int J Numer Methods Eng.* 2010;81(9):1180-1206.

15. Buljak V. *Inverse Analyses with Model Reduction: Proper Orthogonal Decomposition in Structural Mechanics*. Berlin, Germany: Springer Science & Business Media; 2011.

16. Zahr MJ, Avery P, Farhat C. A multilevel projection-based model order reduction framework for nonlinear dynamic multiscale problems in structural and solid mechanics. *Int J Numer Methods Eng.* 2017;112(8):855-881.

17. Drohmann M, Haasdonk B, Ohlberger M. Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation. *SIAM J Sci Comput.* 2012;34(2):A937-A969.

18. Thomas JP, Dowell EH, Hall KC. Three-dimensional transonic aeroelasticity using proper orthogonal decomposition-based reduced-order models. *J Aircr.* 2003;40(3):544-551.

19. Carlberg K, Farhat C, Cortial J, Amsallem D. The GNAT method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows. *J Comput Phys.* 2013;242:623-647.

20. Chinesta F, Ladeveze P, Cueto E. A short review on model order reduction based on proper generalized decomposition. *Arch Comput Methods Eng.* 2011;18(4):395-404.

21. Chinesta F, Ammar A, Cueto E. Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models. *Arch Comput Methods Eng.* 2010;17(4):327-350.

22. Ryckelynck D. A priori hyperreduction method: an adaptive approach. *J Comput Phys.* 2005;202(1):346-366.

23. Allery C, Hamdouni A, Ryckelynck D, Verdon N. A priori reduction method for solving the two-dimensional burgers' equations. *Appl Math Comput.* 2011;217(15):6671-6679.

24. Maday Y, Ronquist EM. A reduced-basis element method. *J Sci Comput.* 2002;17(1):447-459.

25. Quarteroni A, Manzoni A, Negri F. *Reduced Basis Methods for Partial Differential Equations: An Introduction*. Vol 92. Berlin, Germany: Springer; 2015.

26. Dumon A, Allery C, Ammar A. Proper generalized decomposition method for incompressible flows in stream-vorticity formulation. *Europ J Comput Mech.* 2010;19(5-7):591-617. https://doi.org/10.3166/ejcm.19.591-617.

27. Dumon A, Allery C, Ammar A. Proper general decomposition (PGD) for the resolution of Navier-Stokes equations. *J Comput Phys.* 2011;230(4):1387-1407. https://doi.org/10.1016/j.jcp.2010.11.010.

28. Dumon A, Allery C, Ammar A. Simulation of heat and mass transport in a square lid-driven cavity with Proper Generalized Decomposition (PGD). *Numer Heat Transf Part B - Fund.* 2013;63(1):18-43. https://doi.org/10.1080/10407790.2012.724991.

29. Dumon A, Allery C, Ammar A. Proper generalized decomposition method for incompressible Navier-Stokes equations with a spectral discretization. *Appl Math Comput.* 2013;219(15):8145-8162.
30. Falcó A, Hackbusch W. On minimal subspaces in tensor representations. *Found Comput Math*. 2012;12(6):765-803.

31. Kolmogoroff A. Uber die beste Annaherung von Funktionen einer gegebenen Funktionenklasse. *Annals of Mathematics*. 1936;37:107-110.
https://doi.org/10.2307/1968691.

32. Buffa A, Maday Y, Patera AT, Prud'homme C, Turinici G. A priori convergence of the greedy algorithm for the parametrized reduced basis method. *ESAIM Math Modell Numer Anal*. 2012;46(3):595-603.

33. Pinkus A. *N-widths in Approximation Theory*. Berlin/Heidelberg, Germany: Springer-Verlag; 1985.

34. Barrault M, Maday Y, Nguyen NC, Patera AT. An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations. *Comptes Rendus Mathematique*. 2004;339(9):667-672.

35. Chaturantabut S, Sorensen DC. Nonlinear model reduction via discrete empirical interpolation. *SIAM J Sci Comput*. 2010;32(5):2737-2764.

36. Sipp D, de Pando MF, Schmid PJ. Nonlinear model reduction: a comparison between POD-Galerkin and POD-DEIM methods. *Comput Fluids*. 2020;208:104628.

37. Amsallem D, Farhat C. An interpolation method for adapting reduced order models and application to aeroelasticity. *Amer Inst Aeronaut Astronaut*. 2008;46:1803-1813.

38. Amsallem D, Cortial J, Farhat C. Towards real-time computational-fluid-dynamics-based aeroelastic computations using a database of reduced-order information. *AIAA J*. 2010;48(9):2029-2037.

39. Amsallem D, Cortial J, Carlberg K, Farhat C. A method for interpolating on manifolds structural dynamics reduced-order models. *Int J Numer Methods Eng*. 2009;80(9):1241-1258.

40. Degroote J, Vierendeels J, Willcox K. Interpolation among reduced-order matrices to obtain parameterized models for design, optimization and probabilistic analysis. *Int J Numer Methods Fluids*. 2010;63(2):207-230.

41. Mosquera R, Hamdouni A, El Hamidi A, Allery C. POD basis interpolation via inverse distance weighting on Grassmann manifolds. *Discrete Cont Dyn Syst Ser S*. 2019;12:1743-1759.

42. Gousenbourger PY, Samir C, Absil PA. Piecewise-Bézier C1 interpolation on Riemannian manifolds with application to 2D shape morphing. Paper presented at: Proceedings of the 2014 22nd International Conference on Pattern Recognition; 2014:4086-4091; IEEE, Tockholm, Sweden.

43. Arnould A, Gousenbourger PY, Samir C, Absil PA, Canis M. Fitting smooth paths on Riemannian manifolds: endometrial surface reconstruction and preoperative MRI-based navigation. Paper presented at: Proceedings of the International Conference on Geometric Science of Information; 2015:491-498; Springer, New York, NY.

44. Press WH, William H, Teukolsky SA, Saul A, Vetterling WT, Flannery BP. *Numerical Recipes: The Art of Scientific Computing*. 3rd ed. Cambridge, MA: Cambridge University Press; 2007.

45. Neville EH. Iterative interpolation. *J Indian Math Soc*. 1934;20:87-120.

46. Milnor W, Stasheff JD. *Characteristic Classes*. Vol 76. Annals of Mathematics Studies, Princeton, NJ: Princeton University Press; 1974.

47. Edelman A, Arias TA, Smith ST. The geometry of algorithms with orthogonality constraints. *SIAM J Matrix Anal Appl*. 1998;20(2):303-353.

48. Absil PA, Mahony R, Sepulchre R. Riemannian geometry of Grassmann manifolds with a view on algorithmic computation. *Phys Rep Rev Sec Phys Lett*. 2004;80:199-220.

49. Kozlov SE. Geometry of real Grassmannian manifolds. *Zap Nauchn Semin POMI*. 1997;246:108-129.

50. Aitken AC. On interpolation by iteration of proportional parts, without the use of differences. *Proc Edinburgh Math Soc Ser*. 1932;2(3):56-76.

51. Karcher H. Riemannian center of mass and mollifier smoothing. *Commun Pure Appl Math*. 1977;30:509-541.

52. Sirovich L. Turbulence and the dynamics of coherent structures, Parts I-III. *Q Appl Math*. 1987;45:561-590.

53. Holmes PJ, Lumley J, Berkooz J, Mattingly J, Wittenberg R. Low dimensional models of coherent structures in turbulence. *Phys Rev Section Phys Lett*. 1997;4:338-384.

54. Cordier L, Bergmann M. *Proper Orthogonal Decomposition: An Overview*. Lecture Series on Post-Processing of Experimental and Numerical Data. Vol 4. Sint-Genesius-Rode, Belgium: Von Karman Institute for Fluid Dynamics; 2003:1-46.

55. Volkwein S. *Proper Orthogonal Decomposition: Theory and Reduced-Order Modeling*. Konstanz, Germany: University of Konstanz; 2013.

56. Logg A, Wells GN, Hake J. DOLFIN: a C++/Python finite element library. *Automated Solution of Differential Equations by the Finite Element Method*. New York, NY: Springer; 2012:173-225.

57. Boyer F, Fabrie P. *Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models*. Vol 183. New York, NY: Springer; 2013.

58. Ferziger JH, Peric M. *Computational Methods for Fluid Dynamics*. Berlin, Germany: Springer Science & Business Media; 2012.

59. Tallet A, Allery C, Leblond C, Liberge E. A minimum residual projection to build coupled velocity-pressure POD-ROM for incompressible Navier-Stokes equations. *Comm Nonlinear Sci Numer Simulat*. 2015;22:909-932.

60. Glowinski R, Pan TW, Hesla TI, Joseph DD, Periaux J. A fictitious domain method with distributed Lagrange multipliers for the numerical simulation of particulate flow. *Contemp Math*. 1998;218:121-137.

61. Glowinski R, Pan TW, Hesla TI, Joseph DD. A distributed Lagrange multiplier/fictitious domain method for particulate flows. *Int J Multiphase Flow*. 1999;25(5):755-794.
62. Patankar NA, Singh P, Joseph DD, Glowinski R, Pan TW. A new formulation of the distributed Lagrange multiplier/fictitious domain method for particulate flows. *Int J Multiphase Flow*. 2000;26(9):1509-1524.

63. Arnold D, Brezzi F, Fortin M. A stable finite element for the stokes equations. *Calcolo*. 1984;21:337-344.

**How to cite this article:** Mosquera R, El Hamidi A, Hamdouni A, Falaize A. Generalization of the Neville-Aitken interpolation algorithm on Grassmann manifolds: Applications to reduced order model. *Int J Numer Meth Fluids*. 2021;93:2421-2442. https://doi.org/10.1002/fld.4981

APPENDIX. RANDOM CASES FOR THE APPLICATION OF SECTION 4.1

<table><caption>TABLE A1 Random sampling sets $\Lambda_N$ and target parameters $\lambda^*$ for the application of Section 5.1 (Von Karman vortex shedding street) where the parameter $\lambda$ corresponds to the Reynolds number</caption><thead><tr><th>Case</th><th>Sampling Λ<sub>N</sub></th><th>Target λ*</th><th>Case</th><th>Sampling Λ<sub>N</sub></th><th>Target λ*</th></tr></thead><tbody><tr><td>1</td><td>[90, 170, 320, 350, 460]</td><td>200</td><td>51</td><td>[120, 200, 330, 380, 440]</td><td>210</td></tr><tr><td>2</td><td>[85, 190, 370, 470, 500]</td><td>150</td><td>52</td><td>[85, 190, 360, 400, 480]</td><td>250</td></tr><tr><td>3</td><td>[130, 290, 420, 460, 480]</td><td>340</td><td>53</td><td>[150, 200, 400, 470, 480]</td><td>360</td></tr><tr><td>4</td><td>[70, 410, 460, 480, 500]</td><td>390</td><td>54</td><td>[140, 160, 170, 320, 510]</td><td>250</td></tr><tr><td>5</td><td>[70, 80, 160, 400, 420]</td><td>75</td><td>55</td><td>[90, 100, 120, 270, 410]</td><td>150</td></tr><tr><td>6</td><td>[75, 90, 230, 350, 370]</td><td>85</td><td>56</td><td>[110, 140, 270, 360, 420]</td><td>160</td></tr><tr><td>7</td><td>[100, 260, 340, 400, 440]</td><td>250</td><td>57</td><td>[150, 300, 330, 400, 410]</td><td>170</td></tr><tr><td>8</td><td>[85, 160, 200, 380, 510]</td><td>180</td><td>58</td><td>[140, 220, 240, 390, 490]</td><td>180</td></tr><tr><td>9</td><td>[95, 130, 150, 470, 490]</td><td>170</td><td>59</td><td>[70, 80, 90, 400, 480]</td><td>190</td></tr><tr><td>10</td><td>[120, 220, 330, 380, 510]</td><td>140</td><td>60</td><td>[120, 340, 400, 500, 510]</td><td>440</td></tr><tr><td>11</td><td>[85, 100, 350, 400, 440]</td><td>150</td><td>61</td><td>[150, 250, 320, 390, 480]</td><td>210</td></tr><tr><td>12</td><td>[70, 260, 280, 290, 410]</td><td>250</td><td>62</td><td>[180, 200, 270, 330, 500]</td><td>210</td></tr><tr><td>13</td><td>[160, 180, 310, 330, 350]</td><td>200</td><td>63</td><td>[70, 90, 320, 450, 510]</td><td>180</td></tr><tr><td>14</td><td>[90, 110, 350, 380, 480]</td><td>220</td><td>64</td><td>[90, 120, 280, 310, 350]</td><td>270</td></tr><tr><td>15</td><td>[75, 90, 310, 400, 510]</td><td>190</td><td>65</td><td>[70, 90, 140, 170, 390]</td><td>85</td></tr><tr><td>16</td><td>[150, 290, 360, 400, 440]</td><td>190</td><td>66</td><td>[95, 110, 140, 420, 500]</td><td>280</td></tr><tr><td>17</td><td>[190, 300, 350, 450, 460]</td><td>240</td><td>67</td><td>[80, 270, 360, 390, 490]</td><td>160</td></tr><tr><td>18</td><td>[90, 130, 220, 410, 490]</td><td>170</td><td>68</td><td>[80, 85, 130, 190, 260]</td><td>180</td></tr><tr><td>19</td><td>[95, 210, 300, 340, 390]</td><td>240</td><td>69</td><td>[160, 270, 370, 470, 480]</td><td>340</td></tr><tr><td>20</td><td>[80, 240, 250, 340, 490]</td><td>170</td><td>70</td><td>[70, 240, 290, 430, 480]</td><td>85</td></tr><tr><td>21</td><td>[140, 180, 330, 410, 490]</td><td>310</td><td>71</td><td>[150, 270, 280, 440, 490]</td><td>350</td></tr><tr><td>22</td><td>[80, 85, 120, 360, 450]</td><td>90</td><td>72</td><td>[75, 85, 150, 410, 510]</td><td>380</td></tr><tr><td>23</td><td>[100, 180, 330, 400, 450]</td><td>380</td><td>73</td><td>[75, 80, 340, 460, 490]</td><td>290</td></tr><tr><td>24</td><td>[70, 190, 260, 290, 310]</td><td>170</td><td>74</td><td>[85, 250, 330, 410, 480]</td><td>130</td></tr><tr><td>25</td><td>[70, 75, 240, 430, 490]</td><td>380</td><td>75</td><td>[70, 190, 360, 390, 430]</td><td>300</td></tr><tr><td>26</td><td>[80, 340, 350, 420, 500]</td><td>150</td><td>76</td><td>[70, 220, 230, 440, 470]</td><td>95</td></tr><tr><td>27</td><td>[95, 280, 360, 390, 480]</td><td>370</td><td>77</td><td>[200, 390, 410, 440, 480]</td><td>320</td></tr><tr><td>28</td><td>[200, 240, 300, 340, 480]</td><td>320</td><td>78</td><td>[70, 75, 390, 400, 460]</td><td>280</td></tr></tbody></table>
<table><caption>TABLE A1 (Continued)</caption><thead><tr><th>Case</th><th>Sampling Λ<sub>N</sub></th><th>Target λ*</th><th>Case</th><th>Sampling Λ<sub>N</sub></th><th>Target λ*</th></tr></thead><tbody><tr><td>29</td><td>[100, 120, 160, 370, 420]</td><td>240</td><td>79</td><td>[240, 290, 430, 500, 510]</td><td>460</td></tr><tr><td>30</td><td>[90, 140, 270, 390, 460]</td><td>300</td><td>80</td><td>[110, 160, 200, 470, 510]</td><td>360</td></tr><tr><td>31</td><td>[150, 230, 260, 460, 480]</td><td>190</td><td>81</td><td>[120, 150, 240, 280, 380]</td><td>200</td></tr><tr><td>32</td><td>[70, 100, 160, 320, 380]</td><td>75</td><td>82</td><td>[70, 110, 120, 220, 490]</td><td>150</td></tr><tr><td>33</td><td>[80, 250, 350, 400, 500]</td><td>300</td><td>83</td><td>[150, 240, 380, 480, 510]</td><td>170</td></tr><tr><td>34</td><td>[70, 90, 140, 270, 380]</td><td>170</td><td>84</td><td>[140, 210, 240, 450, 460]</td><td>340</td></tr><tr><td>35</td><td>[110, 130, 300, 330, 510]</td><td>250</td><td>85</td><td>[120, 210, 220, 440, 450]</td><td>140</td></tr><tr><td>36</td><td>[260, 340, 350, 440, 460]</td><td>330</td><td>86</td><td>[70, 75, 220, 290, 340]</td><td>100</td></tr><tr><td>37</td><td>[100, 290, 320, 390, 460]</td><td>250</td><td>87</td><td>[75, 150, 160, 240, 300]</td><td>95</td></tr><tr><td>38</td><td>[160, 270, 330, 370, 400]</td><td>350</td><td>88</td><td>[120, 150, 250, 480, 510]</td><td>270</td></tr><tr><td>39</td><td>[90, 190, 230, 420, 490]</td><td>410</td><td>89</td><td>[80, 110, 140, 160, 440]</td><td>85</td></tr><tr><td>40</td><td>[160, 170, 310, 370, 480]</td><td>350</td><td>90</td><td>[85, 200, 340, 490, 510]</td><td>450</td></tr><tr><td>41</td><td>[80, 95, 420, 480, 500]</td><td>260</td><td>91</td><td>[80, 120, 170, 370, 470]</td><td>160</td></tr><tr><td>42</td><td>[120, 190, 230, 490, 500]</td><td>410</td><td>92</td><td>[230, 280, 340, 420, 480]</td><td>300</td></tr><tr><td>43</td><td>[80, 110, 270, 440, 490]</td><td>150</td><td>93</td><td>[100, 220, 250, 390, 450]</td><td>210</td></tr><tr><td>44</td><td>[85, 170, 330, 440, 510]</td><td>430</td><td>94</td><td>[190, 360, 440, 450, 470]</td><td>320</td></tr><tr><td>45</td><td>[85, 140, 240, 350, 390]</td><td>150</td><td>95</td><td>[300, 320, 340, 470, 510]</td><td>390</td></tr><tr><td>46</td><td>[95, 250, 280, 390, 400]</td><td>300</td><td>96</td><td>[70, 75, 220, 250, 450]</td><td>120</td></tr><tr><td>47</td><td>[120, 180, 240, 340, 480]</td><td>170</td><td>97</td><td>[75, 80, 210, 350, 360]</td><td>90</td></tr><tr><td>48</td><td>[140, 200, 320, 440, 510]</td><td>410</td><td>98</td><td>[90, 190, 300, 450, 460]</td><td>370</td></tr><tr><td>49</td><td>[70, 220, 230, 310, 460]</td><td>150</td><td>99</td><td>[75, 85, 95, 350, 500]</td><td>130</td></tr><tr><td>50</td><td>[95, 120, 250, 460, 470]</td><td>390</td><td>100</td><td>[70, 100, 290, 310, 320]</td><td>90</td></tr></tbody></table>