Skip to main content

Apprentissage et calcul scientifique

Emmanuel Franck

Section 17.2 Réduction de dimension pour les problèmes temporelles

Dans cette section on va introduire les méthodes de réduction linéaires pour des EDP temporelles. On va considérer comme modèle l’équation de Burgers visqueuse qui est un bon candidat pour introduire un certains nombre de choses.
Soit \(\rho(t,x)\) solution de
\begin{equation} \left\{\begin{array}{l} \partial \rho +\partial_x f(\rho)=\frac{1}{R_e}\partial_{xx} \rho\\ \rho(0,x)=\rho(x,\alpha)\\ \rho(t,x=0)=r_l, \quad \rho(t,x=L)=r_r \end{array}\right. \tag{17.15} \end{equation}
avec \(\alpha\in [a,b]\) une paramétrisation de la condition initiale et \(R_e\in [c,d]\) le nombre de Reynolds. L’enjeu va être de construire un modèle réduit valide pour toutes valeurs de \(\alpha\) et \(R_e\) dans les intervaux considérés.
Cette équation est parabolique et hyperbolique à la limite \(R_e \rightarrow \infty\text{.}\)
On voit qu’on va essayer de construire un modèle réduit valide sur une gamme de plus en plus large de configurations. On va commencer par introduire plusieurs méthodes de construction des opérateurs de réduction pour les ODE (EDP déjà discrétisées) puis des méthodes de projection pour construire le modèle réduit.

Subsection 17.2.1 Schémas en temps

Plusieurs schémas en temps peuvent être couplés avec les schémas de types volumes finis, éléments finis. On va rapidement rappeler les principales familles qui sont utilisées et qui seront aussi utilisées pour les modèles réduits. On se propose d’écrire les schémas sous une forme particulière qui nous sera utile pour la suite. Ici on résout une ODE du type
\begin{equation} \frac{d \bs{x}}{dt}=\bs{f}(\bs{x}(t))\tag{17.16} \end{equation}
Elle peut être issue de la discrétisation d’une EDP.

Subsubsection 17.2.1.1 Méthodes de Runge-Kutta

Les méthodes de Runge-Kutta sont des méthodes d’ordre élevées qui calculent le le temps courant uniquement à partir du temps précédent. Cependant elle calculent les quantités et les flux à des temps intermédiaires.
Définition 17.8. Méthodes de Runge-Kutta.
La méthode de Runge Kutta s’écrit en deux étapes: une première ou on calcul un résidu solution de:
\begin{equation} \bs{r}_{i}^{n}\left(\bs{w}_{1}^{n}, \ldots, \bs{w}_{s}^{n}\right)=0, \quad i \in \mathbb{N}(s)\tag{17.17} \end{equation}
avec
\begin{equation} \bs{r}_{i}^{n}\left(\bs{w}_{1}, \ldots, \bs{w}_{s}\right) :=\bs{w}_{i}-\bs{f}\left(\bs{x}^{n-1} +\Delta t \sum_{j=1}^{s} a_{i j} \bs{w}_{j}, t^{n-1}+c_{i} \Delta t\right), \quad i \in \mathbb{N}(s)\tag{17.18} \end{equation}
et ensuite on en déduit la solution
\begin{equation} \bs{x}^n=\bs{x}^{n-1}+\Delta t \sum_{i=1}^s b_i \bs{w}_i^n \tag{17.19} \end{equation}

Subsubsection 17.2.1.2 Méthodes temporelles multi-pas

Il s’agit de méthodes qui utilisent l’approximation à plusieurs pas de temps. Ici on va considérer les méthodes explicites.
Définition 17.9. Méthodes linéaires multi-pas.
Une méthode linéaire a \(k\) pas s’écrit en deux étapes: une première ou on calcul un résidu solution de:
\begin{equation} \bs{r}\left(\bs{w}\right)=0\tag{17.20} \end{equation}
avec
\begin{equation} \bs{r}(\bs{w})=\alpha_0 \bs{w} - \Delta t \beta_0 \bs{f}(\bs{w},t) + \sum_{j=1}^{K} \alpha_{j} \bs{x}^{n-j}- \Delta t \sum_{j=1}^{K} \beta_{j} \bs{f}\left(\bs{x}^{n-j}, t^{n-j}\right),\tag{17.21} \end{equation}
avec comme condition \(\sum_{j=0}^K\alpha_j=0\) et la solution est donnée par
\begin{equation} \bs{x}_n=\bs{w} \tag{17.22} \end{equation}
La plus connue est la méthode d’Euler explicite ou \(K=1\) et \(\alpha_0/\beta_0=(1,0)\text{,}\) \(\alpha_1/\beta_1=(1,1)\text{.}\) La méthode d’Euler implicite correspond-elle à \(K=1\) et \(\alpha_0/\beta_0=(1,1)\text{,}\) \(\alpha_1/\beta_1=(1,0)\text{.}\) Pour finir la méthode classe de Crank-Nicolson correspond à \(K=1\) et \(\alpha_0/\beta_0=(1,\frac12)\text{,}\) \(\alpha_1/\beta_1=(1,\frac12)\text{.}\) Les méthodes dites de Adams–Bashforth, Adam-Moulton (explicites) ou BDF (implicite) rentrent dans ce cadre.

Subsection 17.2.2 Méthodes de réduction pour les ODE

On va commencer ici par introduire plusieurs méthodes pour construire l’opérateur. On commence par introduire une ODE général qui sera souvent issue de la discrétisation d’une EDP: On considère le modèle suivant
\begin{equation} \frac{d \bs{x}(t,\bs{\mu})}{dt}= \bs{F}(\bs{x}(t,\bs{\mu}),\bs{\mu})\tag{17.23} \end{equation}
avec \(\bs{x}(t)\in \mathbb{R}^n\text{.}\) Comme dans le cas elliptique, pour construire un modèle réduit, on fait l’hypothèse fondamentale suivante:

Hypothèse fondamentale.

Il existe un sous-espace vectoriel \(\mathcal{S}:=\operatorname{span}\left\{\phi_{i}\right\}_{i=1}^{m} \subseteq \mathbb{R}^{n}\) tel que \(m \ll n\) et tel que
\begin{equation} \bs{x}(t,\bs{\mu}) \approx \bs{x}_{ref}(\bs{\mu})+\Phi \hat{\bs{x}}(t,\bs{\mu})\tag{17.24} \end{equation}
Cela revient à dire que les solutions de notre EDO en grande dimension peuvent être representée par des données vivant dans espace de dimension petit (ici donc un hyperplan de dimension \(m\)). L’enjeu maintenant est de construire l’opérateur linéaire \(\Phi\) qui relie les données en grande dimension et en petite dimension.

Subsubsection 17.2.2.1 La méthode POD

La méthode POD introduite à la section précédente peut être utilisée de la même façon. Seule la matrice des schnapshots va changer. On considère ici la matrice:
\begin{equation} X=[\bs{x}(t_0,\bs{\mu}_1)-\bs{x}_{ref}(\bs{\mu}) ,...,\bs{x}(t_{N_t},\bs{\mu}_1)-\bs{x}_{ref}(\bs{\mu}) ,...,\bs{x}(t_0,\bs{\mu}_{N_{\mu}})-\bs{x}_{ref}(\bs{\mu}) ,...,\bs{x}(t_{N_t},\bs{\mu}_{N_{\mu}})-\bs{x}_{ref}(\bs{\mu})]\tag{17.25} \end{equation}
avec \(X \in \mathcal{M}_{n\times d}(\mathbb{R})\) et \(d=(N_t+1) N_{\mu}\text{.}\)

Subsubsection 17.2.2.2 La méthode DMD

La méthode de la décomposition orthogonale en modes propres n’utilise pas le fait que le matrice de snapshots qu’on cherche à approcher est issue d’une série temporelle. En effet si on permutte les snapshots et casse la structure temporelle le résultat reste inchanger par conséquent cela veut dire que la réduction ne tient pas compte de la structure temporelle du problème. Pour remédier à cela une variante a été introduite appelée: décomposition en modes dynamiques ou DMD. Cette approche est a priori faite pour compresser une trajectoire temporelle (une valeur des paramètres). On va dans l’introduire pour ce type de problème avant de montrer comment traiter le cas plus général. On se donne les matrices de snapshots suivantes:
\begin{equation} X=[\bs{x}(t_0),...,\bs{x}(t_{N_t}) ]\tag{17.26} \end{equation}
Puisque ces snapshots sont issus de la résolution d’une EDO on peut supposer que
\begin{equation*} \bs{x}_{i+1}=A\bs{x}_{i} \end{equation*}
avec \(A\) un opérateur serait en pratique notre schéma de discrétisation si on traite une EDP linéaire. Si l’EDP est nonlinéaire ou si on traite des données réels on ne connais pas a priori \(A \text{.}\) On définit les matrices
\begin{equation*} X_1^{N_t-1}=[\bs{x}(t_0),...,\bs{x}(t_{N_t-1}) ] \end{equation*}
et
\begin{equation*} X_2^{N_t}=[\bs{x}(t_1),...,\bs{x}(t_{N_t}) ] \end{equation*}
Dans ce cas la, on va donc, en premier lieu, construire l’opérateur \(A\) en résolvant un problème de moindre carrés de la forme:
\begin{equation} \underset{Rg(A)=K}{\operatorname{min}}\parallel X_2^{N_t} - A X_1^{N_t-1}\parallel_F^2\tag{17.27} \end{equation}
Dans ce cas notre relation se réécrit
\begin{equation} X_2^{N_t}=AX_1^{N_t-1} +re_{N_t-1}\tag{17.28} \end{equation}
avec \(e_{N_t-1}\) un vecteur nul avec un \(1\) à la dernière composante et \(r\) un résidu qui modélise les comportements qu’on ne peut pas complètement décrire avec \(A\text{.}\) L’idée de la méthode DMD est de chercher une approximation de rang faible de l’opérateur \(A\text{.}\) Il existe plusieurs approches on va introduire l’approche utilisant la SVD.
Définition 17.10.
On se fixe une dimension de compression \(k\text{.}\) On se donne \(X_1^{N_t-1}=U \Sigma V^t\text{.}\) La solution de la méthode DMD \(\Phi_K\in \mathcal{M}_{n \times k}(\mathbb{R})\) est donnée par
\begin{equation*} \Phi_K= U y_K \end{equation*}
avec \(y_K\) les \(K\) premiers vecteurs propres de \(U^tX_2^{N_t}V \Sigma^{-1}\) et \(X_1^{N_t-1}=U \Sigma V^t\text{.}\)
On va pas donner de preuve de ce résultat dans le sens d’un résultat d’optimalité par rapport à la minimisation d’un coup, mais plutôt les principes de construction. On part de la formule
\begin{equation*} X_2^{N_t}=AX_1^{N_t-1} +re_{N_t-1} \end{equation*}
On applique la SVD sur \(X_1^{N_t-1}\) pour obtenir
\begin{equation*} X_2^{N_t}=A U \Sigma V^t +re_{N_t-1} \end{equation*}
et on multiplie par \(U^t\) pour obtenir
\begin{equation*} U^t X_2^{N_t}=U^t A U \Sigma V^t + U^t re_{N_t-1} =U^t A U \Sigma V^t \end{equation*}
car on suppose que \(U^t r=0\text{.}\) En appliquant deux autres multiplication (et en n’oubliant pas que \(V\) est orthogonale) on obtient
\begin{equation*} U^t X_2^{N_t}V \Sigma^{-1} =U^t A U \end{equation*}
La matrice de gauche est similaire a celle de droite puisqu’elle peut être obtenue par changement de base orthogonale à partir de A. On peut en déduire qu’elle partage les mêmes valeurs propres et que \(U\) fois les vecteurs propres de la matrice de gauche donnent les vecteurs propres de \(A\text{.}\) Pour approcher \(A\) par une matrice de rang \(K\) le choix naturel est de prendre les \(K\) premiers vecteurs propres de \(A\) et on voit que donc que cela revient à prendre \(U y_K\) avec \(y_K\) les \(K\) premiers vecteurs de \(U^t X_2^{N_t}V \Sigma^{-1}\text{.}\) Cette méthode nous permet d’avoir un opérateur de réduction pour une trajectoire temporelle, mais maintenant on voudrait l’utiliser dans le cas ou on souhaite compresser notre ODE pour plusieurs valeurs de paramètres. On nomme
\begin{equation*} X_{\mu_i}=[\bs{x}^{\mu_i}(t_0),...,\bs{x}^{\mu_i}(t_{N_t-1}) ] \in \mathcal{M}_{n\times N_t}(\mathbb{R}) \end{equation*}
notre matrice de snapshots pour le paramètre \(\mu_i\) et
\begin{equation*} X=[X_{\mu_1},...,X_{\mu_{N_{\mu}}}]\in \mathcal{M}_{n\times (N_t N_{\mu}}((\mathbb{R})) \end{equation*}
la matrice complète. A partir de là, on applique la POD. La méthode nous donne un opérateur \(\Phi = U_{k_t}\in M_{n,K}(\mathbb{R})\) avec \(U\) les vecteurs propres de \(XX^t\text{.}\) On choisit \(k_t\) de façon à représenter de façon optimale possible le matrice d’origine \(X\text{.}\) A partir de la on définit
\begin{equation*} \tilde{\bs{x}}^{\mu_i}(t_j)= U_n^t \bs{x}^{\mu_i}(t_j) \end{equation*}
Cela revient à appliquer la compression pour réduite la dimension de chaque snapshots temporel. A partir de la on définit
\begin{equation*} \tilde{X}_{\mu_i}= U_n^t X_{\mu_i} \end{equation*}
et
\begin{equation*} X_2=[\tilde{X}_{\mu_1} , ...., \tilde{X}_{\mu_{N_{mu}}} ]^t \end{equation*}
À ce moment-là on a gardé la structure temporelle. En effet dans chaque ligne il y a l’ensemble des snapshots temporels mais pour les variables réduites (par POD) et pour une valeur de paramètre donnée. On peut voir cela comme les snapshots temporelle d’une EDO en petite dimension (puisqu’il y a eu réduction) paramétrique ( dans le sens ou elle est dupliquée pour \(N_{\mu}\) valeurs des paramètres. On a donc au final les snapshots d’une ODE de taille \(k_t N_{\mu}\) ce qui peut être encore assez grand. A partir de la on peut appliquer une méthode DMD qui compresser cette grande ODE en tenant compte la structure temporelle. Pour résumer, on pouvait voir notre matrice de snapshots comme des temps d’une ODE de très grande taille (\(n N_{\mu}\text{,}\) la taille \(n\) pouvant venir d’une discrétisation spatiale). On compresse donc d’abord "spatialement" donc de \(n\) a \(k_t\) puis on compresse par rapport aux paramètres en tenant compte du temps grâce à la méthode DMD (ref "A dynamic mode decomposition extension for the forecasting of parametric dynamical system").
Pour améliorer certaines performances selon le système qu’on regarde, on peut utiliser une variante de la méthode DMD appelée la DMD informé physiquement. L’idée de l’approche est d’utiliser de qu’on connait de la physique fin de construire la matrice \(A\) dans la représentation (17.28). Pour cela, au lieu de résoudre (17.27) on va résoudre
\begin{equation} \underset{A\in \mathcal{M}}{\operatorname{min}}\parallel X_2^{N_t} - A X_1^{N_t-1}\parallel_F^2\tag{17.29} \end{equation}
avec \(\mathcal{M}\) un sous espace ou une sous-variété de matrices qui va dépendre de la physique. La question est maintenant de savoir comment caractériser cette variété \(\mathcal{M}\text{.}\) Prenons l’exemple d’une équation de transport du type:
\begin{equation*} \partial_t u + a\partial_x u =0 \end{equation*}
Si on discretise cette équation on obtiendra un schéma du type
\begin{equation*} U^{n+1}= A U^n \end{equation*}
puisqu’il s’agit d’ un problème de transport la solution va juste être translatée en temps. Si on applique un opérateur de transport \(S_{\phi}\) on a \(A S_{\phi}=S_{\phi} A \text{.}\) En effet qu’on applique une translation avant ou après l’opérateur temporelle c’est équivalent puisque l’opérateur temporelle est lui même un opérateur de translation. On voit donc que pour ce type de problème on peut caractériser \(\mathcal{M}\) de la façon suivante:
\begin{equation*} \mathcal{M}= \left\{ A, \mbox{such that} A S_{\phi}=S_{\phi} A, \forall \phi \in \mathbb{R}\right\} \end{equation*}
Cela est équivalent au matrice circulente. Sur la figure Figure 17.11 on peut voir différent type de sous espace de matrices.
Figure 17.11. Exemple de matrices particulières pouvant être utilisée par la méthode PiDMD (Physic informed DMD). Image issue de [1.24]
Ce type de matrice peut être relié a des équations ou système de type EDP. Des exemples sont données sur la figure Figure 17.12
Figure 17.12. Exemple de matrice particulières associées a des équations classiques. Image issue de [1.24].
Les détails pour résoudre le problème (17.29) sont introduits dans [1.24].

Subsubsection 17.2.2.3 Méthode Gloutonne

La méthode Gloutonne introduite dans le cas elliptique peut aussi s’utiliser dans le cas temporelle. On commence par modifier la matrice de snapshots qui devient:
\begin{equation} X=[\bs{x}(t_0,\bs{\mu}_1),...,\bs{x}(t_{N_t}, \bs{\mu}_1),...,\bs{x}(t_0,\bs{\mu}_{N_{\mu}}),..., \bs{x}(t_{N_t},\bs{\mu}_{N_{\mu}})]\tag{17.30} \end{equation}
Bien que l’algorithme reste sensiblement le même dans me cas temporel, on se propose de le réécrire.
Pour finir la méthode il suffit de savoir calculer la projection. Pour cela on utilise une projection de Galerkin. On écrit
\begin{equation*} \Pi_{\mathcal{B}}(\bs{f})=\sum_{l=1}^K\alpha_f^l \phi_l \end{equation*}
et on cherche à déterminer \((\alpha_f^l,...\alpha_f^K)\) tel que
\begin{equation*} \sum_{l=1}^K\alpha_f^l \phi_l= \bs{f} \end{equation*}
Ensuite on va multiplier cette relation par \((\phi_1,... \phi_l)\text{.}\) La base étant orthonormale on obtient naturellement que
\begin{equation*} \alpha_f^l=(\bs{f},\phi_l) \end{equation*}
et donc
La projection sur la base gloutonne est donnée par
\begin{equation} \Pi_{\mathcal{B}}(\bs{f})=\sum_{l=1}^K (\bs{f},\phi_l) \phi_l\tag{17.31} \end{equation}
On voit donc la phase d’orthonormalisation est très importante pour simplifier les calculs.

Subsubsection 17.2.2.4 Méthode POD-gloutonne

La dernière méthode combine les approches POD et gloutonne. Elle a été introduite dans [1.2] L’idée centrale est d’utiliser une approche gloutonne pour traiter la dépendance aux paramètres et une approche de type POD pour le temps. Voici les étapes de l’algorithme.

Subsection 17.2.3 Construction du modèle réduit

Maintenant que nous avons déterminé notre décodeur \(\Phi_K\) et (l’encodeur qui va avec) ce qui implicitement caractérise le sous-espace de basse dimension sur lequel "vivent" les données, il faut construire le modèle sur ces données réduites. Là encore on va introduire deux approches classiques. Nos méthodes de réduction précédente reviennent à construire \(\tilde{\bs{x}}\in \mathbb{R}^n\) approximation de \(\bs{x}\in \mathbb{R}^n\) satisfaisons l’hypothèse (17.24). Dans la suite on fera l’hypothèse qu’on connait \(\bs{x}_{ref}\) et \(\Phi\text{.}\) Il y a deux approches: les méthodes de projection et les méthodes d’inférence d’opérateur. Ici on va considérer que les méthodes de projection. On parlera des méthodes d’inférences plus tard.

Subsubsection 17.2.3.1 Méthode de Projection: Least-Square Galerkin

La première approche consiste à appliquer au modèle une projection dite de "Galerkin" ou "L2". Pour cela on réécrit le problème global sous la forme:
\begin{equation} \bs{r}(\frac{d\bs{x}}{dt}, \bs{x}, t ; \bs{\mu}):=\frac{d\bs{x}}{dt}-\bs{f}(\bs{x}, t; \bs{\mu})=0, \quad \bs{x}(0 ; \bs{\mu})=\bs{x}_{0}(\bs{\mu}) \tag{17.32} \end{equation}
avec \(\bs{r}: \mathbb{R}^{n} \times \mathbb{R}^{n} \times \mathbb{R}_{+} \times \mathcal{D} \rightarrow \mathbb{R}^{n}\) avec \((\dot{\bs{w}}, \bs{w}, \tau ; \bs{\nu}) \mapsto \bs{r}(\dot{\bs{w}}, \bs{w}, \tau ; \bs{\nu})\) correspond au résidu de notre problème continu en temps. Maintenant l’idée va tout simplement être de remplacer \(\bs{x}\) par son approximation affine (17.24). On pose \(\tilde{\bs{x}}\) les variables reconstruites:
\begin{equation*} \tilde{\bs{x}}(t,\bs{\mu}) = \bs{x}_{ref}(\bs{\mu})+\Phi \hat{\bs{x}}(t,\bs{\mu}) \end{equation*}
On définit aussi un nouveau résidu qu’on va appeler le résidu réduit:
\begin{equation} \hat{\mathbf{r}}(\frac{d\hat{\bs{x}}}{dt}, \hat{\bs{x}}, t ; \bs{\mu}):= \bs{r}\left(\Phi_K \frac{d\hat{\bs{x}}}{dt}, \bs{x}_{ref}+\Phi_K \hat{\bs{x}}, t ; \bs{\mu}\right)\tag{17.33} \end{equation}
Pour obtenir ce résultat on utilise que \(\frac{d \tilde{\bs{x}}}{dt}=\Phi_K\frac{d\hat{\bs{x}}}{dt}\text{,}\) car \(\bs{x}_{ref}(\bs{\mu})\) ne dépend pas du temps.
Définition 17.15. Projection de Galerkin (LSG).
La projection de Galerkin est donnée par la solution du problème
\begin{equation} \dot{\tilde{\bs{x}}}= \underset{\bs{v} \in \operatorname{Vect}(\Phi_k)}{\operatorname{argmin}}\|\mathbf{r}(\bs{v}, \tilde{\bs{x}}, t ; \bs{\mu})\|_{2}^{2}\tag{17.34} \end{equation}
Ce problème est équivalent à
\begin{equation} \dot{\hat{\bs{x}}}= \underset{\hat{\bs{v}} \in \mathbb{R}^m}{\operatorname{argmin}}\|\hat{\mathbf{r}}(\hat{\bs{v}}, \tilde{\bs{x}}, t ; \bs{\mu})\|_{2}^{2}\tag{17.35} \end{equation}
avec \(\hat{\bs{x}}(0 ; \bs{\mu})= \hat{\bs{x}}_{0}(\bs{\mu})=\Phi_K^{T}\left(\bs{x}_{0}(\bs{\mu})-\bs{x}_{r e f}\right)\text{.}\)
On fait deux choix ici: représenter la variable réduite sur l’espace \(\operatorname{vect}(\Phi_k)+\bs{x}_{ref}\) (espace de fonction "trial") et projeter la dérivée en temps sur l’espace \(\operatorname{vect}(\Phi_k)\) (espace de fonction test).
La méthode de Galerkin revient à trouver l’expression de notre dérivée sur variables réduites qui va minimiser le résidu réduit. Si on avait notre résidu qui était nul ce serait équivalent a dire que notre approximation (17.24) ne ferait aucune erreur une fois incorporée dans l’EDO. Ce n’est en pratique pas possible, mais on voit qu’il naturelle d’essayer de minimiser le résidu issu de notre approximation introduite dans l’EDP.
Preuve.
On chercher a minimiser le problème:
\begin{equation*} h(\hat{\bs{v}})=\parallel \Phi_K \hat{\bs{v}}- \bs{f} \left(\bs{x}_{ref}+\Phi_K \hat{\bs{x}}, t ; \bs{\mu}\right)\parallel_2^2 \end{equation*}
On développe:
\begin{equation*} h(\hat{\bs{v}})= \hat{\bs{v}}^{T} \Phi_K^{T} \Phi_K \hat{\bs{v}}-2 \hat{\bs{v}}^{T}\Phi_K^{T} \bs{f}\left(\bs{x}_{ref}+\Phi_K \hat{\bs{x}}, t;\bs{\mu}\right)+\bs{f}\left(\bs{x}_{ref}+\Phi_K \hat{\bs{x}}, t;\bs{\mu}\right)^{T} \bs{f}\left(\bs{x}_{red}+\Phi_K \hat{\bs{x}}, t; \bs{\mu}\right) \end{equation*}
On calcul le gradient
\begin{equation*} \frac{\partial h(\hat{\bs{v}})}{\partial \hat{\bs{v}}}= 2 \Phi_K^{t} \Phi_K \hat{\bs{v}}^{\star}-2 \Phi_K^{t} \bs{f} \left(\bs{x}_{ref}+\Phi_K \hat{\bs{x}}, t;\bs{\mu}\right) \end{equation*}
On résout maintenant \(\frac{\partial h(\hat{\bs{v}})}{\partial \hat{\bs{v}}}=0\) Pour cela on utilise le fait que \(\Phi\) est orthogonale, ce qui donne comme solution \(0= \hat{\bs{v}}^{\star}- \Phi_K^t\bs{f}\left(\bs{x}_{ref}+\Phi_K \hat{\bs{x}}, t;\bs{\mu}\right)\)
Si on considère une ODE linéaire (\(\bs{f}(\hat{\bs{x}}, t ; \bs{\mu})=A(\bs{\mu})\bs{x}\)) on voit qu’il faut calculer les deux termes:
\begin{equation*} \Phi_K^{T} A(\bs{\mu}) \bs{x}_{ref}, \quad \Phi_K^{T} A(\bs{\mu})\Phi_K \end{equation*}
Ces deux termes peuvent être pré-calculé et stocké afin de pouvoir effectuer les calculs sur le modèle réduit. Dans le cas nonlinéaire, la situation est différente et l’abordera par la suite. Il existe une variante à cette projection ou on ne résout pas le problème de minimisation en norme \(\| . \|_2\) mais une norme pondérée. On parle de projection de Galerkin a Poids.
Définition 17.17. Projection de Galerkin à poids (WLSG).
La projection de Galerkin à poids est donnée par la solution du problème
\begin{equation} \dot{\hat{\bs{x}}}= \underset{\hat{\bs{v}} \in \mathbb{R}^{n}_{s}} {\operatorname{argmin}}\|\tilde{\mathbf{r}}(\hat{\bs{v}}, \hat{\bs{x}}, t ; \bs{\mu})\|_{\Theta}^{2}\tag{17.37} \end{equation}
avec \(\hat{\bs{x}}(0 ; \bs{\mu})=\hat{\bs{x}}_{0}(\bs{\mu}) =\Phi_K^{T}\left(\bs{x}_{0}(\bs{\mu})-\bs{x}_{r e f}\right)\) et \(\|x\|_{\Theta} \equiv(x, x)_{\Theta}\) avec \((x, y)_{\Theta} \equiv \sqrt{x^{T} \Theta y}\) et \(\Theta \equiv \Theta^{T / 2} \Theta^{1 / 2}\text{.}\)
Preuve.
on fait les mêmes calculs que pour la méthode LSG.
Maintenant que le modèle est établi on peut le discrétiser avec un schéma temporel de notre choix. On va introduire les schémas en temps de type multi-pas et de Runge-Kutta pur le modèle réduit comme on l’avait fait pour le modèle complet dans la sous-section Subsection 17.2.3.
Définition 17.19. Méthodes de Runge-Kutta pour la méthode LSG.
La méthode de Runge Kutta s’écrit en deux étapes: une première ou on calcul un résidu solution de:
\begin{equation} \bs{r}_{i}^{n}\left(\hat{\bs{w}}_{1}^{n}, \ldots, \hat{\bs{w}}_{s}^{n}\right)=0, \quad i \in \mathbb{N}(s)\tag{17.38} \end{equation}
avec
\begin{equation} \bs{r}_{i}^{n}\left(\hat{\bs{w}}_{1}, \ldots, \hat{\bs{w}}_{s}\right) :=\hat{\bs{w}}_{i}-\Phi_K^{T}\bs{f}\left(\bs{x}_{ref}+\Phi_K\hat{\bs{x}}^{n-1}+ \Delta t \sum_{j=1}^{s} a_{i j}\Phi_K\hat{\bs{w}}_{j}, t^{n-1}+c_{i} \Delta t\right), \quad i \in \mathbb{N}(s)\tag{17.39} \end{equation}
et ensuite on en déduit la solution
\begin{equation} \hat{\bs{x}}^n=\hat{\bs{x}}^{n-1}+\Delta t \sum_{i=1}^s b_i \hat{\bs{w}}_i^n\tag{17.40} \end{equation}
Définition 17.20. Méthodes multi-pas pour la méthode LSG.
Une méthode linéaire a \(s\) pas s’écrit en deux étapes: une première ou on calcul un résidu solution de:
\begin{equation} \bs{r}(\hat{\bs{w}})=0\tag{17.41} \end{equation}
avec
\begin{equation} \bs{r}^n(\hat{\bs{w}})=\alpha_0 \hat{\bs{w}} - \Delta t \beta_0 \Phi_K^{T} \bs{f}\left(\bs{x}_{ref}+\Phi_K \hat{\bs{w}}, t^n ; \bs{\mu}\right) + \sum_{j=1}^{s} \alpha_{j} \hat{\bs{x}}^{n-j}-\Delta t \sum_{j=1}^{s} \beta_{j} \Phi_K^{T} \bs{f} \left(\bs{x}_{ref}+\Phi_K \hat{\bs{x}}^{n-j}, t^{n-j} ; \bs{\mu}\right),\tag{17.42} \end{equation}
avec comme condition \(\sum_{j=0}^s\alpha_j=0\) et la solution est donnée par
\begin{equation} \hat{\bs{x}}_n=\hat{\bs{w}}.\tag{17.43} \end{equation}
Ces formulations reviennent juste à dire qu’on applique nos schémas classiques à l’EDO réduite issue de la projection de Galerkin.

Subsubsection 17.2.3.2 Méthode de projection: Least-Square Petrov-Galerkin (LSPG)

La méthode précédente utilise la minimisation du résidu de l’équation continue pour construire le modèle réduit. Maintenant nous allons utiliser une alternative qui cherche à minimiser le résidu du problème discret en temps (on discrétise avant de construire le modèle réduit). On rappelle fait ici l’hypothèse donc que notre solution approchée satisfait l’hypothèse (17.24).
On se donne \(\bs{r}^n()\) le résidu d’un schéma en temps. Par exemple le résidu des schémas linéaires multi-pas (17.20)-(17.21)-(17.22). L’objectif de la projection de Petrov Galerkin est de minimiser le résidu discret à poids suivant:
\begin{equation} \tilde{\bs{w}}^n=\operatorname{min}_{\bs{w}\in \bs{x}_{ref}+Vect(\phi_1,...,\phi_K)}\parallel A(\bs{w})\bs{r}^n(\bs{w}) \parallel_2^2\tag{17.44} \end{equation}
ce qui est équivalent à
\begin{equation} \hat{\bs{x}}^n=\operatorname{min}_{\hat{\bs{w}}\in \mathbb{R}^k} J(\hat{\bs{w}})= \operatorname{min}_{\hat{\bs{w}}\in \mathbb{R}^k}\parallel A(\bs{x}_{ref}+\Phi_K\hat{\bs{w}})r^n(\bs{x}_{ref}+ \Phi_K\hat{\bs{w}}) \parallel_2^2\tag{17.45} \end{equation}
avec \(A\in \mathcal{M}_{l,n}(\mathbb{R})\) une matrice de poids avec \(l\leq n\text{.}\) On peut évidemment prendre pour \(A\) l’identité ou une restriction au \(l\) premières lignes de l’identité, mais d’autres choix sont aussi possibles.
Preuve.
On commence par calculer le gradient de \(J(\hat{\bs{w}})\) par rapport à \(\hat{\bs{w}}\text{.}\) On pose \(\tilde{\bs{w}}=\bs{x}_{ref}+\Phi_K\hat{\bs{w}}\) On calcul
\begin{equation*} J(\tilde{\bs{w}} +h )= \langle Ar^n(\tilde{\bs{w}} +\bs{h}) ,Ar^n(\tilde{\bs{w}}+\bs{h}) \rangle_{\mathbb{R}^l} \end{equation*}
\begin{equation*} J(\tilde{\bs{w}} +h )= \langle A(r^n(\tilde{\bs{w}}) +\nabla_{\bs{w}}\bs{r}^n \bs{h} ) , A(r^n(\tilde{\bs{w}}) +\nabla_{\bs{w}}\bs{r}^n \bs{h} ) \rangle_{\mathbb{R}^l} \end{equation*}
On développe et ce la donne
\begin{equation*} J(\tilde{\bs{w}} +h )= J(\tilde{\bs{w}})+ \langle A \bs{r}^n , A (\nabla_{\bs{w}}\bs{r}^n) \bs{h} \rangle_{\mathbb{R}^l} \end{equation*}
ce qui se réécrit
\begin{equation*} J(\tilde{\bs{w}} +h )= J(\tilde{\bs{w}})+ (\nabla_{\bs{w}}\bs{r}^n)^t \langle A^t A \bs{r}^n , \bs{h} \rangle_{\mathbb{R}^n} + o(h^2) \end{equation*}
On obtient donc par identification
\begin{equation*} \nabla_{\tilde{\bs{w}}}J(\tilde{\bs{w}}) \bs{h} = \langle (\nabla_{\bs{w}}r^n)^t A^t A \bs{r}^n , \bs{h} \rangle_{\mathbb{R}^n} \end{equation*}
et donc
\begin{equation*} \nabla_{\tilde{\bs{w}}}J(\tilde{\bs{w}}) = (\nabla_{\bs{w}}\bs{r}^n)^t A^t A \bs{r}^n \end{equation*}
Ensuite on utilise que
\begin{equation*} \nabla_{\hat{\bs{w}}}J(\hat{\bs{w}}) =\Phi_k^t \nabla_{\tilde{\bs{w}}}J(\tilde{\bs{w}}) \end{equation*}
ce qui nous donne
\begin{equation*} \nabla_{\hat{\bs{w}}}J(\hat{\bs{w}}) =\Phi_k^t (\nabla_{\bs{w}}r^n)^t A^t A \bs{r}^n =\Psi_K(\hat{\bs{w}}^n) \bs{r}^n \end{equation*}
Ensuite on obtient que \(\nabla_{\bs{w}} J(\hat{\bs{w}})=0\) est équivalent au résultat.
Ce résultat n’est pas immédiatement interprétable. En général, afin de comprendre la différence avec la méthode de Galerkin on propose regarder l’équation continue avec lequel la méthode de moindre carré Petrov-Galerkin pour les schémas multi-pas est consistante.
Preuve.
On applique la projection de Galerkin à l’ODE (17.23). On choisit donc \(\operatorname{vect}(\Phi_k)+\bs{x}_{ref}\) comme espace "trial" et \(\operatorname{vect}(\Psi_k)\) comme espace de fonction test. ON minimise donc
\begin{equation} \dot{\tilde{\bs{x}}}= \underset{\bs{v} \in \operatorname{Vect}(\Psi_k)}{\operatorname{argmin}}\|\mathbf{r}(\bs{v}, \hat{\bs{x}}, t ; \bs{\mu})\|_{2}^{2}\tag{17.48} \end{equation}
Ce qui est équivalent à
\begin{equation} \Phi_k\dot{\hat{\bs{x}}}= \underset{\hat{\bs{v}} \in \mathbb{R}^p}{\operatorname{argmin}}\|\mathbf{r}(\Psi_k \hat{\bs{v}}, \hat{\bs{x}}, t ; \bs{\mu})\|_{2}^{2}\tag{17.49} \end{equation}
Comme dans le cas de LSG on calcule le gradient, on calcul le minimum et on obtient donc
\begin{equation*} \Psi_k (\hat{\bs{x}})^t\Phi_k \frac{d \hat{\bs{x}}}{dt} =\Psi_k (\hat{\bs{x}})^t\bs{f} \left(\bs{x}_{ref}+\Phi_K \hat{\bs{x}}, t ; \bs{\mu}\right) \end{equation*}
donc
\begin{equation} \frac{d \hat{\bs{x}}}{dt} =(\Psi_k (\hat{\bs{x}})^t\Phi_k)^{-1} \Psi_k (\hat{\bs{x}})^t\bs{f} \left(\bs{x}_{ref}+\Phi_K \hat{\bs{x}}, t ; \bs{\mu}\right)\tag{17.50} \end{equation}
Maintenant on va montrer que dans les 3 cas proposés il y a équivalence entre la discrétisation de (17.50) et (17.46) avec (17.47).
  1. Cas 1:
    On considère ici les cas ou les coefficients \(\beta_j\) pour \(j\gt 0\) sont nuls. Lorsqu’on écrit le schéma pour (17.50), le schéma est solution de l’équation
    \begin{equation*} \bs{r}^n(\bs{w})=\alpha_0 \bs{w}-\Delta t\beta_0(\Psi_k (\hat{\bs{w}})^t\Phi_k)^{-1} \Psi_k (\hat{\bs{w}})^t\bs{f} \left(\bs{x}_{ref}+\Phi_K \hat{\bs{w}}, t ; \bs{\mu}\right)+ \sum_{j=1}^{K} \alpha_{j} \bs{x}^{n-j}=0 \end{equation*}
    On multiple par \((\Psi_k (\hat{\bs{w}})^t\Phi_k)\) ce qui donne
    \begin{equation*} \bs{r}^n(\bs{w})=\alpha_0(\Psi_k (\hat{\bs{w}})^t\Phi_k) \bs{w}-\Delta t\beta_0 \Psi_k (\hat{\bs{w}})^t\bs{f} \left(\bs{x}_{ref}+\Phi_K \hat{\bs{w}}, t ; \bs{\mu}\right)+ \sum_{j=1}^{K} \alpha_{j} (\Psi_k (\hat{\bs{w}})^t\Phi_k)\bs{x}^{n-j}=0 \end{equation*}
    On factorise ce qui donne
    \begin{equation*} \bs{r}^n(\bs{w})= (\Psi_k (\hat{\bs{w}})^t\Phi_k)\left[\alpha_0 \bs{w}-\Delta t\beta_0 \Phi_K^t (\hat{\bs{w}})^t\bs{f} \left(\bs{x}_{ref}+\Phi_K \hat{\bs{w}}, t ; \bs{\mu}\right)+ \sum_{j=1}^{K} \alpha_{j} \bs{x}^{n-j}\right]=0 \end{equation*}
    On retrouve donc bien (17.46) avec le résidu donné par (17.42). Cela donne
    \begin{equation*} \Psi_K(\hat{\bs{w}})^t\bs{r}^n(\bs{x}_{ref}+\Phi_K\hat{\bs{w}})=0 \end{equation*}
    On identifie la dérivée du résidu par rapport à \(\bs{w}\)ce qui nous donne le bon \(\Psi_K\text{.}\) Pour finir on utilise que dans le schéma multi-pas \(\hat{bs{x}}^n=\hat{\bs{w}}\) et on obtient
    \begin{equation*} \Psi_K(\hat{\bs{x}}^n)^t\bs{r}^n(\bs{x}_{ref}+\Phi_K\hat{\bs{x}}^n)=0 \end{equation*}
  2. Cas 2 et 3:
    Puisqu’on on considère des équations autonomes le cas \(\bs{f}\) linéaire et le cas \(\beta_0=0\) sont équivalent on a
    \begin{equation*} \Psi_K=\alpha_{0}A^t A \Phi_k \end{equation*}
    A partir de la on applique un schéma linéaire multi pas à l’équation (17.50). Lorsqu’on écrit le schéma pour (17.50), le schéma est solution de l’équation
    \begin{equation*} \bs{r}^n(\bs{w})=\alpha_0 \bs{w}-\Delta t\beta_0(\Psi_k^t\Phi_k)^{-1} \Psi_k (\hat{\bs{w}})^t\bs{f} \left(\bs{x}_{ref}+\Phi_K \hat{\bs{w}}, t ; \bs{\mu}\right)+ \sum_{j=1}^{K} \alpha_{j} \bs{x}^{n-j} - \Delta t \sum_{j=1}^{K} \beta_{j} (\Psi_k^t\Phi_k)^{-1} \Psi_k \bs{f}\left(\bs{x}^{n-j}, t^{n-j}\right) =0 \end{equation*}
    On multiple par \(\Psi_K^t\Phi_k\) on factorise par \(\Psi_K^t\Phi_k\) et on retrouve le produit entre \(\Psi_k\) et le résidu. La fin de la preuve est comme dans le cas 1.
On peut aussi écrire la méthode LSPG pour les schémas de type Runge-Kutta et faire le même lien que précédemment entre version continue et discrète. Les détails ainsi que des estimations d’erreur peuvent se trouver dans la cette référence [1.3].
Pour calculer la solution de notre problème réduit on va donc calculer \(\hat{\bs{x}}^n\) solution de l’équation (17.46). Pour cela on utilise en général une méthode de Newton initialisée avec \(\hat{\bs{x}}^{n-1}\text{.}\) L’évaluation du résidu nécessite d’évaluer la fonction nonlinéaire \(\bs{f}\) et donc sera coûteuse. Comme pour les méthodes LSG il faudra trouver une solution pour rapidement évaluer la nonlinearité.

Conclusion.

Dans ce section, on a défini comment construire un décodeur et un encodeur (réduction), comment construire le modèle réduit (méthodes LSG et LSPG). Cependant il faut actuellement repasser en grande dimension pour calculer la non-linéarité. On verra dans la suite comment faire.
PrevTopNext