Décomposer et itérer pour résoudre un problème

« Diviser pour mieux régner » dans le cas des équations aux dérivées partielles

Hors piste Le 24 décembre 2018  - Ecrit par  Guillaume Legendre, Julien Salomon Voir les commentaires (1)

Une idée souvent utilisée en pratique pour calculer la solution d’un problème mathématique compliqué est de résoudre une succession de problèmes plus simples. On pourrait qualifier cette approche intuitive de « résolution par tâtonnements », mais nous allons voir que ce type de procédé peut être rendu systématique, étudié rigoureusement et s’avérer très efficace.

Il arrive parfois que l’on rencontre un problème d’apparence si complexe que l’on ne sache même pas quelle suite d’étapes, même compliquées, pourrait permettre de le résoudre. Une approche générale existe cependant pour contourner cette difficulté : essayer de trouver au sein du problème un ou plusieurs sous-problèmes ayant pour caractéristiques d’être faciles et peu coûteux à résoudre et de pouvoir être intégrés dans une boucle. On peut alors se lancer dans la résolution itérative de la suite éventuellement infinie de ces sous-problèmes. Considérant la question de la stratégie de mise en œuvre réglée et revenant au point de vue des mathématiques, il reste alors à étudier la convergence d’un tel procédé. Si cette dernière a lieu, elle doit permettre de construire la solution du problème initial.

Nous présentons dans cet article différentes déclinaisons de ce procédé, en nous concentrant sur la résolution numérique, c’est-à-dire au moyen d’un calculateur, de problèmes et la construction de sous-problèmes par décomposition. Nous expliquons également comment le formalisme mathématique offre un cadre efficace pour attaquer la question de la convergence.

La résolution itérative sur des exemples

Le principe auquel on va s’intéresser concerne les situations où l’on sait décomposer un problème en une suite de sous-problèmes pouvant chacun être facilement résolu. Lorsque chaque sous-résolution améliore (en un sens à définir au cas par cas) l’approximation de la solution, on peut espérer, en les répétant de manière cyclique un nombre suffisant de fois, arriver à la solution du problème initial. Quand c’est le cas, on constate que la décomposition en sous-problèmes a permis de simplifier les choses. Le prix à payer est cependant d’avoir construit une procédure itérative (engendrant éventuellement une infinité d’itérations).

Une devinette

Commençons par introduire les ingrédients d’une telle méthode sur un premier exemple. Pour cela, on s’intéresse à la question suivante : comment compléter la phrase

« Cette phrase est composée de ... consonnes et de ... voyelles. »

pour la rendre logiquement correcte ?

En suivant une approche par tâtonnements, on commence par compter les consonnes et les voyelles de la phrase contenant les points de suspension, que l’on complète alors avec les deux nombres obtenus pour arriver à :

« Cette phrase est composée de vingt-six consonnes et de dix-neuf voyelles. ».

Évidemment, cette nouvelle phrase est fausse, mais on peut néanmoins y compter les consonnes et les voyelles et la modifier pour tenter de corriger l’erreur. En répétant cette opération tant que la phrase reste fausse, on obtient alors successivement :

« Cette phrase est composée de trente-six consonnes et de vingt-quatre voyelles. »

« Cette phrase est composée de trente-neuf consonnes et de vingt-six voyelles. »

« Cette phrase est composée de trente-huit consonnes et de vingt-cinq voyelles. »

« Cette phrase est composée de trente-neuf consonnes et de vingt-cinq voyelles. »

et l’on aboutit à un résultat correct au bout de cinq essais. On observe qu’on a effectivement pu répondre à la devinette en considérant à chaque étape deux sous-problèmes très faciles à résoudre, consistant simplement à compter les consonnes et les voyelles de la phrase courante [1]. On est donc bien en présence d’un exemple de méthode par tâtonnements, que l’on peut qualifier plus précisément de méthode d’approximations successives : elle est en effet basée d’une part sur la résolution de problèmes « simples » et d’autre part sur la notion de point fixe, l’argument d’entrée d’une itération (une phrase) étant du même type que l’argument de sortie.

Le principe des méthodes de point fixe

Étant donné une application $g$ d’un ensemble $E$ dans lui-même, un point fixe de $g$ est un élément $\xi$ de $E$ tel que $g(\xi)=\xi$, c’est-à-dire un point laissé invariant par l’application $g$. Les points fixes apparaissent dans une grande variété de problèmes mathématiques (comme le montre l’exemple suivant), et notamment la recherche de solutions d’une équation non linéaire.

Une méthode de base pour le calcul d’un point fixe est de recourir à un procédé itératif du type de celui utilisé plus haut. Connaissant une approximation $x^{(0)}$ du point fixe $\xi$, on construit une suite récurrente $(x^{(k)})_{k\in\mathbb{N}}$ définie par la relation de récurrence
\[ x^{(k+1)}=g(x^{(k)}),\ k=0,1,2,\dots \]
Lorsque l’application $g$ est continue, on peut montrer que si la suite $(x^{(k)})_{k\in\mathbb{N}}$ converge vers un élément $\ell$ de $E$ lorsque le nombre d’itérations $k$ tend vers l’infini, cet élément est nécessairement un point fixe de $g$. Cependant, une telle suite peut ne pas converger, même si $g$ possède un point fixe... Il est alors important de disposer de théorèmes garantissant cette convergence moyennant certaines hypothèses supplémentaires sur la fonction $g$.

Observons que l’on aurait pu choisir de procéder autrement à chaque étape en comptant tout d’abord le nombre de consonnes dans la phrase courante, en insérant ce compte dans la phrase, puis, dans un second temps, en comptant le nombre de voyelles dans la phrase résultante. Dans ce cas, le comptage du nombre de voyelles intervient nécessairement après celui du nombre de consonnes et la suite des phrases construites après chaque cycle de deux comptages est la suivante :

« Cette phrase est composée de vingt-six consonnes et de vingt-et-une voyelles. »

« Cette phrase est composée de trente-huit consonnes et de vingt-sept voyelles. »

« Cette phrase est composée de trente-neuf consonnes et de vingt-cinq voyelles. »

Trois essais ont cette fois suffi pour arriver à une phrase logiquement correcte.

Nous avons donc décrit deux versions d’une technique de résolution de la devinette. La première est souvent qualifiée de parallèle, les comptages des consonnes et des voyelles pouvant être effectués simultanément, c’est-à-dire en parallèle. La seconde version est dite séquentielle, dans le sens où le comptage d’un type de lettre ne peut être entrepris qu’une fois réalisé le comptage de l’autre type de lettre, c’est-à-dire à la suite l’un de l’autre. C’est ce vocabulaire que nous emploierons dans le reste de l’article.

Résolution d’un système linéaire

Passons à présent à un problème mathématiquement plus courant : celui de la résolution d’un système d’équations linéaires, comme le système de trois équations linéaires à trois inconnues suivant :

