Projet Aladin

Précédent : Présentation générale et objectifs Remonter : Rapport activite 1997 Suivant : Grands domaines d'application

• Équations différentielles ordinaires et/ou algébriques
  • Systèmes différentiels ordinaires
  • Équations algébro-différentielles
  • Systèmes hamiltoniens
• Problèmes linéaires et non linéaires
  • Méthodes directes pour les sytèmes linéaires
  • Méthodes itératives pour les systèmes linéaires
  • Accélération de convergence
  • Problèmes non linéaires
• Problèmes aux valeurs propres
  • Méthodes de Davidson
  • Accélération de la méthode d'Arnoldi
  • Pseudo-spectres et dichotomie spectrale
• Encadrement garanti de résultats

Fondements scientifiques

Équations différentielles ordinaires et/ou algébriques

  Mots-clés : systèmes différentiels, systèmes algébro-différentiels, indice 2, systèmes hamiltoniens, méthodes générales linéaires, méthodes symplectiques

Résumé : On cherche à résoudre numériquement des problèmes de valeur initiale pour des systèmes d'équations différentielles ordinaires ou avec contraintes, c'est-à-dire de la forme $$\begin{eqnarray*} \left\{ \begin{array} {ccc} y'(x)&=& f(y(x)), \\ y(x_0) &=& y_... ...), & y(x_0)=y_0 \\ 0 &=& g(y(x)), & z(x_0)=z_0.\end{array}\right.\end{eqnarray*}$$

Dans le cas des équations différentielles ordinaires se pose un problème de coût calcul pour les systèmes de grande taille. Le recours au parallélisme semble alors incontournable, mais il se heurte au caractère intrinsèquement séquentiel des méthodes numériques usuelles. Nos travaux se sont orientés dans des directions orthogonales, suivant que l'on cherche à paralléliser <<à travers l'intervalle d'intégration>> ou <<à travers la méthode>>. Les équations algébro-différentielles posent quant à elles des problèmes spécifiques : en particulier, l'ordre de convergence habituel se trouve notablement diminué du fait de la présence de contraintes algébriques. Nous nous sommes donc attachés à compenser ce phénomène de <<réduction d'ordre>> par des techniques de smoothing. Enfin nous nous intéressons aux systèmes hamiltoniens, qui se présentent sous forme d'équations différentielles ordinaires mais méritent souvent un traitement spécifique. Celui-ci est destiné à préserver certains invariants géométriques qui leur sont attachés.

Systèmes différentiels ordinaires

Ces sont des systèmes qui se posent sous la forme [HNW93,HW96] $$\begin{eqnarray*} \left\{ \begin{array} {ccc} y'(x)&=& f(y(x)), \\ y(x_0) &=& y_0.\end{array}\right.\end{eqnarray*}$$ $y$ et $f$ sont ici des vecteurs de $R^m$, issus typiquement de la discrétisation d'une équation aux dérivées partielles. Dans une telle situation, $m$ peut être très grand et la probabilité que le problème soit <<raide>> très élevée. En clair, il est fort probable que la résolution numérique se heurte à des problèmes de stabilité et qu'il soit nécessaire de recourir à une méthode implicite, donc d'un coût calcul potentiellement prohibitif. Face à cette situation, l'usage des méthodes de différentiation rétrograde s'est généralisé, en raison de leur quasi-optimalité en terme de coût par pas. Il reste que ces méthodes souffrent d'un déficit de stabilité, auquel les codes courants tels DASSL [BCP89] ou VODE remédient par l'utilisation de méthodes d'ordre inférieur à $2$. Les méthodes ``Singly-Implicit Runge-Kutta'' s'affranchissent partiellement de cette difficulté. La propriété fondamentale des méthodes SIRK[Bur78] est la suivante : la matrice de coefficients $A$ possède une seule valeur propre $\lambda$ de multiplicité $s$. Ainsi, si $J$ est la jacobienne $m \times m$ du système différentiel, la décomposition LU (lower-upper) de la matrice $(I_s \otimes I_m - h A \otimes J)$, dont le coût est prédominant dans les formules de passage d'un pas au suivant, peut être évitée et remplacée par la décomposition LU de la matrice $I_m-h \lambda J$. Modulo quelques transformations linéaires, le coût de ces méthodes est alors ramené à un niveau comparable à celui des méthodes multipas et elles sont parfaitement stables.

Il est également possible de contourner cette difficulté en ayant recours au parallélisme : une première technique consiste à reformuler le système sous une forme itérative. Considérons pour simplifier la forme de Picard

$$\begin{eqnarray*} \left\{ \begin{array} {ccc} {y^{(k+1)}}' (x)&=& f(y^{(k)}(x)), \\ y^{(k+1)}(x_0) &=& y_0.\end{array}\right.\end{eqnarray*}$$ Sous cette forme, l'intégration se résume à un problème de quadrature et ne présente plus aucune difficulté. Il est même possible de découpler le système par blocs de composantes et d'affecter chaque bloc à un processeur. On voit aisément que le succès de cette approche est conditionné par la vitesse de convergence du processus, ce qui exclut d'entrée la forme de Picard. Cependant des variantes existent, du type Jacobi, donc implicites, mais néanmoins susceptibles d'être parallélisées. Nous nous sommes quant à nous intéressés à une idée similaire, issue des travaux de A. Bellen et M. Zennaro [BZ89], et à son application aux problèmes dissipatifs [2] et plus récemment aux problèmes de trajectoires de satellite [28].

Une seconde technique consiste à concevoir des méthodes intrinsèquement parallèles possédant des domaines de stabilité plus larges. C'est le cas des méthodes DIMSIM[But93] implicites, dont le format est le suivant

$$\begin{eqnarray*} Y_i &=& h \sum_{j=1}^{s} a_{i,j} f(Y_i) + \sum_{j=1}^{r} u_{i,... ... b_{i,j} f(Y_i) + \sum_{j=1}^{r} v_{i,j} y_j^{[n]}, i=1,\ldots,r.\end{eqnarray*}$$ Les vecteurs $y_i^{[n]}$ et $Y_i$ désignent des approximations de la solution exacte ou de quantités relatives à la solution exacte dont la définition précise importe peu. Toute méthode numérique <<classique>> peut s'écrire sous ce format et son coût par pas est essentiellement déterminé par la forme de la matrice $A$. La première des deux équations ci-dessus constitue en effet un système implicite non-linéaire dont la résolution nécessite là encore la décomposition LU de la matrice $$\begin{eqnarray*} (I_s \otimes I_m - h A \otimes J).\end{eqnarray*}$$ Dans le cas d'une méthode de Runge-Kutta classique (c'est-à-dire non SIRK) à $s$ étapes internes, $A$ est une matrice pleine et le système de dimension $s \times m$. Cette dimension n'est que de $m$ pour les méthodes de différentiation rétrograde. Les méthodes DIMSIM que nous avons considérées possèdent une matrice $A$ diagonale. Le système est alors découplable en $s$ sous-systèmes indépendants de dimension $m$. Le grand avantage de ces méthodes est qu'il est possible de construire des méthodes d'ordre élevé et parfaitement stables [3].

Équations algébro-différentielles

Ce sont des équations du type

$$\begin{eqnarray*} \left\{ \begin{array} {cccc} y'(x) &=& f(y(x),z(x)), & y(x_0)=y_0, \\ 0 &=& g(y(x)), & z(x_0) = z_0,\end{array}\right.\end{eqnarray*}$$ où l'on suppose que $g_yf_z$ est inversible dans un voisinage de la solution exacte $(y(x),z(x))$. Cette hypothèse assure qu'il est possible de retransformer le système algébro-différentiel en un système différentiel pur, et ce en dérivant la contrainte $0=g(y(x))$. En omettant la dépendance en $x$, il vient en effet successivement $$\begin{eqnarray*} 0 &=& g_y(y)y' = g_y(y) f(y,z), \\ 0 &=& g_{yy}(y)(f(y,z),f(y,z)) + g_y(y)f_y(y)f(y,z) + g_y(y)f_z(y,z)z',\end{eqnarray*}$$ d'où il est possible de tirer d'après l'hypothèse faite sur $g_yf_z$ une expression de $z'$. On voit au passage que deux différentiations ont été nécessaires, ce qui caractérise les systèmes de Hessenberg d'indice $2$. Les problèmes d'indice inférieur à $2$ ne posent pas de difficultés particulières, alors que le traitement direct des systèmes d'indice supérieur à $3$ est considéré comme périlleux. Ceci confère aux systèmes d'indice $2$ une importance particulière, reflétée par l'abondance de la littérature portant sur ce cas[HLR89,HW96]. Comme précédemment indiqué, une méthode de Runge-Kutta appliquée à un système d'indice $2$, subit une réduction de son ordre de convergence. Si la méthode possède $s$ étapes internes, l'ordre des méthodes de Runge-Kutta de type Radau IIA, qui sont les mieux adaptées à la situation décrite, est de $2s-1$ pour une équation ordinaire. Il n'est plus que de $s$ pour la composante algébrique ($z$) d'un système d'indice $2$. C'est ce qu'on observe couramment lorsqu'on utilise le code RADAU5 fondé sur ces méthodes.

Systèmes hamiltoniens

Ce sont des systèmes de la forme[SSC94] $$\begin{eqnarray*} \left\{ \begin{array} {cccc} p'(x) &=& -\frac{\partial H}{\par... ...\frac{\partial H}{\partial p}(x), & q(x_0)=q_0,\end{array}\right.\end{eqnarray*}$$ où $H(p,q)$ est une fonction scalaire dite hamiltonienne et $p$ et $q$ des vecteurs de $R^n$. De tels systèmes possèdent un certain nombre de propriétés remarquables, au rang desquelles la conservation de $H$ et de la différentielle $2$-forme $$\begin{eqnarray*} \omega^2 = \sum_{i=1}^{n} dp_i \wedge dq_i\end{eqnarray*}$$ le long des solutions, exprimant en dimension $2$ la conservation des surfaces et en dimension supérieure d'une quantité similaire quoique plus abstraite. On peut également caractériser ces systèmes en exprimant la symplecticité de la fonction flot $\Phi_{x_0,x_1}(p,q)$ qui à un point $(p_0,q_0)$ associe le point $(p(x_1),q(x_1))$ solution du système hamiltonien décrit ci-dessus. $\Phi$ est dite symplectique si et seulement si $$\begin{eqnarray*} {\Phi '}^T \left[ \begin{array} {cc} 0 & I_n \\ -I_n & 0\end{a... ... \left[ \begin{array} {cc} 0 & I_n \\ -I_n & 0\end{array}\right].\end{eqnarray*}$$ Les méthodes dites symplectiques ont été conçues pour conserver certaines quantités de nature géométrique, telle $\omega^2$. Si $\varphi_{x_0,x_0+h}$ désigne la fonction flot numérique associée à une méthode symplectique et $h$le pas d'intégration, alors, quelle que soit la valeur de $h$, $\varphi$vérifie elle-aussi $$\begin{eqnarray*} {\varphi '}^T \left[ \begin{array} {cc} 0 & I_n \\ -I_n & 0\en... ... \left[ \begin{array} {cc} 0 & I_n \\ -I_n & 0\end{array}\right].\end{eqnarray*}$$ Les méthodes de Runge-Kutta symplectiques ont certaines propriétés très séduisantes. Ainsi, si $X$ désigne la longueur de l'intervalle d'intégration et $p$ l'ordre de la méthode considérée, l'erreur $\Delta$entre solution approchée et solution exacte d'un problème hamiltonien, par exemple périodique pour simplifier, varie linéairement par rapport à $X$ : $$\begin{eqnarray*} \Delta \approx C h^p X.\end{eqnarray*}$$ Notons que pour une méthode non symplectique, on a $$\begin{eqnarray*} \Delta \approx C h^p X^2.\end{eqnarray*}$$ Pour des intervalles d'intégration <<astronomiques>> cela représente un avantage considérable.

Problèmes linéaires et non linéaires

  Mots-clés : matrice symétrique, matrice creuse, espace de Krylov, linéarisation, itératif, préconditionnement, déflation

matrice symétrique: les matrices symétriques, telles que $A=A ^t$, sont très fréquentes. Grâce à leurs propriétés, notamment spectrales, la résolution de systèmes symétriques est simplifiée.

matrice creuse: les matrices de très grande taille contiennent en général un grand nombre de coefficients nuls. Lorsque les coefficients non nuls sont à des positions précises, on dit que la matrice est creuse et structurée. Sinon elle est creuse et générale.

espace de Krylov: l'espace engendré par $\{v,Av,\ldots,A ^{m-1}v\}$ est un espace de Krylov. Projeter le problème linéaire sur ce sous-espace permet de se ramener à un problème de petite taille qui approche le problème initial.

linéarisation: en remplaçant localement la fonction $F$ par l'espace tangent, et en résolvant un problème linéaire associé, on trouve une approximation de la solution au problème $F(u)=0$.

Résumé : Un problème linéaire est défini par une matrice $A \in { \mathbb R}^{N \times N}$ et un vecteur $b \in { \mathbb R}^N$ ; on cherche $x \in { \mathbb R}^N$ tel que $$Ax=b.$$ L'entier $N$est l'ordre ou la taille de la matrice.

Ce problème est au coeur de nombreuses applications scientifiques : discrétisation d'équations aux dérivées partielles, linéarisation de problèmes non linéaires, traitement d'images, etc.

Un problème non linéaire est défini par une fonction $F$ d'un domaine $D$ de ${ \mathbb R}^N$ dans ${ \mathbb R}^N$ ; on cherche à résoudre l'équation

$$u \in { \mathbb R}^N,~F(u)=0.$$ Ce problème, très général, se pose après discrétisation d'équations différentielles par des schémas implicites (voir module ), après discrétisation d'équations aux dérivées partielles non linéaires, etc.

Le projet étudie principalement les méthodes itératives dites de Krylov pour résoudre les systèmes non linéaires. La recherche porte sur des méthodes permettant d'accélérer la convergence. Dans le cas non linéaire, le projet étudie les méthodes itératives de linéarisation et leur couplage avec la résolution des problèmes linéaires induits.

Méthodes directes pour les sytèmes linéaires

Si $N$ est assez petit ($N \leq 5000$ environ), le système $Ax=b$ se résout par une méthode directe basée sur la factorisation de Gauss avec pivot $PA=LU$, où $P$ est une matrice de permutation liée à la stratégie de pivot qui assure la stabilité numérique, $L$ est triangulaire inférieure et $U$ est triangulaire supérieure. Cette méthode est précise et stable numériquement mais sa complexité, mesurée par $O(N^3)$ opérations flottantes et $O(N^2)$ variables flottantes en mémoire, est un frein à son utilisation pour $N$ grand.

Pour réduire le coût mémoire, il faut alors exploiter la structure creuse de la matrice en ne stockant que peu, voire pas, de coefficients nuls [DER86]. Il existe de nombreux types de stockage creux : bande, profil, compressé par lignes, par colonnes, par diagonales, etc. Mais la factorisation de Gauss induit du remplissage dans les facteurs $L$ et $U$. Des techniques de renumérotation ont pour objectif de limiter ce remplissage : minimisation de la largeur de bande ou du profil, degré minimum, etc. Dans le cas symétrique où $A=L L^T$,la complexité devient alors $O(N \times d^2)$ opérations flottantes et $O(NZ(L))$ variables flottantes, où $NZ(L)$ est le nombre de coefficients non nuls dans le facteur $L$ et $d=NZ(L)/N$ est le nombre moyen par ligne. Typiquement, $NZ(L)=O(N^{3/2})$ pour une discrétisation de problèmes 2D et $NZ(L)=O(N^{5/3})$ pour des problèmes 3D. Une autre difficulté est de concevoir une version parallèle. Les méthodes multifrontales sont souvent efficaces. Il est aussi possible d'exploiter la structure creuse pour dégager un parallélisme intrinsèque à gros grain [6].

Méthodes itératives pour les systèmes linéaires

  Pour $N$ très grand ($N \geq 10000$ environ), le volume mémoire des méthodes directes est souvent prohibitif. De plus en plus, elles sont remplacées par des méthodes itératives. Les méthodes stationnaires de type relaxation ont l'inconvénient de converger lentement et pour seulement certaines classes de matrices. Les méthodes de projection sont plus générales et plus robustes.

De plus, ces méthodes sont intéressantes parce que la matrice $A$ n'est utilisée qu'à travers l'opérateur produit matrice-vecteur $w=Av$. Il est donc possible d'utiliser un stockage compressé, voire de ne pas stocker $A$.Un autre avantage est l'utilisation de méthodes "matrix-free ", où le produit matrice-vecteur est recalculé ou approché à chaque occurence. De plus, cette opération se parallélise assez bien.

Dans cette famille, les méthodes de projection de Krylov sont les plus étudiées actuellement [Saa95b]. Elles sont définies par le sous-espace de Krylov

$${\cal K}_m(A,r_0) = vect\{r_0,Ar_0,\ldots,A ^{m-1}r_0\},$$ construit itérativement, par une matrice $B$ souvent symétrique définie positive et par deux conditions : la condition d'espace $$x_{m+1}-x_m \in K_m$$ et la condition de Petrov-Galerkin $$(B(x-x_m),v)=0, ~~\forall v \in {\cal K}_m.$$ Ces méthodes de Krylov sont polynomiales, en effet $$x-x_m=R_m(A)(x-x_0).$$ On en déduit qu'elles convergent en au plus $N$ itérations mais l'objectif est d'obtenir une bonne approximation en beaucoup moins d'itérations. Si $A$ est diagonalisable sous la forme $A=XDX^{-1}$ avec $D$ diagonale contenant les valeurs propres de $A$, alors $R_m(A)=XR_m(D)X ^{-1}$ et il suffit d'étudier $R_m(D)$.Cette propriété permet de relier les méthodes de Krylov aux outils mathématiques manipulant les polynômes.

Il faut distinguer, comme pour le cas des méthodes directes, les systèmes symétriques et non symétriques.

Les systèmes symétriques sont les plus faciles à résoudre. Les méthodes de Krylov peuvent dans ce cas construire ${\cal K}_m$ et calculer $x_{m+1}$ à l'aide de récurrences courtes, d'où une faible complexité. Si de plus $A$ est définie positive, la méthode du gradient conjugué permet d'allier récurrences courtes et minimisation :

$$\Vert x-x_m\Vert _A \leq \Vert x-x_0+v\Vert _A,~\forall v \in {\cal K}_m.$$ Les systèmes non symétriques sont plus difficiles à résoudre. Une approche possible est de se ramener au cas symétrique défini positif en résolvant l'équation normale $A^t A x= A ^t b$ ou l'équation $A A ^t (A ^{-t}x)=b$. Cette solution est robuste mais coûteuse puisque chaque itération requiert à la fois le produit par $A$ et par $A ^t$.Hormis cette méthode, il existe deux grandes classes de méthodes, soit avec récurrences courtes, soit avec minimisation, les deux propriétés étant incompatibles ici. La méthode GMRES (Generalized Minimum Residual), qui est très utilisée pour sa robustesse et son efficacité, impose une propriété de minimisation $$\Vert r_m\Vert _2=\min_{v \in {\cal K}_m}{\Vert r_0-Av\Vert _2}$$ mais la construction de ${\cal K}_m$ a une complexité en $O(mNZ(A)+m^2N)$ opérations flottantes et $O(NZ(A)+mN)$ variables flottantes. Un moyen de limiter ce coût est de redémarrer l'algorithme toutes les $m$ itérations, toutefois la convergence n'est plus garantie. Le choix du paramètre $m$ s'avère très délicat.

Accélération de convergence

  Pour le gradient conjugué comme pour GMRES, la vitesse de convergence dépend des valeurs propres de la matrice (le spectre). Préconditionner la matrice, c'est-à-dire résoudre $$M_1 A M_2 (M_2 ^{-1} x) = M_1 b,$$ avec $M_1$ et $M_2$ inversibles, permet d'accélérer la convergence grâce à un spectre de $M_1 A M_2$ plus favorable [Bru95]. Chaque itération est alors plus coûteuse puisqu'elle implique, outre le produit $Av$, les produits $M_1v$ et $M_2v$ (il est bien sûr hors de question de stocker la matrice pleine $M_1 A M_2$).

Le préconditionnement diagonal $M_1=D, M_2=I$ a un coût faible en $O(N)$opérations flottantes et variables flottantes. La parallélisation en est aisée mais l'efficacité est parfois réduite.

La factorisation incomplète est définie par

$$M_2=I,~M_1=(L_1 U_1)^{-1}~~{\mbox avec}~~A=L_1 U_1 + R$$ et $R$ est choisi implicitement par le taux de remplissage dans $L_1$ et $U_1$. Ce préconditionnement est en général efficace mais est en contrepartie coûteux et peu parallèle (on retombe sur les inconvénients des méthodes directes).

Les préconditionnements polynomiaux sont également assez coûteux bien qu'ils soient parallèles puisque seul le produit $w=Av$ intervient.

Une autre approche est un préconditionnement par déflation, défini par

$$M_2=I ~~et~~ M_1=I-U U ^t + U (U ^t A U)^{-1} U ^t$$ où $U$ est une base orthonormée d'un sous-espace invariant (en pratique une approximation). Des variantes peuvent être construites autour de la même idée. Si $U$ est exactement invariant, la matrice préconditionnée restreinte à $vect(U)$ est l'identité et son spectre sur $vect(U) ^{\perp}$ est le spectre de $A$ privé des valeurs propres associées à $U$ de sorte que le problème est ramené à une résolution sur $vect(U) ^{\perp}$. Il faut donc choisir $U$ associé aux petites valeurs propres qui freinent la convergence. L'approximation de $U$ se fait grâce à des relations avec les méthodes de Lanczos et d'Arnoldi (voir module ). Le produit $w=M_1 v$ est basé sur des opérations de type BLAS2 qui se parallélisent bien.

Augmenter l'espace des solutions avec un sous-espace de Krylov calculé précédemment ou une approximation d'un sous-espace invariant est une autre accélération possible, assez similaire à la précédente. Leur comparaison est faite dans [Saa95a] et dans [5].

Enfin les méthodes par blocs ont un effet semblable au précédent et ont l'intérêt de recourir à des opérateurs de type BLAS2 ou BLAS3.

Problèmes non linéaires

  La méthode de Newton est souvent utilisée pour résoudre le problème $F(u)=0$. C'est une méthode de linéarisation qui s'écrit $$u_{k+1}=u_k - J(u_k) ^{-1} F(u_k)$$ où $J(u_k)$ est le Jacobien de $F$ en $u_k$. Les méthodes de type Newton modifié utilisent un Jacobien approché et les méthodes de type Newton-inexact résolvent de façon approchée le système linéaire $J(u_k) x_k = -F(u_k)$[OR70],[DS83].

Alors que la convergence locale de Newton est quadratique, donc très rapide, celles de Newton modifié et Newton inexact sont linéaires, au mieux super-linéaires donc plus lentes. Par contre, la complexité par itération est en général moindre, ce qui peut compenser le plus grand nombre d'itérations.

Les méthodes de Newton-Krylov, de type inexact, résolvent le système linéaire à l'aide d'une méthode de projection sur un espace de Krylov (voir section ). La difficulté réside dans un bon choix des critères de convergence et dans un préconditionnement efficace (voir section ).

D'autre part, pour garantir la convergence globale, il faut combiner la méthode de Newton avec une technique de ``backtracking'' ou de continuation par exemple.

Problèmes aux valeurs propres

  Mots-clés : valeurs propres,Davidson,Chebychev,Faber,pseudo-spectres,dichotomie spectrale, intervalles


Davidson: Méthode de calcul de valeurs propres adaptée aux grandes matrices creuses symétriques qui peut être vue comme une méthode de Lanczos accélérée à chaque étape par un pas d'une méthode de Newton inexacte pour corriger l'estimation courante du vecteur propre. Elle nécessite une matrice de préconditionnement. La méthode peut se décliner en méthode de Davidson ou de Jacobi-Davidson.
Chebychev, Faber: Deux types de polynômes définis sur un domaine $D$ convexe du plan complexe (respectivement ellipse et polygone) bornés par l'unité sur $D$ et qui croissent très rapidement dès que l'on s'éloigne du contour vers l'extérieur de $D$.
Pseudo-spectre: Ensemble des valeurs propres de toutes les matrices voisines d'une matrice donnée où le voisinage est défini par la norme euclidienne. Permet de connaître la sensibilité aux perturbations des valeurs propres d'une matrice.
Dichotomie spectrale: Procédé de partitionnement de l'ensemble des valeurs propres d'une matrice qui contient aussi l'évaluation d'un critère de fiabilité sur le résultat du partitionnement. Contient généralement aussi le calcul de projecteurs spectraux (projecteurs définis par la diagonalisation de la matrice).

Résumé : Le problème standard de valeur propre consiste pour une matrice $A \in { \mathbb R}^{N\times N}$ donnée à trouver tous les couples $(\lambda , x) \in { \mathbb R} \times { \mathbb R}^N$ (ou seulement une partie d'entre-eux) qui vérifient : $$Ax= \lambda x .$$ Le problème généralisé est défini par deux matrices $A,~B \in { \mathbb R}^{N\times N}$ et les couples doivent alors vérifier : $$Ax= \lambda B x .$$ L'amélioration des méthodes de calcul de valeurs propres porte essentiellement sur la recherche d'accélérateur de la convergence à appliquer dans les méthodes itératives adaptées aux grandes matrices creuses. Le projet travaille depuis plusieurs années sur la méthode de Davidson pour le cas des matrices symétriques, méthode qu'il adapte maintenant au calcul des plus petites valeurs singulières de matrices. Dans le cas des matrices non symétriques (ou non hermitiennes), les travaux ont porté sur la recherche d'accélérateur polynômiaux (polynômes de Chebychev et de Faber).

La mesure de la sensibilité des valeurs propres d'une matrice aux perturbations est un problème à résoudre lorsque l'on doit localiser les valeurs propres d'une matrice imprécisément connue. On peut définir l'ensemble des perturbations possibles par la norme euclidienne des matrices ou par des matrices d'intervalles.

Méthodes de Davidson

On se place dans le cas du calcul de quelques valeurs propres extrêmes d'une grande matrice symétrique $A$ (ou hermitienne complexe). La méthode de Davidson est une méthode de sous-espace car elle génère pas-à-pas un système orthonormé de vecteurs $V_m$ sur lequel on projette le problème initial afin d'obtenir des approximations des valeurs et vecteurs propres cherchés. Soit $(\lambda , z)$ un couple d'éléments propres de la matrice $H_m={V_m}^tAV_m$, alors $(\lambda , x=V_mz)$ est une approximation d'un couple d'éléments propres de $A$. Contrairement à la méthode de référence Lanczos, les sous-espaces engendrés ne sont pas des sous-espaces de Krylov (voir la définition dans ), car on incorpore à chaque pas une correction du vecteur $x$ par un procédé de Newton : on recherche $y$ petit tel que $y \bot x$ et $x+y$ soit un vecteur propre de $A$. En négligeant les termes du deuxième ordre par rapport à $\Vert y \Vert$, on est ramené à résoudre $$\begin{eqnarray*} r &=& (\lambda I - A) y \\ \mbox{où } r &=& Ax-\lambda x \mbox{ et } y \bot x ~\end{eqnarray*}$$ Les méthodes de Davidson consistent à résoudre approximativement la première équation, soit en remplaçant $A$ par une matrice $M$ plus facile à manipuler, soit en effectuant quelques pas d'une méthode de type gradient conjugué. Lorsque l'on résout directement dans l'orthogonal de $x$, on obtient une méthode de type Jacobi-Davidson. Le comportement de la méthode de Davidson a été étudié dans [4] tandis que la méthode de Jacobi-Davidson est décrite dans [SVdV96]. Ces méthodes sont de nettes améliorations de la méthode de Lanczos dans le cas du calcul de valeurs propres de petites valeurs absolues. C'est pourquoi le projet les a considérées pour calculer les plus petites valeurs singulières d'une grande matrice en les appliquant à la matrice ${A}^tA$. L'étude théorique en est faite dans [7] et l'application au calcul de pseudo-spectres dans [22].

Accélération de la méthode d'Arnoldi

Nous étudions le cas d'une matrice $A$ de grande taille réelle non symétrique ou complexe non hermitienne, dont par conséquent les valeurs propres peuvent être complexes. La méthode d'Arnoldi [Cha88,Saa92] construit pas à pas une base orthonormée $V_m$ de l'espace de Krylov ${\cal K}_m(A,v)$ et approche les couples propres cherchés par ceux de la matrice $H_m={V_m}^tAV_m$ qui est une matrice de Hessenberg de petite taille. Appliquée à une matrice symétrique, cette méthode serait la méthode de Lanczos. Comme la dimension $m$ de la base est limitée par les capacités de stockage des données, il est nécessaire de redémarrer périodiquement le procédé avec la meilleure approximation courante du vecteur propre cherché et c'est à cet endroit que s'insèrent les techniques d'accélération. Afin d'accoître les composantes du nouveau vecteur initial $v$ dans les directions des vecteurs propres que l'on recherche, on calcule $p(A)v$ où $p$ est un polynôme dont le module discrimine les valeurs propres cherchées des autres. Il s'agit donc de définir une courbe convexe qui entoure les valeurs propres à rejeter et elles seulement. Pour cela on choisit généralement une ellipse et le polynôme de Chebychev qu'elle définit. Si la matrice est réelle, l'ellipse est choisie symétrique par rapport à l'axe des réels. Dans le cas complexe, le placement de l'ellipse est plus délicat. Il est souvent préférable de construire un polygone et de calculer le polynôme de Faber associé, polynôme obtenu à partir de la transformation de Schwarz-Christopher qui fait correspondre l'extérieur du polygone avec l'extérieur du disque unité. Le projet a travaillé ces trois approches et les a appliquées à la méthode d'Arnoldi par blocs [11,23]. Nous avons systématiquement choisi dans ce domaine d'exprimer les algorithmes par des versions par blocs car elles permettent de meilleures efficacités informatique et numérique.

Pseudo-spectres et dichotomie spectrale

Dans de nombreuses applications, les valeurs propres d'une matrice sont calculées afin de décider si tout le spectre est ou non inclus dans une partie donnée du plan complexe (demi-plan des complexes à partie réelle négative, disque unité, etc.) Malheureusement, la matrice n'est définie qu'à une précision donnée (au maximum la précision de l'arithmétique flottante) et le conditionnement (voir la définition en ) de certaines valeurs propres peut être assez élevé pour que le résultat ne soit pas sûr. On désire donc non calculer les valeurs propres d'une unique matrice mais plutôt en déterminer les zones de variation lorsque la matrice varie dans un voisinage donné de la matrice initiale. On considère ici le cas où le voisinage est défini par la norme euclidienne mais il peut aussi être défini par une matrice d'intervalles (cf ).

La notion de pseudo-spectre a été introduite simultanément mais de manière indépendante par Godunov [GKK91] et Trefethen [Tre91]. Pour une matrice $A$ et un paramètre de précision $\epsilon$ donnés, le pseudo-spectre $\Lambda _{\epsilon} (A)$ est l'ensemble des valeurs propres de toutes les matrices $A+\Delta$ où $\Vert \Delta \Vert \le \epsilon \Vert A \Vert$. Sa caractérisation repose sur l'équivalence :

$$\lambda \in \Lambda _{\epsilon} (A) \Leftrightarrow \sigma_{\min}(A - \lambda I) \le \epsilon \Vert A \Vert$$ où $\sigma_{\min}$ représente la plus petite des valeurs singulières. Les pseudo-spectres sont représentés dans le portrait spectral de la matrice par des lignes de niveaux correspondant à différentes valeurs de $\epsilon$.La construction de ces objets se fait habituellement par le calcul de $\sigma_{\min}(A - \lambda I)$ pour $\lambda$ parcourant une grille définie dans le plan complexe. Comme cette approche est coûteuse en volume de calcul et non totalement fiable, on recherche actuellement des méthodes de continuation qui permettent de suivre directement les lignes de niveaux. Dans tous les cas, ces méthodes reposent sur une utilisation massive du calcul de la plus petite valeur singulière d'une matrice complexe, calcul pour lequel le projet a développé des algorithmes parallèles.

La dichotomie spectrale peut se réaliser grâce à la détermination d'un pseudo-spectre mais il existe aussi une direction de recherche qui consiste à définir des projecteurs spectraux correspondant à des sous-espaces invariants bien séparés, un indicateur évaluant la qualité de la séparation. Le cas de la dichotomie par un cercle ou un axe a été traité par Malyshev [Mal93]. Les cas de dichotomie par une ellipse ou un polygone s'y ramènent. Les procédures obtenues sont très fiables mais si coûteuses qu'elles ne peuvent pour le moment être envisagées que pour des matrices de très petite taille.

Encadrement garanti de résultats

  Mots-clés : conditionnement, intervalles, algèbre linéaire

conditionnement: quantité qui mesure le taux de variation d'un résultat en fonction de l'amplitude des perturbations des données. Cette quantité dépend généralement des normes choisies. Lorsque la résolution du problème correspond au calcul d'une fonction différentiable, le conditionnement est égal à la norme de la matrice jacobienne aux données.

intervalles: l'arithmétique d'intervalles est une extension de l'arithmétique classique qui est mise en oeuvre dans les programmes par le recours à une bibliothèque. Malheureusement, elle n'est généralement pas suffisante car la décorrélation des calculs qu'elle entraîne amène à des surestimations inexploitables. Il est donc nécessaire de définir des procédures spécialisées d'encadrement pour chaque type de problème. [9]

Résumé : Il existe des situations où l'on recherche un encadrement sûr du résultat d'un calcul pour garantir qu'on ne se trouvera pas dans une situation d'échec, malgré l'imprécision des calculs effectués sur ordinateur et l'imprécision des données. C'est le cas par exemple en géométrie algorithmique où des arbres de décision reposent sur des tests arithmétiques. Le projet recherche des procédures d'encadrement pour les problèmes classiques de l'algèbre linéaire.

Lorsque les données d'un problème sont imprécises, elles sont définies par des intervalles. La première solution qui consiste en une application directe de l'arithmétique d'intervalles conduit à une très large surestimation du résultat. Pour traiter ces problèmes, deux grandes classes de méthodes ont donc été développées [Her94]:

- les méthodes itératives : elles consistent à se ramener à des problèmes de type point fixe : si une fonction continue $f$ vérifie $f(I) \subset I$ où I est un compact convexe de ${ \mathbb R}^N$ alors l'équation $f(x)=x$ a au moins une solution dans $I$. Ces méthodes nécessitent habituellement un préconditionnement du système, ce qui produit souvent une large surestimation du résultat.

- les méthodes directes : elles utilisent la structure particulière d'un problème pour obtenir une résolution exacte. Elles nécessitent souvent des conditions d'application très strictes, telles que des conditions de signe sur les vecteurs manipulés.

Les méthodes que nous développons se situent dans une catégorie intermédiaire. Elles consistent à obtenir une description de l'ensemble résultat sous la forme d'une condition nécessaire et suffisante qui ne fasse pas intervenir d'intervalles ; par contre la condition fait intervenir des valeurs absolues et des non linéarités dont on se débarrasse en l'affaiblissant. Nous obtenons alors une approximation du résultat (un sur-ensemble, ou un sous-ensemble, selon les cas) sous la forme d'un polyèdre convexe. Finalement, l'analyse des perturbations au premier ordre nous permet d'obtenir de bons points de départ pour les algorithmes classiques de programmation linéaire. Un traitement numérique des résultats est encore nécessaire, pour obtenir des encadrements sûrs malgré l'utilisation de l'arithmétique flottante.

Nous avons développé nos méthodes d'encadrement pour la résolution de systèmes linéaires et pour la résolution de problèmes aux valeurs propres.

Précédent : Présentation générale et objectifs Remonter : Rapport activite 1997 Suivant : Grands domaines d'application