Projection-free approximation of flows of harmonic maps with quadratic constraint accuracy and variable step sizes

1. Problématique et Contexte

Le papier traite de l'approximation numérique des écoulements de cartes harmoniques vers une sphère (le flot de chaleur des cartes harmoniques). Ce problème mathématique consiste à trouver des champs vectoriels $u$ minimisant l'énergie de Dirichlet sous la contrainte de longueur unitaire $|u|=1$ presque partout.

Les équations d'Euler-Lagrange associées sont non linéaires et incluent un multiplicateur de Lagrange. Une approche courante pour résoudre ces équations est d'utiliser des flots de gradient (comme le flot de chaleur $L^2$ ou $H^1$).

Défis numériques principaux :

• Préservation de la contrainte : La solution numérique doit respecter la condition $|u|=1$. Les méthodes classiques nécessitent souvent des étapes de projection (renormalisation) qui peuvent imposer des restrictions sur la taille du pas de temps ou la triangulation.
• Stabilité énergétique : Il est crucial que la méthode numérique préserve la décroissance de l'énergie (stabilité inconditionnelle).
• Précision de la contrainte : L'erreur de violation de la contrainte (écart à $|u|=1$) doit être contrôlée avec une haute précision.
• Adaptabilité : La capacité à gérer des pas de temps variables est essentielle pour accélérer la convergence vers les points stationnaires ou pour affiner la résolution près des singularités (explosions finies ou données initiales non lisses).

2. Méthodologie Proposée

Les auteurs développent une méthode linéairement implicite et sans projection (projection-free) basée sur la formulation faible d'Alouges.

A. Le schéma $(\theta, \mu)$

L'article introduit une famille de méthodes paramétrées par $\theta$ et $\mu$ ($0 < \theta \le 1$, $0 \le \mu \le 1$) pour la discrétisation temporelle.

• Formulation : À chaque étape $n$, on cherche une différence arrière $d_t u^n$ telle que :
  1. Elle satisfait une condition de longueur unitaire linéarisée : $d_t u^n \cdot \hat{u}_{n-1+\mu} = 0$.
  2. Elle satisfait l'équation discrétisée : $(d_t u^n, v)_\star + (\nabla u_{n-1+\theta}, \nabla v) = 0$ pour tout $v$ orthogonal à $\hat{u}_{n-1+\mu}$.
• Paramètres clés :
  • $\hat{u}_{n-1+\theta}$ est une interpolation linéaire (valeur à $t_{n-1} + \theta \tau$).
  • $\hat{u}_{n-1+\mu}$ est une extrapolation linéaire basée sur les deux pas précédents (valeur à $t_{n-1} + \mu \tau$).
• Cas particuliers :
  • $\theta=1, \mu=0$ : Méthode d'Euler implicite linéaire (ordre 1 pour la contrainte).
  • $\theta=1/2, \mu=1/2$ : Méthode du point milieu (proposée comme la méthode optimale).
  • $\theta=1, \mu=1/2$ : Méthode d'Euler modifiée.

B. Extension aux pas de temps variables

Une contribution majeure est l'extension de cette méthode à des pas de temps variables ($\tau_n$).

• Les auteurs dérivent des lois d'énergie discrètes et des bornes d'erreur de contrainte pour le cas non uniforme.
• Ils proposent une stratégie de mise à jour du pas de temps (basée sur le rapport de pas successifs $s_n$) pour équilibrer les termes d'erreur, permettant d'accélérer la convergence vers l'état stationnaire tout en maintenant la stabilité.
• Un mécanisme de contrôle adaptatif est également introduit : le pas de temps est réduit si la dynamique du flot s'accélère (dérivée temporelle croissante) et augmenté sinon.

3. Contributions Clés et Résultats Théoriques

1. Stabilité Énergétique Inconditionnelle :
   La méthode est prouvée comme étant inconditionnellement stable en énergie pour tout choix de $\theta \ge 1/2$. Une identité d'énergie discrète est établie, garantissant la décroissance de l'énergie de Dirichlet à chaque itération.