\[ \left\{\begin{aligned} 4\,x_1+ 2\,x_2 -x_3&=5\\ x_1+ 2\,x_2+x_3&=4\\ x_1+x_2+4\,x_3&=6 \end{aligned}\right. \]

dont une solution est $x_1=x_2=x_3=1$, mais cela n’est pas ce qui nous intéresse ici. Une méthode classique pour déterminer la solution de ce problème de manière « directe », c’est-à-dire de manière systématique et avec un nombre fini d’opérations arithmétiques (addition, soustraction, multiplication ou division), est donnée par le fameux procédé d’élimination de Gauss (1809), mais il est aussi possible de considérer une technique qualifiée d’« itérative » : la méthode de Jacobi [Jac1845]. Pour expliquer son fonctionnement, on commence par réécrire le système ci-dessus en isolant une inconnue par équation de la manière suivante :

\[ \left\{\begin{aligned} x_1&=\frac{1}{4}\left(5 - 2x_2+x_3\right)\\ x_2&=\frac{1}{2}\left(4 - x_1- x_3\right)\\ x_3&=\frac{1}{4}\left(6 - x_1- x_2\right) \end{aligned}\right. \]

et l’on introduit ensuite trois suites réelles récurrentes $(x_1^{(k)})_{k\in\mathbb{N}}, (x_2^{(k)})_{k\in\mathbb{N}}$ et $(x_3^{(k)})_{k\in\mathbb{N}}$ définies, à partir de valeurs $x_1^{(0)}$, $x_2^{(0)}$ et $x_3^{(0)}$ données, par les formules de récurrence suivantes :

\[ \begin{aligned} x_1^{(k+1)}&=\frac{1}{4}\left(5-2\,x_2^{(k)} +x_3^{(k)}\right)\\ x_2^{(k+1)}&=\frac{1}{2}\left(4-x_1^{(k)} -x_3^{(k)}\right)\\ x_3^{(k+1)}&=\frac{1}{4}\left(6-x_1^{(k)}-x_2^{(k)}\right) \end{aligned}\ k=0,1,2,\dots \]

L’indice entier $k$ correspond au nombre d’itérations effectuées dans le processus itératif et, bien qu’il soit mathématiquement destiné à tendre vers l’infini, il restera, comme on va le voir, fini en pratique.

Le calcul des termes des suites est particulièrement aisé en raison de la décomposition faite sur le système et l’on peut tester le procédé avec une simple calculatrice en choisissant (par exemple) comme valeurs de départ $x_1^{(0)}=x_2^{(0)}=x_3^{(0)}=0$. Le tableau ci-dessous contient les valeurs des approximations numériques successivement obtenues (on n’indique que cinq chiffres après la virgule) :

Itération $k$ $x_1^{(k)}$ $x_2^{(k)}$ $x_3^{(k)}$
1 1,25000 2,00000 1,50000
2 0,62500 0,62500 0,68750
3 1,10940 1,34380 1,18750
4 0,87500 0,85156 0,88672
5 1,04590 1,11910 1,06840
6 0,95752 0,94287 0,95874
7 1,01820 1,04190 1,02490
8 0,98529 0,97842 0,98497
9 1,00700 1,01490 1,00910
10 0,99483 0,99195 0,99452
11 1,00270 1,00530 1,00330
12 0,99817 0,99702 0,99801
13 1,00100 1,00190 1,00120
14 0,99934 0,99890 0,99927
15 1,00040 1,00070 1,00040
16 0,99976 0,99960 0,99974
17 1,00010 1,00030 1,00020
18 0,99991 0,99985 0,99990

On observe que les valeurs des suites tendent vers la valeur $1$ à la précision affichée, qui correspond effectivement à la valeur de $x_1$, $x_2$ et $x_3$ satisfaisant le système d’équations introduit plus haut. On dit que la méthode itérative converge.

On peut encore remarquer que, à chaque étape, les calculs respectifs des termes des suites à une étape donnée ne font intervenir que des valeurs obtenues à l’étape précédente : on peut donc les calculer de manière simultanée. On est en présence d’une méthode parallèle.

Comme on l’a fait pour le problème de la devinette, on peut trouver une variante séquentielle à la méthode de Jacobi. Il s’agit de la méthode de Gauss-Seidel [Gau1823,Sei1874], dont les formules de récurrence sont :

\[ \begin{aligned} x_1^{(k+1)}&=\frac{1}{4}\left(5-2x_2^{(k)}+x_3^{(k)}\right)\\ x_2^{(k+1)}&=\frac{1}{2}\left(4 -x_1^{(k+1)} -x_3^{(k)}\right)\\ x_3^{(k+1)}&=\frac{1}{4}\left(6 -x_1^{(k+1)} -x_2^{(k+1)}\right) \end{aligned}\ k=0,1,2,\dots \]

En inspectant ces relations, on peut remarquer que le calcul de $x_2^{(k+1)}$ ne peut maintenant se faire qu’une fois la valeur de $x_1^{(k+1)}$ obtenue. De la même manière, le calcul de $x_3^{(k+1)}$ n’est possible que si les valeurs de $x_1^{(k+1)}$ et de $x_2^{(k+1)}$ sont connues. La méthode est donc bien séquentielle. Pour l’exemple considéré, elle conduit aussi à la solution du système. En utilisant la même initialisation que précédemment, on obtient en effet les approximations numériques suivantes (là encore, on n’indique que cinq chiffres après la virgule) :

Itération $k$ $x_1^{(k)}$ $x_2^{(k)}$ $x_3^{(k)}$
1 1,25000 1,37500 0,84375
2 0,77344 1,19141 1,00879
3 0,90649 1,04236 1,01279
4 0,98202 1,00260 1,00385
5 0,99966 0,99825 1,00052
6 1,00101 0,99923 0,99994
7 1,00037 0,99985 0,99995
8 1,00006 1,00000 0,99999
9 1,00000 1,00001 1,00000
10 1,00000 1,00000 1,00000

La convergence observée se trouve même être plus rapide que celle de la méthode de Jacobi [2].

L’usage de l’une ou l’autre ces deux méthodes se généralise facilement à tout système linéaire possédant autant d’équations que d’inconnues et une unique solution. Il est cependant important de remarquer que la convergence de tels procédés n’est a priori pas acquise. Pour le voir, ordonnons différemment les équations du système linéaire considéré plus haut :

\[ \left\{\begin{aligned} x_1+ 2\,x_2+x_3&= 4\\ x_1+x_2+4\,x_3&= 6\\ 4\,x_1+ 2\,x_2 -x_3&=5 \end{aligned}\right. \]

Cette opération ne change bien évidemment pas la solution du système. Appliquons à présent la méthode de Jacobi à la résolution du système réordonné à partir de l’initialisation $x_1^{(0)}=x_2^{(0)}=x_3^{(0)}=0$. Les approximations numériques successivement obtenues sont alors

Itération $k$ $x_1^{(k)}$ $x_2^{(k)}$ $x_3^{(k)}$
1 4 6 -5
2 -3 22 23
3 -63 -83 27
4 143 -39 -423
5 505 1555 489
6 -3595 -2455 5125
7 -211 -16899 -19295
8 53097 77397 -34647
9 -120143 85497 367177
10 -538167 -1348559 -309583

On observe que le comportement de la suite est très différent de ce qu’il était dans le cas précédent. On n’approche effectivement plus la solution, les valeurs absolues des quantités calculées tendant même vers l’infini avec le nombre d’itérations effectuées. Le même phénomène se produit avec la méthode de Gauss-Seidel.

Dans ces conditions, comment savoir que la méthode fonctionnera ou pas ? La réponse est donnée par le domaine des mathématiques portant le nom d’analyse numérique, qui est dédié à l’étude des diverses méthodes numériques proposées depuis des siècles pour résoudre des problèmes de natures variées. Il existe ainsi des théorèmes fixant des conditions sur le système linéaire garantissant la convergence vers la solution des suites d’approximations construites par les méthodes de Jacobi ou de Gauss-Seidel.

Un point sur les méthodes itératives linéaires de résolution de système linéaire

Les méthodes de Jacobi et de Gauss-Seidel décrites ci-dessus sont deux exemples emblématiques d’une classe de méthodes, dites linéaires, de résolution de système linéaire, basée sur un principe de point fixe. Expliquons brièvement le fonctionnement de ces méthodes d’un point de vue matriciel.

Supposons que l’on souhaite résoudre un système de $n$ équations linéaires à $n$ inconnues de la forme $AX=B$, où $A$ est une matrice inversible, $B$ une matrice colonne donnée et $X$ une matrice colonne ayant pour coefficients les inconnues du système.

Pour définir une méthode linéaire, il suffit d’écrire la matrice du système comme la somme de matrices $A=M-N$, avec pour contraintes que la matrice $M$ de cette décomposition soit inversible et qu’un système linéaire ayant $M$ pour matrice soit, d’un point de vue calculatoire, « facile » à résoudre. Le procédé itératif définissant un telle méthode étant en effet basé sur la formule de récurrence

\[ MX^{(k+1)}=NX^{(k)}+B,\ k=0,1,2,\dots \]

qui correspond à la résolution d’un tel système linéaire, il est de toute première importance que cette dernière soit bien plus aisée que celle du système initial $AX=B$.

Dans ce cadre, la méthode de Jacobi fait le choix pour $M$ de la matrice diagonale extraite de $A$ (un système linéaire diagonal se résolvant trivialement), la méthode de Gauss-Seidel celui de la partie triangulaire inférieure de $A$ (un système triangulaire étant systématiquement et facilement résoluble par substitutions successives) [3]. D’autres choix pertinents existent, donnant lieu à d’autres méthodes.

La convergence d’une méthode itérative linéaire est assurée si les valeurs propres de la matrice $M^{-1}N$, appelée matrice d’itération de la méthode, ont toutes un module strictement plus petit que $1$. Suivant les caractéristiques (symétrique définie positive, dominance des coefficients diagonaux...) de la matrice $A$ du système à résoudre, il existe des résultats garantissant que c’est le cas. La convergence d’une méthode donnée est alors d’autant plus rapide que le plus grand module de ses valeurs propres (le rayon spectral de la matrice) est petit. Dans l’exemple traité plus haut, la matrice d’itération de la méthode de Jacobi possède un rayon spectral qui vaut environ $0,60355$, alors que celui de la matrice d’itération de la méthode de Gauss-Seidel est égal à $0,25$ expliquant la convergence plus rapide de cette dernière méthode. En revanche, pour le cas de divergence également présenté, le rayon spectral de la matrice d’itération de la méthode de Jacobi vaut environ $3,8836$ et est donc plus grand que $1$.

En pratique, les méthodes itératives constituent une alternative aux méthodes directes comme l’élimination de Gauss, ces dernières pouvant s’avérer trop coûteuses en termes d’espace mémoire pour être mises en œuvre lorsque la taille du système à résoudre est grande. Elles s’avèrent aussi plus résilientes à l’accumulation des erreurs d’arrondi (inévitable lorsque que les calculs sont effectués dans une arithmétique en précision finie) que certaines méthodes directes.

Le lecteur pourra être interpellé par le fait qu’un nombre infini d’opérations arithmétiques est a priori nécessaire pour calculer la solution du système. En pratique, ceci n’est jamais un problème : les calculs étant effectués dans une arithmétique en précision finie, une solution approchée suffisamment précise sera indistinguable de la solution exacte telle que représentée sur la machine et l’on peut donc mettre fin au procédé après un nombre fini d’itérations, par exemple dès que l’approximation de la solution n’est plus modifiée significativement ou plus simplement dès qu’une précision attendue est atteinte.

« Diviser pour régner »

Nous venons de voir comment un procédé itératif et une décomposition adaptée d’un problème peuvent donner lieu à des méthodes de résolution efficaces. Cette manière d’attaquer la résolution d’un problème complexe est aujourd’hui classique et le lecteur averti aura peut-être reconnu le principe fondateur du paradigme « diviser pour régner », qui est à la base de nombreuses méthodes de traitement de données en informatique.

Décomposition de domaine et décomposition de frontière

Nous allons maintenant montrer comment ces approches itératives peuvent être mises à profit pour résoudre certaines équations aux dérivées partielles. Pour cela, on exploite de manière cruciale la linéarité des équations du problème considéré : au travers de ce que l’on appelle un principe de superposition, cette propriété rend en effet possible l’écriture de la solution du problème sous la forme d’une somme de solutions de sous-problèmes qui lui sont liés.

La méthode de Schwarz

JPEG - 48.4 ko
Hermann Schwarz

Considérons tout d’abord une méthode inventée en 1870 par Hermann Schwarz [4] pour démontrer l’existence d’une solution au problème de Dirichlet posé dans un domaine consistant en l’union d’un carré et d’un disque se chevauchant partiellement [Sch1870]. Dans ce problème, l’inconnue est une fonction satisfaisant une équation aux dérivées partielles appelée équation de Laplace [5].

Schwarz eut l’idée d’utiliser le fait que des solutions du problème de Dirichlet dans des domaines de géométrie simple (comme le carré et le disque) étaient explicitement connues de longue date [6]. Il construisit ainsi une solution du problème sur tout domaine décomposable en sous-domaines « simples » en alternant des résolutions sur chacun de ces derniers et en combinant les solutions obtenues.

Le principe de Dirichlet et la méthode de Schwarz

Les travaux de Schwarz portaient sur la recherche d’un outil permettant de démontrer le principe de Dirichlet dans des domaines de forme compliquée. Ce principe permet d’établir un lien entre la solution du problème de Dirichlet et un problème de minimisation. La solution du problème de Dirichlet est la fonction $u$ satisfaisant d’une part l’équation de Laplace
\[ \Delta u=0 \]
dans un domaine ouvert borné $\Omega$ de $\mathbb{R}^d$ (avec $d=2$ dans l’exemple de Schwarz), et d’autre part la contrainte de bord
\[ u=g \]
sur la frontière $\partial\Omega$ de $\Omega$, avec $g$ une fonction donnée.
On peut montrer que $u$ minimise la quantité
\[ \int_\Omega|\nabla v(x)|^2\,\mathrm{d}x, \]
appelée énergie de Dirichlet, parmi toutes les fonctions $v$ deux fois différentiables vérifiant la condition aux limites ci-dessus.

PNG - 6.4 ko
Le domaine $\Omega$ considéré par Schwarz.

Schwarz utilisa un domaine $\Omega$ bidimensionnel auquel il associa une décomposition en deux sous-domaines de géométrie simple et présentant un recouvrement : un disque $\Omega_1$ et un carré $\Omega_2$ (voir figure ...). Pour prouver l’existence d’une solution au problème de Dirichlet dans $\Omega$, il introduisit une méthode itérative utilisant uniquement des solutions du problème de Dirichlet dans $\Omega_1$ et $\Omega_2$, qui peuvent être obtenues au moyen de séries de Fourier.

L’algorithme définissant cette méthode peut être résumé de la façon suivante. Partant d’une estimation, notée $u_2^{(0)}$, de la solution au niveau de l’interface $\Gamma_1=\partial\Omega_1\cap\Omega_2$, on calcule de manière successive les fonctions $u_1^{(k+1)}$ et $u_2^{(k+1)}$, $k=0,1,2,\dots$, solutions respectives de
\[ \left\{ \begin{aligned} \Delta u_1^{(k+1)}&=0\text{ dans }\Omega_1,\\ u_1^{(k+1)}&=g\text{ sur }\partial\Omega_1\cap\partial\Omega,\\ u_1^{(k+1)}&=u_2^{(k)}\text{ sur }\Gamma_1=\partial\Omega_1\cap\Omega_2, \end{aligned} \right. \quad \left\{ \begin{aligned} \Delta u_2^{(k+1)}&=0\text{ dans }\Omega_2,\\ u_2^{(k+1)}&=g\text{ sur }\partial\Omega_2\cap\partial\Omega,\\ u_2^{(k+1)}&=u_1^{(k+1)}\text{ sur }\Gamma_2=\partial\Omega_2\cap\Omega_1. \end{aligned} \right. \]
On peut voir que le recouvrement des deux sous-domaines garantit que la seconde condition aux limites de chacun des problèmes a bien un sens. On observe également que la résolution des problèmes se fait nécessairement par alternance (Schwarz utilise d’ailleurs dans sa présentation l’analogie d’un système de pompe à vide manuelle), donc de manière séquentielle.

Au moyen du principe du maximum, Schwarz montra que les valeurs des différences $u_i^{(k+1)}-u_i^{(k)}$ sur l’interface $\Gamma_i$, $i=1,2$, décroissent en valeur absolue géométriquement au fur et à mesure que l’entier $k$ augmente, ce qui lui permit de déduire que les suites $(u_1^{(k)})_{k\in\mathbb{N}^*}$ et $(u_2^{(k)})_{k\in\mathbb{N}^*}$ sont convergentes, de limites respectives
\[ u_1=u_1^{(1)}+(u_1^{(2)}-u_1^{(1)})+(u_1^{(3)}-u_1^{(2)})+\dots\text{ et }u_2=u_2^{(1)}+(u_2^{(2)}-u_2^{(1)})+(u_2^{(3)}-u_2^{(2)})+\dots \]
par simple sommation de série. Observant enfin que les fonctions $u_1$ et $u_2$ coïncident sur $\Gamma_1$ et $\Gamma_2$, il conclut qu’elles sont identiques dans le recouvrement $\Omega_1\cap\Omega_2$ et qu’elles correspondent par conséquent à des restrictions de la solution du problème de Dirichlet dans le domaine $\Omega$.

PNG - 174.6 ko
Illustration de la convergence de la méthode de Schwarz pour le problème de Dirichlet : représentation des lignes de niveau de l’approximation numérique (calculée avec le logiciel Freefem++) de la solution après 1 (en haut à gauche), 2 (en haut à droite), 3 (en bas à gauche) et 4 (en bas à droite) itérations du procédé.

Si la motivation initiale de cette technique de décomposition de domaine était avant tout théorique, sa généralisation à la résolution d’autres problèmes que ceux liés à l’équation de Laplace et son étude mathématique réalisée avec des outils modernes d’analyse fonctionnelle lui firent par la suite connaître un renouveau en tant que méthode numérique à partir des années 1960. Elle est à ce jour, suite à de nombreuses extensions, un sujet de recherche actif en analyse numérique et en calcul scientifique.

La méthode des réflexions

GIF - 80.1 ko
Séquence d’images de sphères de polymères qui sédimentent dans un tube étroit (10 fois le diamètre des sphères).

La méthode des réflexions est une consœur de la méthode de Schwarz, qui reste cependant moins connue de la communauté mathématique que cette dernière. En revanche, elle est depuis plusieurs décennies présente dans les ouvrages traitant de micro-hydrodynamique, domaine de la physique traitant des écoulements à petite échelle.

Elle fut introduite en 1911 par le physicien polonais Marian Smoluchowski pour étudier le phénomène de sédimentation en déterminant les vitesses respectives de $n$ particules solides de petite taille en suspension dans un fluide visqueux et soumises à des efforts extérieurs (typiquement l’action de la gravité) [Smo1911]. Le problème mathématique sous-jacent, parfois appelé problème de mobilité [7], fait alors intervenir un système d’équations aux dérivées partielles portant le nom d’équation de Stokes stationnaire.

Le problème de mobilité en micro-hydrodynamique

La micro-dynamique s’intéresse à l’évolution de solides placés dans un fluide, d’échelles petites devant les grandeurs caractérisant l’écoulement. Lorsque le fluide est visqueux et que le mouvement des solides, ou particules, est très lent, les effets inertiels et instationnaires sont négligeables et l’on se trouve en présence de ce qu’on appelle un écoulement à faible nombre de Reynolds. Dans ce cas, on a coutume de considérer que l’écoulement est régi par l’équation de Stokes stationnaire,
\[ -\mu\,\Delta\vec{u}+\vec{\nabla}p=\vec{0}\text{ dans }\mathbb{R}^3\setminus\bigcup_{i=1}^n\overline{O_i}, \]
dans laquelle $\vec{u}$ et $p$ sont respectivement le champ de vitesse (un vecteur à trois composantes) et le champ de pression (un scalaire) de l’écoulement, $\mu$ est la viscosité dynamique du fluide, l’opérateur $\Delta$ est le laplacien, l’opérateur $\vec{\nabla}$ est le gradient et le domaine ouvert $O_i$ désigne l’espace occupé dans le fluide (supposé emplir tout l’espace) par la $i$e des $n$ particules.

Pour définir le problème complètement, il faut ajouter encore ajouter quelques équations :

  • une condition d’incompressibilité
    \[ \text{div}\,\vec{u}=0\text{ dans }\mathbb{R}^3\setminus\bigcup_{i=1}^n\overline{O_i}, \]
    dans laquelle l’opérateur $\text{div}$ est la divergence,
  • une condition aux limites d’adhérence du fluide, traduisant le fait que la vitesse du fluide au bord d’une particule est égale à la vitesse de cette particule,
    \[ \begin{equation}\label{vitesse rigide} \forall\vec{x}\in\partial O_i, \vec{u}(\vec{x})=\vec{V_i}+\vec{\Omega_i}\wedge(\vec{x}-\vec{x_i}),\ i=1,\dots,n, \end{equation} \]
    où les vecteurs $\vec{V_i}$, $\vec{\Omega_i}$ et $\vec{x_i}$ sont respectivement la vitesse de translation, la vitesse de rotation et le centre de masse de la $i$e particule et où $\partial O_i$ désigne le bord de la $i$e particule,
  • une équation rendant compte des forces extérieures, de la forme
    \[ \int_{\partial O_i}\frac{\partial\vec{u}}{\partial\vec{n}}\,\mathrm{d}\sigma=\vec{F_i},\ i=1,\dots,n. \]
    où le vecteur $\vec{F_i}$ est la résultante des forces appliquées à la $i$e particule et le vecteur $\vec{n}$ est le vecteur unitaire normal au bord $\partial O_i$, orienté vers l’intérieur de la particule,
  • une condition décrivant le comportement asymptotique du champ de vitesse à l’infini, prenant la forme
    \[ \lim_{\|\vec{x}\|\to+\infty}\vec{u}(\vec{x})=\vec{0} \]
    si l’on suppose le fluide initialement au repos.

Ce système d’équations possède une unique solution, et résoudre mathématiquement ce que l’on nomme le problème de mobilité revient à trouver le champ $\vec{u}$ et les vecteurs $\vec{V_i}$ et $\vec{\Omega_i}$, $i=1,\dots,n$, les vecteurs $\vec{F_i}$ et les positions $\vec{x_i}$ étant connus.

Pour une particule isolée sphérique de rayon $a$ (correspondant au cas $n=1$ et au choix d’une sphère de rayon $a$ pour $O_1$ dans le problème décrit ci-dessus), la loi de friction de Stokes fournit une solution explicite. En effet, on a
\[ \vec{V_1}=-\frac{\vec{F_1}}{6\pi\mu a}, \]
ce qui signifie que la particule acquiert un mouvement de translation, opposé et proportionnel à la résultante des forces appliquées.

JPEG - 52.8 ko
Marian Smoluchowski

Dans ses investigations, Smoluchowski souhaitait arriver par le biais de la modélisation mathématique à quantifier l’effet sur l’écoulement des interactions hydrodynamiques entre les particules solides, phénomène observé expérimentalement. Or, si la solution du problème à une seule particule était connue depuis les travaux de Stokes sous la forme d’une expression plus ou moins simple mais néanmoins explicite [8], il ne disposait d’aucun outil de calcul permettant de traiter le cas d’un nombre arbitraire de particules. Inspiré par la méthode des images en électromagnétisme, il proposa un procédé itératif d’approximation de la solution utilisant uniquement la résolution de problèmes ne faisant chacun intervenir qu’une seule particule.

Pour cela, son idée était de partir d’un écoulement ne contenant aucune particule et de le corriger progressivement de manière à ce que la vitesse du fluide soit tour à tour égale à la vitesse du bord de chacune des particules. Dans le cadre de la mécanique des fluides, cette correction s’apparente à la réaction (ou réflexion) d’une particule placée dans un écoulement donné, ce qui explique le nom donné à la méthode.

JPEG - 62.7 ko
Gennadi Michailowitsch Golusin

Chaque nouvelle correction sur une particule vient perturber la vitesse sur toutes les autres particules et l’on est donc amené à parcourir l’ensemble du système de particules considéré de manière cyclique et récurrente pour tenter de reconstruire les interactions hydrodynamiques absentes des problèmes à une particule qui sont résolus.

Ce faisant, la méthode des réflexions décrite par Smoluchowski suit un principe itératif séquentiel permettant de résoudre un problème posé dans un domaine dont la frontière possède plusieurs composantes d’un seul tenant bornées (en d’autres termes, le domaine présente des « trous » et est dit non simplement connexe), modélisant les bords de particules, d’obstacles ou plus généralement d’« objets » présents. Une variante parallèle de la méthode fut introduite indépendamment par Golusin en 1934 [Gol1934] et appliquée à la résolution d’un problème de Dirichlet dans un tel domaine.

La méthode des réflexions appliquée à la résolution d’un problème de Dirichlet

Tentons de faire un parallèle entre la méthode de Schwarz décrite plus haut et la méthode des réflexions. Considérons pour cela un problème classique d’électrostatique consistant à déterminer le potentiel électrique dans un domaine ouvert $\Omega$ contenant deux objets conducteurs, respectivement représentés par les sous-domaines ouverts $O_1$ et $O_2$, au bord desquels le potentiel est connu et respectivement donné par des fonctions $g_1$ et $g_2$. Dans ce cas, le potentiel est solution d’un problème de Dirichlet prenant alors la forme du système d’équations suivant :
\[ \left\{ \begin{aligned} \Delta u&=0\text{ dans }\Omega\setminus(\overline{O_1}\cup\overline{O_2}),\\ u&=g_1\text{ sur }\partial O_1,\\ u&=g_2\text{ sur }\partial O_2,\\ \end{aligned} \right. \]
que l’on complète par une condition de Dirichlet homogène sur le bord du domaine $\Omega$ (que l’on suppose borné pour éviter les complications techniques),
\[ u=0\text{ sur }\partial\Omega. \]

Pour la résolution de ce problème, la méthode des réflexions introduite par Smoluchowski construit deux suites récurrentes $(u_1^{(k)})_{k\in\mathbb{N}^*}$ et $(u_2^{(k)})_{k\in\mathbb{N}^*}$, les fameuses « réflexions », dont les termes sont les solutions respectives de
\[ \left\{ \begin{aligned} \Delta u_1^{(1)}&=0\text{ dans }\Omega\setminus\overline{O_1},\\ u_1^{(1)}&=g_1\text{ sur }\partial O_1,\\ u_1^{(1)}&=0\text{ sur }\partial\Omega, \end{aligned} \right. \quad \left\{ \begin{aligned} \Delta u_2^{(1)}&=0\text{ dans }\Omega\setminus\overline{O_2},\\ u_2^{(1)}&=g_2-u_1^{(1)}\text{ sur }\partial O_2,\\ u_2^{(1)}&=0\text{ sur }\partial\Omega. \end{aligned} \right. \]
et, pour $k\geq1$,
\[ \left\{ \begin{aligned} \Delta u_1^{(k+1)}&=0\text{ dans }\Omega\setminus\overline{O_1},\\ u_1^{(k+1)}&=-u_2^{(k)}\text{ sur }\partial O_1,\\ u_1^{(k+1)}&=0\text{ sur }\partial\Omega, \end{aligned} \right. \quad \left\{ \begin{aligned} \Delta u_2^{(k+1)}&=0\text{ dans }\Omega\setminus\overline{O_2},\\ u_2^{(k+1)}&=-u_1^{(k+1)}\text{ sur }\partial O_2,\\ u_2^{(k+1)}&=0\text{ sur }\partial\Omega. \end{aligned} \right. \]
On voit que le problème pour la fonction $u_2^{(k+1)}$ fait intervenir la fonction $u_1^{(k+1)}$, la méthode est donc séquentielle par nature. L’approximation fournie par la méthode après un nombre $\ell$ de cycles est obtenue en sommant les restrictions à $\Omega\setminus(\overline{O_1}\cup\overline{O_2})$ des termes des deux suites :
\[ \forall\ell\in\mathbb{N}^*,\ u^{(\ell)}(x)=\sum_{k=1}^\ell\left(u_1^{(k)}(x)+u_2^{(k)}(x)\right)\text{ pour }x\in\Omega\setminus(\overline{O_1}\cup\overline{O_2}). \]
On peut facilement montrer que si la série correspondante converge, c’est bien vers la solution $u$ du problème de Dirichlet. Tout l’enjeu de l’étude de la méthode est donc de prouver cette convergence. En attendant un tel résultat, on peut revenir sur un propriété des trois suites introduites ci-dessus. On peut en effet montrer, par un simple raisonnement par récurrence, qu’elles sont telles que
\[ \forall \ell\in\mathbb{N}^*,\ u^{(\ell)}+\sum_{j=1}^iu_j^{(\ell+1)}=g_i\text{ sur }\partial O_i,\ i=1,2. \]
Une réflexion est donc à un champ qui corrige les conditions au bord d’un objet - et d’un objet seulement - en utilisant le principe des images. Cette correction est nécessaire car les réflexions sur les autres objets ont introduit des erreurs au cours du cycle.

PNG - 109.5 ko
Illustration du principe de la méthode des réflexions pour un problème de Dirichlet à deux objets : la solution du problème (à gauche) est obtenue de manière constructive en superposant deux contributions, l’une liée à l’objet $O_1$ (au centre), l’autre liée à l’objet $O_2$ (à droite), respectivement données par les sommes $\sum_{k=1}^{+\infty}u_1^{(k)}$ et $\sum_{k=1}^{+\infty}u_2^{(k)}$ des réflexions correspondantes (approximations numériques calculées avec le logiciel Freefem++).

Dans la variante de Golusin, on procède de façon similaire et les termes respectifs des deux suites sont solutions de
\[ \left\{ \begin{aligned} \Delta u_1^{(1)}&=0\text{ dans }\Omega\setminus\overline{O_1},\\ u_1^{(1)}&=g_1\text{ sur }\partial O_1,\\ u_1^{(1)}&=0\text{ sur }\partial\Omega, \end{aligned} \right. \quad \left\{ \begin{aligned} \Delta u_2^{(1)}&=0\text{ dans }\Omega\setminus\overline{O_2},\\ u_2^{(1)}&=g_2\text{ sur }\partial O_2,\\ u_2^{(1)}&=0\text{ sur }\partial\Omega. \end{aligned} \right. \]
et, pour $k\geq1$,
\[ \left\{ \begin{aligned} \Delta u_1^{(k+1)}&=0\text{ dans }\Omega\setminus\overline{O_1},\\ u_1^{(k+1)}&=-u_2^{(k)}\text{ sur }\partial O_1,\\ u_1^{(k+1)}&=0\text{ sur }\partial\Omega, \end{aligned} \right. \quad \left\{ \begin{aligned} \Delta u_2^{(k+1)}&=0\text{ dans }\Omega\setminus\overline{O_2},\\ u_2^{(k+1)}&=-u_1^{(k)}\text{ sur }\partial O_2,\\ u_2^{(k+1)}&=0\text{ sur }\partial\Omega. \end{aligned} \right. \]
Pour un problème à deux objets comme celui étudié, seul le second problème est modifié (ou les $n-1$e derniers problèmes d’un cycle si $n$ objets étaient considérés) et fait maintenant intervenir une fonction calculée lors du cycle précédent : la méthode est cette fois parallèle. Pour cette méthode, les réflexions sont telles que l’on a
\[ \forall \ell\in\mathbb{N}^*,\ u^{(\ell)}+u_i^{(\ell+1)}=g_i\text{ sur }\partial O_i,\ i=1,2. \]

Une caractéristique de la méthode des réflexions la rend intéressante d’un point de vue pratique. Dans cette méthode, la frontière du domaine sur lequel est posé le problème est en effet décomposée en tenant compte des bords respectifs des objets. On parle de décomposition de frontière en composantes connexes. Nous avons ainsi vu que les sous-problèmes résolus au cours d’un cycle de la méthode ne font intervenir qu’un objet et donc une seule composante du bord à la fois. Ceci est particulièrement adapté à l’utilisation de formulations du problème basées sur la représentation intégrale de sa solution. Cette technique se sert du fait que, pour certains problèmes, la valeur de la solution en tout point intérieur du domaine est donnée par celle d’une intégrale ne faisant intervenir que ses valeurs (et celles de sa dérivée normale) au bord du même domaine. En se servant d’une telle représentation, il est possible de reformuler le système d’équations aux dérivées partielles du problème considéré en un système équivalent d’équations intégrales, chacune étant relative à l’une des composantes de la frontière. Dans ce cas particulier, l’utilisation de la méthode des réflexions consiste à résoudre itérativement le système obtenu en considérant une à une ces équations intégrales, comme on l’avait fait avec les équations algébriques du système linéaire en début d’article, et en procédant comme dans la méthode de Gauss-Seidel (pour la version de Smoluchowski) ou de Jacobi (pour la version de Golusin).

On notera pour finir qu’à l’époque de son introduction, bien antérieure à l’apparition des premiers calculateurs électroniques, la méthode des réflexions permettait, pour une valeur du nombre d’objets $n$ de quelques unités et en n’effectuant que deux ou trois cycles au total, de produire un développement en série tronqué de la solution en accord avec les résultats expérimentaux disponibles en micro-hydrodynamique. Elle a connu ces dernières années diverses applications numériques mais aussi théoriques, puisqu’elle a servi d’outil de démonstration constructive d’existence de solution. Comme c’est souvent le cas pour une « bonne » méthode, elle a même été réinventée au début des années 2000 pour résoudre des problèmes de diffraction multiple [Bal1997,Bal2004].

L’intérêt de la formalisation et de l’analyse mathématique

Les méthodes de Schwarz et des réflexions présentent certaines similarités, comme un caractère itératif ou le recours à une décomposition géométrique du domaine, mais ce ne sont pas les mêmes méthodes. L’un des points forts des mathématiques est de proposer des formalisations qui permettent de s’affranchir des cas particuliers, ouvrant la voie à une analyse générale s’appliquant dans de nombreuses situations. Pour comprendre comment ce principe permet de relier les deux méthodes, concentrons-nous sur une interprétation abstraite de la méthode des réflexions faisant appel à deux outils déjà connus : la méthode des projections orthogonales d’une part et la méthode des projections alternées d’autre part.

La première idée consiste à voir la fonction solution de chaque problème aux limites résolu comme le résultat de la projection d’une fonction liée aux données du problème. On a en effet dit plus haut qu’une réflexion modifiait les valeurs de l’approximation courante de la solution pour que cette dernière vérifie la condition aux limites imposée au bord de l’objet associé. On peut aisément se convaincre qu’une telle opération est idempotente, c’est-à-dire que l’effectuer une ou plusieurs fois conduit au même résultat, et qu’elle est par conséquent une projection, au sens algébrique du mot. En concevant ainsi les choses, on ramène la méthode de Smoluchowski à l’application cyclique de projecteurs, chacun étant associé à un objet, à une approximation initiale de la solution. La méthode entre alors dans le cadre de la méthode des projections alternées et il ne reste qu’à faire appel aux résultats déjà existants pour cette dernière pour prouver sa convergence.

GIF - 281.1 ko
John von Neumann

Détaillons davantage ce dernier point. On sait depuis un travail de John von Neumann dans les années 1930 [Neu1949] que la projection orthogonale sur l’intersection de deux sous-espaces peut être obtenue en projetant orthogonalement sur un sous-espace, puis sur l’autre de façon répétée. On peut illustrer ce fait de la manière suivante. Supposons que l’on veuille déterminer le point d’intersection, noté $I$, de deux droites $(d1)$ et $(d2)$ du plan en ne disposant que d’une approximation $I^{(0)}$ de ce point et en ne sachant que projeter orthogonalement un point quelconque du plan sur l’une ou l’autre de ces droites. Dans ce cas, on projette $I^{(0)}$ sur $(d1)$, puis sur $(d2)$ pour obtenir le point $I^{(1)}$, on réitère le procédé avec le point $I^{(1)}$ et ainsi de suite...

PNG - 6.6 ko
Détermination du point d’intersection de deux droites sécantes par la méthode des projections alternées

On se convainc aisément, en observant la figure ci-contre, que la suite $(I^{(k)})_{k\in\mathbb{N}}$ des points construits de cette manière converge vers le point $I$.

Le procédé de von Neumann a été généralisé à un nombre fini de projecteurs par Israel Halperin dans les années 1960 [Hal1962] et s’applique à d’autres espaces que le plan géométrique et d’autres sous-espaces que des droites géométriques, puisqu’on peut l’utiliser dans des espaces fonctionnels de dimension infinie. Par ce biais, on a trouvé un moyen de ramener la preuve de convergence de la méthode des réflexions à celle de l’orthogonalité d’un projecteur. Ajoutons que si cette approche ne fonctionne pas avec la version de Golusin du procédé, il est néanmoins possible de modifier légèrement cette dernière pour parvenir à la même conclusion.

Cette technique de formalisation s’avère suffisamment générale pour avoir été utilisée par Pierre-Louis Lions vers la fin des années 1980 afin de donner une preuve variationnelle de convergence de... la méthode de Schwarz [Lio1988].

L’orthogonalité des projections pour un problème de Dirichlet

Dans le cas du problème de Dirichlet à deux objets décrit plus haut, le caractère orthogonal des projecteurs intervenant dans la méthode des réflexions est simple à établir. Pour le voir, il faut cependant introduire une formulation variationnelle du problème.

En effet, qui dit orthogonalité suppose que l’espace (vectoriel) dans lequel on travaille est muni d’un produit scalaire. Ici, la structure requise va naturellement apparaître en posant le problème sous une forme variationnelle. Celle-ci est obtenue en multipliant l’équation du problème à résoudre par une fonction, dite fonction test, choisie dans un espace fonctionnel adéquat et en intégrant par parties sur le domaine $\Omega$. Dans le cas considéré, l’espace fonctionnel en question est l’espace de Sobolev $H^1_0(\Omega)$, c’est-à-dire l’ensemble des fonctions définies sur le domaine $\Omega$, qui s’annulent sur sur le bord $\partial\Omega$, sont de carré intégrable et dont les dérivées partielles (au sens des distributions) sont également de carré intégrable.

Rappelons que les réflexions sont définies par des problèmes aux limites, dont chacun est posé à l’extérieur d’un des objets $O_i$, avec $i$ égal à $1$ ou $2$, de la forme :
\[ \left\{ \begin{aligned} \Delta u&=0\text{ dans }\Omega\setminus\overline{O_i},\\ u&=v\text{ sur }\partial O_i,\\ u&=0\text{ sur }\partial\Omega, \end{aligned} \right. \]
où $v$ est une fonction donnée définie sur $\Omega$. On peut étendre de façon naturelle les solutions de ces sous-problèmes à l’intérieur de l’objet concerné en posant simplement :
\[ u=v\text{ dans }O_i. \]
La formulation variationnelle associée aux sous-problèmes est alors : étant donné une fonction $v$ dans $H^1_0(\Omega)$, trouver $u$ dans $H^1_0(\Omega)$ vérifiant $u=v$ dans $\overline{O_i}$ et telle que
\[ \int_\Omega\nabla u(x)\cdot\nabla w(x)\,\mathrm{d}x=0 \]
pour toute fonction $w$ dans $H^1_0(\Omega)$ vérifiant de plus $w=0$ dans $\overline{O_i}$, où le symbole $\nabla$ désigne l’opérateur gradient et $\cdot$ désigne le produit scalaire euclidien entre deux vecteurs. On remarque que cette formulation fait intervenir la forme (bilinéaire et symétrique)
\[ (u,v)\mapsto\int_\Omega\nabla u(x)\cdot\nabla v(x)\,\mathrm{d}x, \]
qui est effectivement un produit scalaire sur $H^1_0(\Omega)$.

On définit les projecteurs associés aux réflexions comme les applications $P_i$, $i=1,2$, de $H^1_0(\Omega)$ dans lui-même, telles que $P_i(v)=u$, où $u$ est la solution du sous-problème, posé sous forme variationnelle, ci-dessus.

Il est tout d’abord clair que ces applications sont idempotentes. Considérons ensuite l’image $V_i$ de l’application $P_i$, qui est le sous-espace de $H^1_0(\Omega)$ constitué des fonctions qui sont harmoniques, c’est à dire dont à laplacien nul, à l’extérieur de l’objet $O_i$. Il est bien connu depuis les travaux de Weyl en théorie du potentiel [Wey1940] que l’espace $H^1_0(\Omega)$ peut s’écrire comme la somme directe orthogonale (relativement au produit scalaire que l’on a introduit) du sous-espace $V_i$ et d’un sous-espace $M_i$, constitué des fonctions de $H^1_0(\Omega)$ s’annulant dans et au bord de l’objet $O_i$. Il en découle que les applications $P_i$, $i=1,2$, sont des projecteurs orthogonaux.

Décomposition de domaine, diversité des techniques et nouveaux enjeux

Les différentes méthodes que nous venons de passer en revue se servent toutes d’une décomposition pour résoudre de manière itérative un problème considéré comme trop compliqué ou trop coûteux pour être abordé directement. Elles permettent en pratique d’obtenir des approximations de la solution du problème, car on ne peut généralement pas réaliser la totalité des itérations conduisant à la solution, en particulier lorsque cette dernière n’est obtenue qu’à l’issue d’un nombre infini de ces itérations. La convergence du procédé itératif utilisé est alors certifiée
par l’analyse mathématique, la profondeur des outils mathématiques requis par celle-ci étant bien sûr liée à la nature du problème à résoudre.

Les méthodes de décomposition constituent en soi un champ de recherche complet, engageant aujourd’hui de nombreux chercheurs au sein d’une communauté très active. Les outils mathématiques qui leur sont liés sont variés : la modélisation qui conduit des phénomènes considérés à des systèmes d’équations, l’analyse fonctionnelle pour l’étude de ces équations, l’analyse numérique pour celles des algorithmes de résolution, le calcul scientifique pour leur mise-en-œuvre pratique... On voit que la décomposition de domaine offre en fait une porte d’entrée de choix vers les mathématiques appliquées !

Certains défis actuels concernant la résolution numérique des équations aux dérivées partielles sont apportés par de nouveaux types de décomposition. La décomposition en temps en est un exemple : paralléliser par rapport au temps la résolution d’une équation d’évolution ajoute comme difficulté le fait que l’écoulement du temps est profondément séquentiel, dans le sens où il est difficile de prévoir l’état d’un système dynamique à un certain instant sans avoir connaissance de son état à des instants antérieurs. Un autre exemple de recherche émergente dans ce domaine est le couplage de méthodes de décomposition avec d’autres procédures itératives, issues de l’optimisation par exemple.

Références

[Bal1997]
M. Balabane et V. Tirel, Boundary decomposition for Helmholtz and Maxwell equations 1 : disjoint sub-scatterers, Asymptotic Anal., 38(1), pp. 281—286, 1997.

[Bal2004]
M. Balabane, Boundary decomposition for Helmholtz and Maxwell equations 1 : disjoint sub-scatterers, Asymptotic Anal., 38(1), pp. 281—286, 2004.

[Gau1823]
C. F. Gauss, Lettre du 26 décembre adressée à Christian Ludwig Gerling, 1823.

[Jac1845]
C. G. J. Jacobi, Ueber eine neue Auflösungsart der bei der Methode der kleinsten Quadrate vorkommenden lineären Gleichungen, Astronom. Nachr., 22(20), pp. 297—306, 1845.

[Gol1934]
G. M. Golusin, Auflösung einiger ebenen Grundaufaben der mathematischen Physik im Fall der Laplaceschen Gleichung und mehrfachzusammenhängender Gebiete, die durch Kreise begrenzt sind, Mat. Sb., 41(2), pp. 246—276, 1934.

[Hal1962]
I. Halperin, The product of projection operators, Acta Sci. Math. (Szeged), 23(1-2), pp. 96—99, 1962.

[Lio1988]
P.-L. Lions, On the Schwarz alternating method. I, First international symposium on domain decomposition methods for partial differential equations, SIAM, pp. 1—42, 1988.

[Neu1949]
J. von Neumann, On rings of operators. Reduction theory, Ann. Math. (2), 50(2), pp. 401—485, 1949.

[Sei1874]
P. L. von Seidel, Über ein Verfahren die Gleichungen, auf welche die Methode der kleinsten Quadrate führt, sowie lineare Gleichungen überhaupt, durch successive Annäherung aufzulösen, Abh. Kgl. Bayer Akad. Wiss. Math. Phys. Kl., 11(3), pp. 81—108, 1874.

[Sch1870]
H. A. Schwarz, Ueber einen Grenzübergang durch alternirendes Verfahren, Vierteljahresschr. Naturforsch. Ges. Zürich, 15(3), pp. 272—286, 1870.

[Smo1911]
M. Smoluchowski, Über die Wechselwirkung von Kugeln, die sich in einer zähen Flüssigkeit bewegen, Bull. Int. Acad. Sci. Cracovie, Cl. Sci. Math. Nat., Sér. A Sci. Math., pp. 28—39, 1911.

[Wey1940]
H. Weyl, The method of orthogonal projection in potential theory, Duke Math. J., 7(1), pp. 411—444, 1940.

Post-scriptum :

Les auteurs tiennent à remercier chaleureusement Maxime Chupin, pour son aide et ses suggestions, les relecteurs d’Images des mathématiques Cidrolin, Denis Chadebec, Didier Roche et alainfa, dont les remarques et commentaires ont permis d’améliorer cet article, et enfin Franck Boyer, pour ses apports, son travail d’édition et sa patience.

Article édité par Franck Boyer

Notes

[1Un tel procédé itératif ne converge pas toujours. Que constate-t-on quand on l’applique à la phrase suivante « Cette phrase contient ... consonnes et ... voyelles. » ?

[2Le fait que la convergence de la méthode de Gauss-Seidel soit plus rapide que celle de la méthode de Jacobi ne signifie pas pour autant que cette dernière est moins efficace en pratique. Dans l’exemple choisi, on a besoin de dix itérations de la méthode de Gauss-Seidel (soit encore un total de trente résolutions d’équation linéaire) pour arriver à une solution approchée ayant cinq chiffres exacts après la virgule, alors que la méthode de Jacobi en nécessite vingt-quatre (c’est-à-dire un total de soixante-douze résolutions d’équation linéaire). Cependant, si l’on dispose, comme c’est aujourd’hui couramment le cas, de plusieurs unités de calcul (machines ou plus simplement processeurs), on peut faire résoudre simultanément à des unités différentes chacune des équations linéaires intervenant dans une itération de la méthode de Jacobi. On parle alors de parallélisation des calculs, expliquant le vocable « parallèle » utilisé pour qualifier la méthode de Jacobi. Une telle parallélisation est impossible dans le cas de la méthode de Gauss-Seidel, pour laquelle les résolutions au cours d’un itération s’enchaînent nécessairement. Par conséquent, en négligeant certains éléments comme le temps de communication entre unités de calcul, une itération de la méthode de Jacobi dure trois fois moins de temps qu’une itération de la méthode de Gauss-Seidel. Si l’on prend donc soin de la paralléliser, la méthode de Jacobi s’avère en pratique 20% plus rapide que la méthode de Gauss-Seidel sur le cas traité, malgré un nombre d’itérations plus que double.

[3La matrice $A$ étant supposée inversible, il existe toujours une numérotation des inconnues et des équations du système faisant que la matrice $M$ de la décomposition est elle aussi inversible pour l’un ou l’autre de ces choix.

[4On parle bien du mathématicien ayant donné son nom à une célèbre inégalité.

[5L’opérateur différentiel apparaissant dans cette équation est le laplacien, défini en coordonnées cartésiennes par $\Delta u=\sum_{i=1}^n\frac{\partial^2 u}{\partial{x_i}^2}$.

[6Le problème posé dans le carré avait en effet été résolu par Joseph Fourier en 1807 et celui posé dans le disque par Siméon Denis Poisson en 1815.

[7On aurait pu tout aussi bien chercher à déterminer les forces appliquées connaissant les vitesses des particules. On parle dans ce cas d’un problème de résistance.

[8Cette expression fait en général apparaître un développement en série, mais, dans certains cas particuliers d’intérêt (comme lorsque la particule est sphérique), ce développement ne contient qu’un nombre fini de termes.

Partager cet article

Pour citer cet article :

Guillaume Legendre, Julien Salomon — «Décomposer et itérer pour résoudre un problème» — Images des Mathématiques, CNRS, 2018

Crédits image :

img_18877 - Los Alamos National Laboratory
img_17432 - Author Unknown photographer Source scanned from Polish « Problemy » monthly, DEcember 1967, page 726
img_17433 - Author Zipfel, Ludwig, Fotograf Copyright holder ETH-Bibliothek Zürich, Bildarchiv Date and time of data generation 06:22, 25 April 2007
img_17453 - © Laboratoire de Physique des Solides
img_18334 - http://www.math.spbu.ru/analysis/history/goluzin/index.html

Commentaire sur l'article

  • Décomposer et itérer pour résoudre un problème

    le 28 décembre 2018 à 16:53, par Pierre Cami

    Décomposer et réitérer pour résoudre un problème complexe
    Nous allons nous attaquer à la conjecture de Collatz dite aussi conjecture de Syracuse ou problème 3x+1.
    On défini la transformation de Collatz de la façon suivante :
    on part d’un nombre entier positif non nul x(i), si x(i) est pair x(i-1)=x(i)/2, si x(i) est impair x(i-1)=3*x(i)+1.
    La conjecture de Collatz dit que on aura toujours une fin x(ij)=4, x(ij)=2, x(i-j)=1 puis répétition du cycle trivial 4, 2, 1.
    C’est un chercheur portugais Oliveira e Silva qui détient depuis 2009 le record de vérification de la conjecture par ordinateur, il a vérifié sa validité pour tout les nombres x(i) < 5*10^18, mais la conjecture n’a toujours pas été démontrée..
    Si la conjecture de Collatz est vérifiée on part de x(i) et on obtient après i transformations x(0)=1, mais on n’a pas à connaître la valeur de i.
    Autrement dit si on peut à partir de x(0)=1, ensemble unique, obtenir par un chemin unique tout nombre impair quelconque préalablement défini, on aura fait la preuve de la validité de la conjecture de Collatz si le chemin est suivi selon la transformation inverse de Collatz.
    On peut simplifier la transformation de Collatz sans en changer le résultat en constatant que si on part de x(i) pair il est de la forme (2*n-1)*2^k avec n et k de 1 à N et x(i-1) est impair=2*n-1, x(i-2) est pair et égal à 6*n-2 soit( 2*m-1)*2^j et ainsi de suite.
    A un nombre pair succède un nombre impair et à un nombre impair un nombre pair, donc on peut considérer uniquement les nombres impairs.
    Tout nombre impair succède à un prédécesseur impair unique ( tout comme tout nombre pair) quand on a défini le point de départ x(i).
    Si on ne défini pas de point de départ et qu’on s’intéresse aux prédécesseurs possibles d’un nombre impair il en va tout autrement.
    Un nombre impair multiple de trois n’a aucun prédécesseur impair possible, il doit donc être le premier x(i) de la suite sinon il sera le résultat de la division de x(i)=(6*n-3)*2^k par 2^k.
    Par contre tout nombre impair non multiple de trois peut avoir une infinité de prédécesseurs possibles dans une suite de Collatz.
    Prenons les nombres impairs de la forme 6*n-5 pour n de 1 à N, le plus petit prédécesseur est égal à ((6*n-5)*4^1-1)/3 et tous les autres prédécesseurs possibles sont ((6*n-5)*4^j-1)/3 pour j de 1 à N.
    Prenons les nombres impairs de la forme 6*n-1 pour n de 1 à N, le plus petit prédécesseur est égal à ((6*n-1)*2^1-1)/3 et tous les autres prédécesseurs possibles sont ((6*n-5)*2^(2*j-1)-1)/3 pour j de 1 à N.

    A(L) = 6*(L-1)+1 , L de 1 à N
    A(n) = 6*n-5 pour n de 1 à N
    T(n, j) = ((6*n-5)*4^j-1) pour n et j de 1 à N, soit le tableau suivant

    n A(n) j 1 2 3 4 5 6 7 8
    1 1 1 5 21 85 341 1365 5461 21845
    2 7 9 37 149 597 2389 9557 38229 52917
    3 13 17 69 277 1109 4437 17749 70997 283989
    4 19 25 101 405 1621 6485 25941 103765 415061
    5 25 33 133 533 2133 8533 34133 136533 546133
    6 31 41 165 661 2645 10581 42325 169301 677205
    7 37 49 197 789 3157 12629 50517 202069 808277
    8 43 57 229 917 3669 14677 58709 234837 939349
    9 49 65 261 1045 4181 16725 66901 267605 1070421
    10 55 73 293 1173 4693 18773 75093 300373 1201493
    11 61 81 325 1301 5205 20821 83285 333141 1332565
    12 67 89 357 1429 5717 22869 91477 365909 1463637
    13 73 97 389 1557 229 24917 99669 398677 1594709
    14 79 105 421 1685 6741 26965 107861 431445 1725781
    15 85 113 453 1813 7253 29013 116053 464213 1856853
    16 91 121 485 1941 7765 31061 124245 496981 1987925
    17 97 129 517 2069 8277 33109 132437 529749 2118997
    18 103 137 549 2197 8789 35157 140629 562517 2250069

    A(n) = 6*n-1 pour n de 1 à N
    T(n, j) = ((6*n-1)*2^(2*j-1)-1)/3 pour n et j de 1 à N, soit le tableau suivant :

    n A(n) j 1 2 3 4 5 6 7 8
    1 1 3 13 53 213 853 3413 13653 54613
    2 7 7 29 117 469 1877 7509 30037 120149
    3 13 11 45 181 725 2901 11605 46421 185685
    4 19 15 61 245 981 3925 15701 62805 251221
    5 25 19 77 309 1237 4949 19797 79189 316757
    6 31 23 93 373 1493 5973 23893 95573 382293
    7 37 27 109 437 1749 6997 27989 111957 447829
    8 43 31 125 501 2005 8021 32085 128341 513365
    9 49 35 141 565 2261 9045 36181 144725 578901
    10 55 39 157 629 2517 10069 40277 161109 644437
    11 61 43 173 693 2773 11093 44373 177493 709973
    12 67 47 189 757 3029 12117 48469 193877 775509
    13 73 51 205 821 3285 13141 52565 210261 841045
    14 79 55 221 885 3541 14165 56661 226645 906581
    15 85 59 237 949 3797 15189 60757 243029 972117
    16 91 63 253 1013 4053 16213 64853 259413 1037653
    17 97 67 269 1077 4309 17237 68949 275797 1103189
    18 103 71 285 1141 4565 18261 73045 292181 1168725

    On peut vérifier que les deux tableaux T(n, j) ci-dessus contiennent tous les nombres impairs une fois et une fois seulement.
    Définissons U(0) ensemble à un seul élément égal à l’unité 1.
    Définissons U(1) l’ensemble des prédécesseurs possibles de 1 dans une suite de Collatz, on trouve U(1) = (4^n-1)/3 avec N termes pour n de 1 à N.
    Les valeurs U(1) sont sur la première ligne du premier tableau ci-dessus, on remarque que 1 est présent en première position, cycle trivial oblige.
    Définissons U(2) l’ensemble des prédécesseurs possibles de U(1) dans une suite de Collatz
    Si U(1) est 1 modulo 6 alors U(2)=(U(1)*4^n-1)/3, si U(1) est 5 modulo 6 U(2)=(U(1)*2^(2*n-1)-1)/3 cela pour chaque U(1) et n de 1 à N.
    Définissons U(i+1) l’ensemble des prédécesseurs possibles de U(i) dans une suite de Collatz.
    Si U(i) est 1 modulo 6 alors U(i+1)=(U(i)*4^n-1)/3, si U(i) est 5 modulo 6 U(i+1)=(U(1)*2^(2*n-1)-1)/3 cela pour chaque U(i) et n de 1 à N.
    Tous les nombres impairs seront présent une fois et une fois seulement dans l’ensemble des ensembles U(i) pour i de 1 à N ce qui démontre la validité de la conjecture de Collatz et apporte la preuve recherchée.
    On a démontré qu’il existe une trajectoire unique pour atteindre un nombre impair quelconque en partant de 1 et il n’est pas nécessaire de connaître cette trajectoire pour en faire la preuve, la longueur i de la trajectoire elle est dépendante de la cible impaire choisie.
    Merci si vous avez lu jusqu’au bout, en attente de vos commentaires.
    Pierre CAMI

    Répondre à ce message

Laisser un commentaire

Forum sur abonnement

Pour participer à ce forum, vous devez vous enregistrer au préalable. Merci d’indiquer ci-dessous l’identifiant personnel qui vous a été fourni. Si vous n’êtes pas enregistré, vous devez vous inscrire.

Connexions’inscriremot de passe oublié ?

Suivre IDM