Les
Limites
Une deuxième raison
est la complexité des problèmes aux dérivées partielles. Pour s’en faire une
idée, réfléchissons à ce que représente la résolution numérique d’un problème
aux limites. Très schématiquement, on sera amené à résoudre un système de N
équations à N inconnues, où N est de l’ordre de grandeur de (1/h )n ; ici, h
est le pas de discrétisation (dont dépendra la précision de la solution) et n
le nombre de variables indépendantes. C’est dire que, même avec une
approximation très grossière pour un problème simple, N se comptera en
centaines pour deux variables et en milliers pour trois variables
indépendantes. On n’oubliera pas, pour estimer ces ordres de grandeur, que le
volume des calculs croît plus vite que N.
Une troisième
raison vient donner toute son importance à la deuxième; c’est dans cette
théorie que, plus encore que dans beaucoup d’autres branches des mathématiques,
les situations où on dispose de solutions explicites sont rares. Il en est de
même des solutions semi-explicites sous forme de séries.
L’analyse numérique
des équations aux dérivées partielles n’est pas née avec les ordinateurs, tant
s’en faut; la situation est plus complexe. De toute façon, l’introduction de
l’analyse numérique comme branche particulière des mathématiques est un fait
relativement récent, qui exprime son extrême spécialisation.
Un des premiers
travaux mathématiques sur les équations aux dérivées partielles – le mémoire de
Daniel Bernoulli publié en 1753 – contient deux procédés d’approximation de la
solution. L’un est celui des séries trigonométriques; l’autre consiste à
remplacer la corde vibrante par un nombre fini de masses ponctuelles reliées
par un fil élastique sans masse, ce qui donne une approximation du type
différences finies. Il est vrai que, au XVIIIe siècle, ces procédés étaient
considérés comme des outils de démonstration et non comme des méthodes
numériques. Quand Fourier, au début du XIXe siècle, reprend la méthode des
séries trigonométriques (cf. J. FOURIER, séries TRIGONOMÉTRIQUES), il
l’applique à l’équation de la chaleur qu’il vient de trouver. Il est alors
pleinement conscient de l’importance du calcul numérique des solutions que
cette méthode permet. La résolution des équations aux dérivées partielles sans
l’aide d’un ordinateur est maintenant limitée soit à des cas très particuliers,
soit à l’approximation grossière de problèmes simples.
Une autre
particularité de l’analyse numérique des équations aux dérivées partielles,
plus précisément l’approximation des problèmes d’évolution, est l’apparition de
types spécifiques d’instabilités numériques que nous examinerons dans le
chapitre 2 (Problèmes d’évolution ).
Un problème aux
dérivées partielles n’est posé de façon précise que lorsqu’on a fixé l’espace
de fonctions dans lequel on cherche la solution. Autant que possible, c’est la
nature physique du problème qui fixe cet espace et, avec lui, la norme pour
laquelle doit être évaluée et limitée l’erreur. Dans de nombreux problèmes,
l’intégrale de Dirichlet:
représente, à un
coefficient près, l’énergie du système étudié. On est alors amené à chercher
les solutions d’énergie finie , c’est-à-dire qui appartiennent à l’espace de
Sobolev H1 (cf. équations aux DÉRI- VÉES PARTIELLES - Théorie linéaire, chap.
4). Les évaluations d’erreur les plus significatives sont, en général, celles
qui concernent la norme correspondante.
1.
Problèmes aux limites stationnaires (elliptiques)
Différentes
méthodes de résolution approchée
On peut distinguer
trois familles:
– les
développements en séries de fonctions propres;
– les
approximations par différences finies ;
– les méthodes
d’éléments finis.
À la première famille
appartiennent les résolutions classiques de l’équation des cordes vibrantes et
de l’équation de la chaleur par développement en série trigonométrique de la
solution. Ce sont des méthodes puissantes qui ont été pendant longtemps un des
principaux outils de la physique mathématique. Leur grand défaut est de ne
s’appliquer qu’à des situations extrêmement spéciales, celles où on connaît les
fonctions propres de l’opérateur étudié.
Le principe des
méthodes de différences finies est simple. Presque par définition, les
expressions dérivée avancée, etdérivée retardée, sont, pour h petit, des
approximations de b u /b x. De même, le développement de Taylor montre que et
sont des approximations respectivement de cela permet d’approcher tout
opérateur aux dérivées partielles du second ordre. Si on se donne un réseau
rectangulaire, c’est-à-dire l’ensemble des points l’approximation ne fait
intervenir que les valeurs de u aux points du réseau. Notant u* jk la valeur
approchée cherchée de u (jh 1, kh 2), l’équation aux dérivées partielles est
remplacée par un système fini d’équations.
Pour fixer les
idées, considérons le problème de Dirichlet et prenons h 1 = h 2 = h . Le
système qui détermine les valeurs approchées de la solution est toujours
employées, les méthodes de différences finies ne font actuellement guère
l’objet d’études systématiques. Leur développement a probablement été stoppé
par la généralisation des méthodes d’éléments finis que nous allons examiner
maintenant.
Principe des
méthodes d’éléments finis
Commençons par un
exemple très simple, la déformation d’une tige élastique fixée à un bout et
soumise à une force longitudinale F à l’autre. Imaginons cette tige composée de
N petits ressorts accrochés bout à bout (ce sont les "éléments
finis"). Au repos, le i -ème ressort va du point x i-1 au point x i , avec
x 0 = 0 et x N = L, longueur de la tige. Appelons u i le déplacement de x i .
Nous supposons les ressorts linéaires, c’est-à-dire qu’ils suivent une loi de
Hooke, avec une constante que nous prendrons égale à k i / (x i _ x i-1 ), en
ne supposant donc pas la tige homogène. Appliquons alors le principe des
travaux virtuels à un déplacement de l’extrémité du i -ème ressort (point qui,
au repos, se trouve en x i ). On obtient ce système nous permet de calculer une
solution approchée du problème. (Dans la première de ces équations, le premier
terme est le travail effectué sur le i -ème ressort, le second sur le (i +
1)-ème ressort; on obtient l’équation en sommant élément par élément.)
Nous sommes ainsi
parvenus à une solution approchée sans tenir compte de l’équation
différentielle, vérifiée par la solution exacte. Ce point est typique de la
première phase de l’histoire des méthodes d’éléments finis qui ont surtout été
développées par des ingénieurs sur la base de considérations physiques où
l’équation aux dérivées partielles n’apparaissait pas. C’est la raison pour
laquelle ces méthodes sont restées quelque temps cachées aux mathématiciens.
Mais la situation avait ses inconvénients, en particulier un gaspillage
d’énergie du fait qu’il fallait réinventer une même méthode pour chaque champ
nouveau d’application.
Du point de vue
mathématique, les méthodes d’éléments finis sont une sous-famille des méthodes
de Ritz-Galerkin. Pour les problèmes variationnels, ces méthodes consistent à
remplacer l’espace V des fonctions admissibles par un de ses sous-espaces VN
dit "espace d’approximation ".
Si VN est de
dimension finie N et que les (fi ) en sont une base, où les coefficients u j
sont donnés par le système de N équations obtenues en faisant v = fj dans le
problème variationnel.
En pratique, on
commence en général par linéariser l’équation. Lorsqu’elle est non linéaire, on
est amené à une méthode itérative dont chaque pas est la résolution d’un
problème aux dérivées partielles linéaires. Les calculs sont évidemment
d’autant plus longs.
Dans une première
approche, on peut caractériser les méthodes d’éléments finis parmi les méthodes
de Ritz-Galerkin par le choix de cet espace VN. L’ouvert où se pose le problème
est muni d’une triangulation (ou d’un découpage en quadrangles) dont les
triangles peuvent d’ailleurs être curvilignes; ce sont les éléments. Une
fonction appartient alors à VN si sa restriction à chaque élément est d’un type
donné (en général un polynôme de degré assez petit) et si elle vérifie des
conditions de raccordement.
Dans l’exemple très
simple ci-dessus, les éléments sont les N intervalles [x i-1 , x i ]; une
fonction u appartient à VN si, sur chacun de ces intervalles, elle coïncide
avec une fonction affine et si, de plus, elle est continue (condition de
raccordement). Une telle fonction est déterminée par ses valeurs aux points x i
; ces points sont appelés des nœuds.
Une des
caractéristiques des méthodes d’éléments finis, au niveau de la mise en œuvre numérique,
est le balayage par éléments. Si on désigne par EI les éléments, une des
intégrales de la formule (4).
On calcule alors
ces intégrales pour chaque élément. Les fj sont nulles en dehors d’un petit
nombre d’éléments et il y en a donc peu à calculer. On ajoute aux équations
concernées les intégrales à mesure qu’on les calcule et on passe ensuite à
l’élément suivant.
Avantages
des méthodes d’éléments finis
Les intégrales
figurant dans (5) ne sont différentes de 0 que si les supports de fi et fj se coupent.
Il en résulte qu’un petit nombre seulement des inconnues u k figure dans chaque
équation (comme dans les méthodes de différences finies).
Le principal
avantage des méthodes d’éléments finis est leur flexibilité: on peut les
appliquer de façon quasi automatique à tout problème mis sous forme
variationnelle. On peut aussi varier le type d’éléments choisi et on a intérêt
à utiliser dans chaque problème un type d’éléments adapté. On a surtout une
très grande souplesse dans le découpage en éléments, que l’on pourra faire
beaucoup plus serré dans certaines parties du domaine si la solution y présente
des irrégularités (fig. 1).
C’est aussi un
avantage des méthodes d’éléments finis de partir d’une formulation
variationnelle, quand il en existe une, adaptée à l’aspect physique du
problème. Cela permet en particulier d’avoir un contrôle de l’erreur sous une
forme physiquement significative.
Enfin, les méthodes
d’éléments finis permettent aussi le calcul approché des premières fonctions et
valeurs propres.
Aspects
économiques
On peut écrire un
programme d’éléments finis pour la plupart des problèmes aux limites qui
interviennent dans les applications. La flexibilité de ces méthodes d’un côté
et les progrès de la technique des ordinateurs de l’autre ont permis de
développer des jeux de sous-programmes plus ou moins universels, c’est-à-dire
permettant pratiquement de résoudre presque tout problème variationnel.
De nombreux jeux de
sous-programmes ont été commercialisés. Les plus connus sont ASKA, édité par
l’université de Stuttgart sous la direction d’Argyris, et NASTRAN, produit par
la N.A.S.A. Il est bon de connaître quelques ordres de grandeur: la longueur
des programmes se compte par milliers d’instructions Fortran, leur mise au
point a demandé un millier d’hommes-années de travail et leur prix varie de 500
000 à 1 million de francs. Pour comprendre le succès de ces programmes, il faut
savoir que, dans le budget d’un centre de calcul, la rémunération des hommes
pèse maintenant beaucoup plus lourd que le matériel.
Toute médaille a
son revers. Il faut une bonne semaine à un chercheur bien formé pour commencer
à utiliser un de ces programmes géants; il est alors moins efficace qu’un
programme modeste écrit sur le tas, et il faut encore des mois pour le prendre
en main. D’autre part, quoique susceptibles de traiter les problèmes les plus
variés, les grands programmes d’éléments finis perdent beaucoup de leur
efficacité quand on s’éloigne du type d’application sur lequel ils ont été mis
au point (en général le calcul des structures).
Revenons, pour
conclure sur cette présentation sommaire des éléments finis, sur le rapport
entre l’analyse numérique des équations aux dérivées partielles et les
ordinateurs. Les problèmes posés par la résolution numérique de ce type de problème
ont plus évolué de la première génération d’ordinateurs à la situation actuelle
que de d’Alembert à cette première génération incluse. Les problèmes de volume
de mémoire, d’erreurs d’arrondi et de programmation ont constitué tour à tour
le goulet d’étranglement et le centre de l’attention, selon les changements de
la technologie et l’appréhension des techniques par les ingénieurs.
L’utilisation de gros ordinateurs effectuant des calculs en parallèle et la
construction de microsystèmes adaptés à certains problèmes d’équations aux
dérivées partielles ont permis des développements récents.
2.
Problèmes d’évolution
Généralités
On peut toujours
formuler ces problèmes de la manière suivante: Trouver une fonction u vérifiant
où u 0 est une fonction (ou une distribution) donnée et A un opérateur aux
dérivées partielles en x , complété par des conditions aux limites (problème
mixte) ou non (problème de Cauchy).
On peut, en
particulier, mettre sous cette forme le problème de Cauchy pour l’équation des
ondes:
où B est un
opérateur elliptique du second ordre; il suffit de prendre bu/bt comme fonction
inconnue auxiliaire.
Les procédés
d’approximation des problèmes aux limites ci-dessus reviennent tous à remplacer
A par un opérateur approché Ah . On est ainsi ramené au problème aux valeurs
initiales pour le système différentiel ordinaire on peut alors puiser dans les
techniques de l’analyse numérique des systèmes différentiels, en commençant par
la méthode de Runge-Kutta (cf. équations DIFFÉRENTIELLES, chap. 7) et les méthodes
linéaires multipas.
Toutefois, ces
systèmes présentent des particularités. Ils sont compliqués et à un grand
nombre de dimensions, ce qui rend malaisé l’emploi de méthodes sophistiquées; en
fait même illusoire, du fait que le système lui-même n’est qu’une
approximation.
Mais surtout, il ne
faut pas oublier que h est un paramètre "destiné à tendre vers zéro".
Ce point est source de difficultés que nous nous proposons de préciser
maintenant.
Instabilité de
discrétisations linéaires
Considérons encore
un problème très simple, que l’on peut d’ailleurs résoudre plus efficacement
par développement en série trigonométrique.
Une plaque d’un
matériau homogène est plongée dans l’eau bouillante jusqu’à obtention partout
de la température 1000, puis, à l’instant t = 0 dans de la glace fondante. La
température u ne dépend donc que du temps et de la variable d’espace
transversale à la plaque. Avec des unités de longueur et de temps adaptées, ce
problème satisfait aux équations. Approchons ce problème en divisant
l’intervalle [0, 1] en N sous-intervalles égaux et en appliquant la méthode des
différences finies par rapport à la variable d’espace et la méthode d’Euler, de
pas t, par rapport au temps. Désignant par u j, k la solution approchée au
point (j t, k /N), qui permet de calculer, sans aucune difficulté, la solution
approchée pour les valeurs successives de j .
Prenons d’abord t =
0,004 et N = 10. Au bout de dix pas de temps (t = 0,04), la solution approchée
est donnée par le tableau 1: elle est exacte à 10 près, malgré la rusticité de
l’approximation et la discontinuité dans la donnée initiale.
Dans l’espoir
d’améliorer encore cette précision, faisons maintenant N = 20; on trouve alors
le tableau 2, sur lequel on a reporté les résultats partiels des trois derniers
pas de temps. Trois constatations s’imposent à l’examen de ces résultats:
– la "solution
approchée" obtenue ne présente aucune ressemblance avec la solution
exacte;
– elle est très
grande en valeur absolue, et de plus en plus quand la variable t augmente
(instabilité);
– elle oscille,
tant en x qu’en t .
Un retour sur le
système (10) conduit à incriminer le facteur tN2 qui est passé de 0,4 à 1,6
quand on a doublé N.
La méthode d’Euler
appliquée au système (8) donne:
Si v est un vecteur
propre de Ah correspondant à une valeur propre lh , il apparaît une solution du
système discrétisé qui croît exponentiellement avec j si l’on n’a pas dans
notre exemple, cette condition donne, à un terme près qui tend vers 0 quand N
tend vers l’infini, et ce type de comportement est général dans les problèmes
paraboliques.
Ce phénomène ne se
présente pas si, au lieu de la méthode d’Euler, on emploie le schéma implicite
.On peut alors prendre des pas de temps beaucoup plus grands. Le prix à payer
pour ce gain est que, à chaque pas de temps, il faut résoudre un
"grand" système pour calculer u *j+1 . On peut montrer que toutes les
méthodes stables sont implicites; souvent, ces méthodes sont plus avantageuses
que les méthodes explicites.
3.
Analyse numérique des problèmes hyperboliques
On a vu que les
problèmes hyperboliques possèdent les propriétés suivantes:
a ) vitesse finie
de propagation;
b ) propagation des
singularités dans le cas linéaire;
c ) apparition,
dans le cas non linéaire, de singularités, interaction entre deux singularités,
propagation dans les intervalles entre ces événements.
Les méthodes
numériques relatives à ces problèmes doivent prendre en compte ces propriétés.
Le fait, en particulier, que la solution se propage a conduit à privilégier les
méthodes de différences finies par rapport aux méthodes d’éléments finis, car
on suit plus facilement l’évolution de la solution en parcourant, selon sa
propagation, les points du maillage. Il n’existe encore aucun résultat
systématique à plus d’une dimension d’espace; ou complétés par la donnée
initiale u (x , 0) = u 0(x ).
Dans (11) et (12),
u est un vecteur; dans (11), A est une matrice à valeurs propres réelles et
distinctes. Le système (12) est un système hyperbolique (cf. chap. 1 in
équations aux DÉRIVÉES PARTIELLES - Équations aux dérivées partielles non
linéaires).
Désignons par h et
t deux paramètres destinés à tendre vers zéro, et soit u n i , i ª Z et n ª N
une approximation de la solution au point (ih , n t). On remplace la dérivée
par rapport au temps par l’expression et la dérivée A(bu /bx ), ou (b/bx ) (F(u
)), par une expression de la forme. Par exemple, on peut remplacer (b/bx ) (F(u
)) .
L’équation
approchée s’écrit sous la forme (13) qui définit un schéma explicite à (2r + 1)
points.
Pour simplifier,
nous nous limiterons à un schéma à trois points. Dans ce cas, on obtient, aussi
bien pour la solution de (11) que pour celle de (12), une expression de la
forme. La fonction G devra être choisie de telle sorte que la méthode soit
stable, consistante, et, dans le cas non linéaire, conduise à une solution
respectant la condition d’entropie. Si ces trois conditions, que nous allons
expliciter et commenter, sont remplies, on peut, dans certains cas (il reste
encore beaucoup de problèmes mathématiques ouverts), prouver la convergence de
la solution approchée vers la solution de (11) ou (12).
La stabilité
signifie que la suite (u n i ) est bornée, pour une norme convenable,
indépendamment de h et de t. La consistance signifie que (14) ressemble bien à
l’équation initiale; elle se vérifie en montrant que, pour toute solution
régulière u (x , t ) de l’équation initiale, on a où a et b désignent les
ordres de la méthode. La condition d’entropie a été discutée au chapitre
Équations aux dérivées partielles non linéaires. La différence entre la notion
de stabilité et celle de consistance apparaît très bien sur la discrétisation
de l’équation sont tous deux consistants d’ordre 1 en h et t; par contre, le
schéma (15) ne peut convenir. En effet, on détermine u i n+1 à partir de u n i
et u n i+1 et donc finalement en fonction de la donnée initiale sur
l’intervalle [ih , (i + n )h ], ce qui est en contradiction avec la formule
explicite (Puisque le schéma (15) ne convient pas, c’est qu’il présente des
défauts de stabilité, que nous allons mettre en évidence en prenant une donnée
initiale égale à 0 pour x D 0 et à 1 pour x O 0.
Le schéma (15)
donne alors qui est complètement différent de la solution exacte pour a tn P ih
, mais borné, et on voit que, si le rapport t/h est fixé, u n 0 tend toujours
vers _ d avec n . Un raisonnement simple d’analyse fonctionnelle permet de
déduire de cet exemple l’existence de l’instabilité avec des données initiales
régulières. On notera (fig. 2)que, dans ce schéma, u n i est déterminé par les
données appartenant à l’intervalle [ih , (i + n )h ] qui ne contient pas le
point auquel se réduit le domaine de dépendance du point (n t, ih ).
Une condition
nécessaire de stabilité est donc liée à la direction de la propagation. Mais il
y en a une autre, liée à la vitesse de cette propagation; où les lk sont les
valeurs propres de la matrice gradU F (x) et où x parcourt l’enveloppe convexe
de l’ensemble des u n i ; dans le cas linéaire, les lk sont donc les valeurs
propres de la matrice fixe A. La condition (19) a été dégagée vers les années
1930 par Courant, Friedrichs et Lewy, d’où son nom de condition C.F.L. Nous la
supposerons désormais vérifiée. Elle ne tient pas compte de la direction de la
vitesse ni de la condition d’entropie dans le cas non linéaire. Pour remédier à
cet inconvénient, on utilise le fait que les solutions entropiques des
équations hyperboliques d’évolution (linéaires ou non) peuvent s’obtenir comme
limites de solutions du système parabolique.
On peut donc, au
lieu de schémas du type (15) ou (16), utiliser des schémas de la forme où les
coefficients a n i +1/2 et a n i _1/2 sont déterminés en fonction des u n i .
Les derniers termes de (20) simulent une viscosité, c’est-à-dire une discrétisation
de l’opérateur eb2u /bx 2. Par exemple, en prenant pour a n i +1/2 = 1/2 pour
toutes les valeurs de i et n , l’équation (20) s’écrit et on retrouve bien,
dans le second membre de (21), la discrétisation usuelle du laplacien. Ce
schéma porte le nom de Lax-Friedrichs et date des années 1950.
De même, dans le
cas scalaire, si FH(u ) est une fonction positive (dans l’exemple linéaire
(11), FH(u ) = a ), on peut choisir.Dans ce cas, l’équation (20) devient qui
n’est autre que le schéma retardé. Sous la condition C.F.L.,
ce schéma est
stable et consistant; il porte le nom de schéma de Godounov et date aussi des
années 1950.
Plus les
coefficients sont grands, meilleure est la stabilité mais plus faible est la
consistance. Cela se traduit (pour le schéma de Lax-Friedrichs par exemple) par
un lissage trop prononcé des chocs (fig. 3). La tendance actuelle est d’essayer
de développer des schémas dont les coefficients de viscosité artificielle
soient choisis de manière à n’être actifs que très près des chocs pour ne pas
détruire la précision en dehors de ceux-ci, tout en assurant la stabilité le
long de ces chocs.
Une analyse à peu
près complète de ces problèmes a été faite dans le cas scalaire, même s’il y a
plusieurs variables d’espace. Pour les systèmes hyperboliques, la situation est
beaucoup plus compliquée et a donné lieu à un développement particulièrement
intéressant que l’on va, pour terminer, détailler un peu.
Cette méthode, due
à Glimm (1966) repose sur une analyse de solutions particulières du problème.
On commence par regarder le problème de Riemann, qui est la résolution de
l’équation.
Avec une donnée
initiale constante pour x S 0 et pour x O 0 et présentant une discontinuité en
x = 0. Une interprétation physique "naïve" est la suivante: on
considère un tube infiniment long séparé en deux parties par une membrane
située au point x = 0. Dans les deux régions x S 0 et x O 0, se trouve le même
gaz, dans des états stationnaires de pression, vitesse et température homogènes.
À l’instant t = 0, on brise la membrane et on observe l’évolution du phénomène.
Mathématiquement, ce problème ressemble à la recherche des solutions
élémentaires pour les systèmes hyperboliques linéaires à coefficients
constants. Comme dans ce cas, la solution est invariante par homothétie et ne
dépend donc que de x = x /t . Comme le problème possède des propriétés de
vitesse finie de propagation, la solution reste égale à u + pour x O At et à u
- pour x S _At (où A O 0 est choisi assez grand). Dans la région |x | S At ,
elle se décompose en un nombre fini d’ondes (dites élémentaires) de choc ou de
détente. Par une méthode purement algébrique (utilisation éventuelle du
théorème des fonctions implicites), on peut, sous des hypothèses convenables,
prouver l’existence d’une solution satisfaisant la condition d’entropie. On
peut également donner un algorithme permettant de la calculer numériquement.
Par exemple, la résolution de l’équation (scalaire) de Burger avec la condition
initiale u (x , 0) = 1 si x S 0 et u (x , 0) = 0 si x O 0, conduit à une
solution on a une onde de choc qui se propage à la vitesse x Ht = 1/2. De même,
l’équation (23), avec la donnée initiale a une solution définie par les
formules suivantes on a ainsi une onde de détends on remarque ensuite que le
point de discontinuité x = 0 ne joue pas de rôle particulier; on peut, bien
entendu, résoudre un problème de Riemann avec une discontinuité au point x = M.
La solution sera alors constante sur les droites x _ M = l t . Cela conduit au
procédé d’approximation suivant. On introduit un pas de temps t et un pas
d’espace h (comme dans les méthodes de différences finies) et on construit des
solutions exactes, vérifiant la condition d’entropie, dans les bandes Rx Z [n
t, (n + 1)t]; on suppose que pour t = n t, n pair, la solution a été choisie
constante sur les intervalles ]2ih , 2(i + 1)h [, et on note u n 2i +1 sa
valeur. On résout au point 2ih le problème de Riemann avec donnée à gauche
égale à u n 2i-1 et donnée à droite égale à u n 2i+1 . Supposons maintenant que
t est choisi assez petit pour que la solution reste égale à u n 2i-1 sur
l’intervalle (2i _ 1)h Z ]n t, (n + 1)t[ et égale à u n 2i+1 sur l’intervalle
(2i + 1)h Z ]n t,(n + 1)t[. Ces conditions permettent de raccorder continûment
cette solution avec celle qui est calculée aux points 2(i _ 1)h et 2(i + 1)h .
On a obtenu ainsi une fonction u n qui est solution exacte du système et qui
vérifie la condition d’entropie. On va ensuite réitérer le procédé à l’instant
(n + 1)t (rappelons que n + 1 est impair), en partant d’une solution constante
sur les intervalles ](2i _ 1)h , (2i + 1)h [, en résolvant à nouveau des
problèmes de Riemann. Pour ce faire, il faut définir les constantes à l’instant
(n + 1)t. Cette détermination se fait en interpolant la solution obtenue à
l’instant précédent. On a donc l’algorithme suivant, décrit à partir de
l’instant n t, avec n pair pour fixer les idées (fig. 5).
1. Résolution de
problèmes de Riemann aux points 2ih à l’instant n t.
2. Interpolation de
la solution obtenue à l’instant (n + 1)t par des fonctions constantes sur les
intervalles ](2i _ 1)h , (2i + 1)h [.
3. Résolution de
problèmes de Riemann à l’instant (n + 1)t aux points (2i + 1)h .
4. Interpolation de
la solution obtenue à l’instant (n + 2)t par des fonctions constantes sur ]2ih
, 2(i + 1)h [.
5. Revenir en 1,
avec n augmenté de 2 (fig. 5).
Il faut maintenant
examiner la stabilité et la consistance de ce procédé.
La stabilité est
facile pour une équation scalaire. Par contre, elle n’a été démontrée pour les
systèmes qu’en 1966, par Glimm. Dans le cas des systèmes, la démonstration
suppose de plus qu’il s’agit d’une "vraie" interpolation,
c’est-à-dire que la fonction interpolée prend effectivement la valeur constante
d’interpolation (nous verrons que ce n’est pas toujours le cas).
La démonstration de
Glimm présente un grand intérêt. Il prouve que la solution approchée est
uniformément bornée dans les espaces de fonctions à variation bornée, en
faisant une analyse très fine de l’interaction entre les ondes engendrées par
les différents problèmes de Riemann. Cette interaction se mesure en
introduisant une fonctionnelle qui tient compte à la fois des interactions
entre ondes proches et ondes lointaines. Cette démarche ressemble – bien
qu’aucun lien rigoureux ne soit fait à ce niveau – à l’analyse donnée par
Kolmogorov (1941) du spectre de la turbulence (et des équations de
Navier-Stokes dans les problèmes non linéaires).
Une fois obtenue la
stabilité, il convient maintenant de regarder la consistance . Elle est visiblement
réalisée dans les bandes ouvertes ]n , (n + 1)t[, puisqu’on y a pris la
solution exacte. Mais on aurait tort d’en conclure que la méthode est toujours
consistante.
Pour s’en
convaincre, il suffit d’appliquer la méthode de Glimm à l’exemple en choisissant
t/h = l S 1 et en interpolant chaque fois par la valeur de la solution aux
points 2ih ou (2i + 1)h selon que l’on résout à un temps correspondant à n pair
ou impair. Un calcul direct montre que l’on obtient alors une solution qui,
pour h tendant vers 0, t/h = l, converge vers la fonction cette fonction n’a
rien à voir avec la solution cherchée, ne serait-ce qu’à cause de la présence
de l.
Pour surmonter
cette difficulté, Glimm a eu l’idée de choisir les points d’interpolation de
manière aléatoire, pour obtenir la consistance de la méthode. Cette idée a été
développée ensuite numériquement par Chorin (1976) qui a remarqué qu’il
suffisait de choisir comme interpolé u k h (n t, 2ih + jn h ), où (jn ) est une
suite déterministe de points équidistribués dans l’intervalle ] _ 1, 1[. Les
résultats expérimentaux de Chorin furent ensuite complétés par une
démonstration de Liu. Enfin, Colella et Chorin déterminèrent un choix optimal
pour la suite (jn ). C’est un bel exemple d’interaction entre la démonstration
de théorèmes et leur mise en œuvre numérique. On peut, pour terminer, remarquer
que la méthode de Glimm ne diffère des schémas aux différences de
Lax-Friedrichs et Godounov que par le procédé d’interpolation. Ainsi, pour
obtenir le schéma de Lax-Friedrichs, on choisit à l’instant (n + 1)t, n pair,
la valeur moyenne qui n’est pas obligatoirement une valeur de la fonction,
celle-ci étant discontinue.
Comme le second
membre de (25) est l’intégrale d’une solution exacte de
dans la bande ]n t,
(n + 1)t[ Z Rx , on peut intégrer cette équation (26) sur le rectangle ](2i _
1)h , (2i + 1)h [ Z ]n t, (n + 1)t[ (fig. 5), grâce à la vitesse finie de
propagation; la fonction u n (x , t ) reste constante sur les côtés (2i _ 1)h Z
]n t, (n + 1)t[ et (2i + 1)h Z ]n t, (n + 1)t[. On obtient ainsi la relation ce
qui est bien le schéma de Lax-Friedrichs (fig. 5). L’équation sous la forme
(27) permet le calcul aux points pairs du maillage à partir des points impairs
mais on peut, bien entendu, compléter cette équation par celle qui est obtenue
en calculant les points d’indice impair à partir des points d’indice pair.
Cette démarche a l’avantage de ne pas différencier dans le calcul les valeurs
approchées aux points 2ih et (2i + 1)h . Comme il a déjà été dit, on ne sait
pas, pour les systèmes, prouver la stabilité des schémas du type (28), mais les
expériences numériques donnent, pour les systèmes, des résultats de stabilité
et de consistance analogues à ceux qui sont établis pour les équations.
De façon générale,
on peut dire qu’il n’existe pas, pour les systèmes hyperboliques, de méthode
générale, contrairement aux situations elliptique ou parabolique. Cette
différence se traduit par l’absence de codes comparables à NASTRAN et autres.
Les propriétés particulières de chaque problème jouent un rôle essentiel.
L’interaction entre la réflexion sur la structure du problème et les essais
numériques est fondamentale. Les travaux s’orientent vers l’obtention de
schémas à haute précision en dehors des chocs et vers des méthodes d’analyse de
certaines solutions particulières de problèmes pluridimensionnels; il s’agit de
comprendre et de calculer la progression d’une onde de choc en dimension
supérieure à un.
On peut conclure
cette présentation de l’analyse numérique des équations aux dérivées partielles
en insistant sur la fécondité de l’interaction entre l’étude numérique et le
travail plus purement théorique. Cette fécondité a été vivifiée par
l’apparition des ordinateurs, puis par la puissance qu’ils ont acquise.