2. Précision Quadratique de la Contrainte :
   Contrairement à la méthode d'Euler implicite standard (erreur de contrainte linéaire en $\tau$), la méthode du point milieu ($\theta=\mu=1/2$) atteint une précision quadratique ($O(\tau^2)$) pour la violation de la contrainte $|u|^2 - 1$, sous une condition de régularité discrète faible (bornitude des différences secondes).

   • Cela signifie que l'erreur de contrainte est contrôlée indépendamment du nombre d'itérations et ne dépend que de la taille du pas de temps.

3. Gestion des Pas de Temps Variables :
   Le papier démontre que la stabilité énergétique et le contrôle de la violation de la contrainte sont préservés même avec des pas de temps variables. Cela permet d'utiliser des techniques d'accélération (augmenter $\tau$ quand la solution se stabilise) sans compromettre la robustesse du schéma.

4. Comparaison avec les méthodes existantes :
   La méthode proposée est comparée à la méthode d'Euler implicite linéaire et à la méthode BDF2 (Backward Differentiation Formula d'ordre 2).

   • La méthode BDF2 nécessite une stabilité "G-stable" complexe qui rend son extension aux pas de temps variables difficile.
   • La méthode $(\theta, \mu)$ offre une alternative plus flexible et simple à implémenter avec des performances similaires ou supérieures.

4. Résultats Numériques

Des expériences numériques en 2D (éléments finis linéaires) valident la théorie :

• Convergence :
  • Pour le flot de gradient $H^1$, la méthode du point milieu montre une convergence d'ordre 2 pour l'erreur de contrainte ($L^1$ et $L^\infty$), confirmant la théorie.
  • Pour le flot de gradient $L^2$, la convergence d'ordre 2 est observée, mais avec une phase pré-asymptotique plus longue.
• Précision de la contrainte :
  • La méthode du point milieu est nettement plus précise que la méthode BDF2 (erreur environ 3 fois plus faible pour le même pas de temps).
  • La méthode d'Euler modifiée ($\theta=1, \mu=1/2$) se comporte comme la méthode du point milieu en termes de précision de la contrainte.
• Singularités :
  • Dans le cas de données initiales générant des singularités (explosion de la dérivée), la méthode conserve sa stabilité. L'erreur de contrainte converge, bien que l'ordre de convergence puisse être légèrement inférieur à 2 en raison de la perte de régularité de la solution exacte.
• Pas de temps variables :
  • L'utilisation de pas de temps variables (augmentant au cours du temps) permet de réduire considérablement le nombre d'itérations nécessaires pour atteindre un état stationnaire (gain de ~20% dans les tests) tout en maintenant une haute précision, en particulier pour le flot $H^1$.
  • Le contrôle adaptatif permet de réduire automatiquement le pas de temps lors des phases de dynamique rapide (près des singularités) et de l'augmenter ensuite.

5. Signification et Impact

Ce travail est significatif pour plusieurs raisons :

• Efficacité Algorithmique : Il propose une méthode qui évite les projections coûteuses (renormalisation) tout en garantissant une précision élevée de la contrainte géométrique.
• Flexibilité : La capacité à utiliser des pas de temps variables sans perdre la stabilité inconditionnelle est une avancée majeure par rapport aux méthodes BDF classiques pour ce type de problème géométrique.
• Applicabilité : Bien que développé pour les cartes harmoniques, le cadre théorique s'applique à une large classe de problèmes géométriquement contraints, notamment en magnétisme (équation LLG), en théorie des cristaux liquides et en théorie de la flexion.
• Robustesse : La méthode est robuste face aux singularités et aux données initiales non lisses, ce qui est crucial pour les simulations physiques réalistes.

En résumé, les auteurs ont conçu un schéma numérique linéairement implicite, sans projection, inconditionnellement stable et d'ordre 2 pour la contrainte, capable de s'adapter dynamiquement à la complexité de la solution, offrant ainsi un outil supérieur pour la simulation des écoulements de cartes harmoniques.