Méthodes numériques appliquées à la conception par

éléments finis
David Dureisseix

To cite this version:

David Dureisseix. Méthodes numériques appliquées à la conception par éléments finis. DEA. ENS
cachan - Université de Montpellier 2, 2008, pp.98. �cel-00351713v1�

HAL Id: cel-00351713

[Link]
Submitted on 10 Jan 2009 (v1), last revised 3 Jan 2014 (v2)

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
 Méthodes numériques appliquées
à la conception par éléments finis

David Dureisseix

28/09/2008
 TABLE DES MATIÈRES

INTRODUCTION 1

1 RÉSOLUTION DES SYSTÈMES LINÉAIRES 5

1.1 Conditionnement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.1.1 Définition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.1.2 Influence des perturbations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Techniques de stockage et coût . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2.1 Paramètres influents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2.2 Coût . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3 Méthodes directes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3.1 Méthode d’élimination de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3.2 Factorisation de Cholesky . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3.3 Factorisation de Crout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3.4 Lien avec la condensation de Schur . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3.5 Autres techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.4 Quelques particularités en éléments finis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.4.1 Conditions aux limites en déplacement . . . . . . . . . . . . . . . . . . . . . . . . 14
1.4.2 Relations linéaires entre degrés de liberté . . . . . . . . . . . . . . . . . . . . . . 16
1.4.3 Système singulier : semi défini positif . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.5 Méthodes itératives stationnaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.5.1 Méthode de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.5.2 Méthode de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.5.3 Méthode de relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.5.4 Méthode d’Uzawa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.6 Méthodes itératives non-stationnaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.6.1 Méthode générique de projection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.6.2 Méthode du gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.6.3 Méthode du gradient conjugué . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.6.4 Préconditionnement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
1.6.5 Autres méthodes itératives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
1.7 Critères d’arrêt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2 SYSTÈMES AUX VALEURS PROPRES 25

2.1 Préliminaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.1.1 Rappels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.1.2 Domaine d’emploi des différentes méthodes . . . . . . . . . . . . . . . . . . . . 26
2.2 Méthode de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.3 Méthode des puissances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.4 Méthode des itérations inverses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.5 Méthode des itérations inverses avec décalage spectral . . . . . . . . . . . . . . . . . . 30
2.5.1 Principe et convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.5.2 Exemples d’utilisation dans un code éléments finis . . . . . . . . . . . . . . . 31
2.6 Transformation de Householder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.6.1 Factorisation QR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.6.2 Détermination des valeurs propres par la méthode QR . . . . . . . . . . . . . 32
2.6.3 Forme de Hessenberg . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
 2.7 Exemple d’utilisation dans un code éléments finis . . . . . . . . . . . . . . . . . . . . . 34
2.7.1 Réduction dans un sous-espace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.7.2 Technique itérative . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.8 Méthode de Lanczos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3 RÉSOLUTION DE SYSTÈMES NON-LINÉAIRES PARTICULIERS 37

3.1 Non-linéarités matérielles de type plasticité . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.1.1 Formulation du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.1.2 Résolution pour un incrément de temps . . . . . . . . . . . . . . . . . . . . . . . 38
3.1.3 Méthode de la plus grande pente . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.1.4 Méthodes de gradient conjugué . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.1.5 Méthode de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.1.6 Méthode de Newton-Raphson modifiée . . . . . . . . . . . . . . . . . . . . . . . 39
3.1.7 Mises à jour de type Quasi-Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

4 RÉSOLUTION DE SYSTÈMES DIFFÉRENTIELS 43

4.1 Équation scalaire du premier ordre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.1.1 Problème de Cauchy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.1.2 Ré-écriture du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.1.3 Schémas explicites à un pas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.1.4 Quelques problèmes pratiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.1.5 A-stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.1.6 Méthode des trapèzes généralisée . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.1.7 Exemple de méthode à pas multiples : Adams . . . . . . . . . . . . . . . . . . . 48
4.1.8 Exemple de méthode de prédiction correction . . . . . . . . . . . . . . . . . . . 50
4.1.9 Méthode de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.1.10 Autres méthodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.2 Systèmes différentiels du premier ordre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.2.1 Emploi des méthodes précédentes . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.2.2 Choix d’un schéma : influence de l’ordre d’intégration . . . . . . . . . . . . . 54
4.3 Systèmes différentiels du deuxième ordre . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.3.1 Méthode de Newmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.3.2 Méthode de Gear implicite à deux pas . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.3.3 θ-méthode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.3.4 α-HHT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.3.5 Autres méthodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.4 Intégration de modèles de comportement . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.4.1 Exemple de modèle de plasticité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.4.2 Méthodes de type retour radial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.4.3 Exemple de modèle de viscoplasticité . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.5 Cas de la dynamique non régulière . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

RÉFÉRENCES 69

A FACTORISATION DE CHOLESKY - ORTHOGONALISATION D’UN SOUS-ESPACE 73

B BORNAGE PAR LES VALEURS PROPRES ÉLÉMENTAIRES 75

B.1 Système aux valeurs propres généralisé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
B.2 Système aux valeurs propres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

C TECHNIQUE D’ACCÉLÉRATION D’AITKEN 77

D REPÈRES BIOGRAPHIQUES 79

INDEX 85
 INTRODUCTION

Ce document de travail (avec renvoi à des ressources de références s’il y a besoin d’approfon-
dir) est utile à ceux qui sont intéressés ou amenés à utiliser des méthodes de simulation (es-
sentiellement par éléments finis) dans le but de dimensionner des structures mécaniques.
Il ne s’agit pas pour autant d’un document sur les éléments finis, puisqu’on suppose que
les problèmes ont déjà été formulés dans ce cadre, et que l’on se place en aval, lorsque la
résolution des problèmes éléments finis est en jeu. On se ramènera ainsi fréquemment aux
techniques effectivement utilisées dans des codes de calcul industriels. Pour cela, on sera
amené à utiliser des méthodes classiques d’algèbre linéaire, de calcul matriciel, adaptées
aux calculs sur ordinateurs (algorithmique).
Selon les problèmes qui se posent, on est amené à utiliser un grand nombre d’outils dif-
férents mais ce document ne se veut pas exhaustif. Il ne traite en particulier pas les cas sui-
vants :
– la recherche de racine de f (U) = 0 ;
– les problèmes de minimum non quadratiques (bisection, point fixe, Newton, Newton-
Raphson, quasi-Newton, sécante, BFGS) ;
– les problèmes non linéaires de type grands déplacement ou contact unilatéral.
Structuré en quatre chapitres, il aborde les techniques de résolution de grands systèmes
linéaires (de façon directe ou itérative), de problèmes de recherche de valeurs et vecteurs
propres, et enfin de systèmes différentiels.
Pour un aperçu moins détaillé mais plus large des méthodes numériques, le lecteur pourra
se reporter à [1].
Il semble aussi utilise de préciser que ce document est en perpétuelle évolution et que
de ce fait, son organisation ainsi que sa table des matières peuvent afficher quelque incohé-
rence.

En cas d’erreur
Les différents algorithmes et techniques présentés ne sont pas garantis contre les erreurs et
autres coquilles. . . De la même façon, tout document est perfectible, et celui-ci ne fait pas
exception à la règle. Aussi, si vous en trouvez, signalez-les : c’est aussi une méthode itéra-
tive. . . pour réduire le nombre de bugs. J’ai déjà eu des retours qui m’ont permis d’améliorer
la chose. Merci donc en particulier à Martin Genet.
Tout lecteur ayant des remarques et propositions constructives est amicalement invité à
envoyer un retour à l’auteur par courriel à [Link]@[Link].
2 INTRODUCTION

Notations
Un vecteur est représenté par une lettre romane minuscule (comme x), une matrice par une
lettre romane majuscule (comme A), un scalaire par une lettre grecque minuscule (comme
α). La composante i du vecteur x est notée x i (c’est un scalaire) et la composante i , j de la
matrice A, Ai ,j : c’est un scalaire également.
La colonne j de A, par exemple, sera classiquement notée A1:n,j . Les vecteurs de la base
dans laquelle les différentes matrices sont écrites seront notés ei .
À chaque fois, l’itéré numéro i est désigné par une mise en exposant (comme x(i ) ) et
en général, une même lettre désignera en fait un même emplacement en mémoire ; par
exemple, x(i ) écrase x(i −1) en mémoire.

Pour approfondir
Pour des renseignements généraux complémentaires sur les techniques employées, il est
utile de ce référer par exemple à [2, 3]. . . et pour les algorithmes et leur implantation, aux
librairies de netlib. . .
De façon générale, la bibliographie présentée regroupe deux types de références : les
anciennes, voire très anciennes, en tant que base historique, et les plus récentes, en tant
qu’ouvrages contemporains. Les biographies utilisées se sont largement inspirées des bases
de données disponibles aux adresses suivantes :
– [Link]
– [Link]
– [Link]

Cadre de travail
Actuellement, et pour ce qui nous intéresse, nous avons tout intérêt à utiliser des librairies
scientifiques existantes de façon intensive : on ne peut prétendre à reprogrammer de façon
plus efficace et plus robuste des algorithmes standards que des spécialistes de ce domaine
et qui travaillent sur ces librairies depuis un grand nombre d’années. Il est cependant néces-
saire de connaître le principe d’un certain nombre de ces méthodes, leur domaine d’emploi,
leur limite et d’avoir une idée des cas à problème pour aider au choix d’une méthode. C’est
dans cette direction que se place ce document.
Parmi ces librairies, en grande partie du domaine public, certaines émergent ou com-
mencent à émerger comme des standards de fait. Par exemple, la librairie LAPACK (Linear
Algebra PACKage), [4] :
[Link]
qui succède en grande partie aux anciennes libraries EISPACK, pour la résolution de sys-
tèmes aux valeurs propres :
[Link]
et LINPACK pour la résolution de systèmes linéaires :
[Link]
 3

Sa description est celle de la FAQ (Frequently Asked Questions) :

LAPACK provides routines for solving systems of simultaneous linear equations,
least-squares solutions of linear systems of equations, eigenvalue problems, and
singular value problems. The associated matrix factorizations (LU, Cholesky, QR,
SVD, Schur, generalized Schur) are also provided, as are related computations
such as reordering of the Schur factorizations and estimating condition numbers.
Dense and banded matrices are handled, but not general sparse matrices. In all
areas, similar functionality is provided for real and complex matrices, in both
single and double precision.
Citons aussi sa version parallélisée, ScaLAPACK (Scalable LAPACK, [5]) :
The ScaLAPACK project is a collaborative effort involving several institutions (Oak
Ridge National Laboratory, Rice University, University of California, Berkeley, Uni-
versity of California, Los Angeles, University of Illinois, University of Tennessee,
Knoxville), and comprises four components : dense and band matrix software
(ScaLAPACK), large sparse eigenvalue software (PARPACK and ARPACK), sparse
direct systems software (CAPSS and MFACT), preconditioners for large sparse
iterative solvers (ParPre). The ScaLAPACK (or Scalable LAPACK) library includes
a subset of LAPACK routines redesigned for distributed memory MIMD parallel
computers. It is currently written in a Single-Program-Multiple-Data style using
explicit message passing for interprocessor communication. It assumes matrices
are laid out in a two-dimensional block cyclic decomposition. Like LAPACK, the
ScaLAPACK routines are based on block-partitioned algorithms in order to mi-
nimize the frequency of data movement between different levels of the memory
hierarchy. (For such machines, the memory hierarchy includes the off-processor
memory of other processors, in addition to the hierarchy of registers, cache, and
local memory on each processor.) The fundamental building blocks of the Sca-
LAPACK library are distributed memory versions (PBLAS) of the Level 1, 2 and
3 BLAS, and a set of Basic Linear Algebra Communication Subprograms (BLACS)
for communication tasks that arise frequently in parallel linear algebra computa-
tions. In the ScaLAPACK routines, all interprocessor communication occurs wi-
thin the PBLAS and the BLACS. One of the design goals of ScaLAPACK was to have
the ScaLAPACK routines resemble their LAPACK equivalents as much as possible.
D’un niveau d’utilisation plus bas, et utilisées dans les librairies précédentes, on trouve aussi
des utilitaires intéressants dans les librairies ARPACK :
[Link]
ARPACK is a collection of Fortran77 subroutines designed to solve large scale ei-
genvalue problems. The package is designed to compute a few eigenvalues and
corresponding eigenvectors of a general n × n matrix A. It is most appropriate
for large sparse or structured matrices A where structured means that a matrix-
vector product w ← Av requires order n rather than the usual order n 2 floating
4 INTRODUCTION

point operations. This software is based upon an algorithmic variant of the Ar-
noldi process called the Implicitly Restarted Arnoldi Method (IRAM). When the
matrix A is symmetric it reduces to a variant of the Lanczos process called the
Implicitly Restarted Lanczos Method (IRLM). These variants may be viewed as
a synthesis of the Arnoldi/Lanczos process with the Implicitly Shifted QR tech-
nique that is suitable for large scale problems. For many standard problems, a
matrix factorization is not required. Only the action of the matrix on a vector is
needed.

et BLAS (Basic Linear Algebra Subprograms) :

[Link]

The BLAS (Basic Linear Algebra Subprograms) are high quality “building blo-
ck” routines for performing basic vector and matrix operations. Level 1 BLAS do
vector-vector operations, Level 2 BLAS do matrix-vector operations, and Level 3
BLAS do matrix-matrix operations. Because the BLAS are efficient, portable, and
widely available, they are commonly used in the development of high quality li-
near algebra software, LINPACK and LAPACK for example.

Les développements concernant ces librairies sont toujours en cours, en particulier les ver-
sions creuses SPBLAS (SParse BLAS) et parallèles PBLAS (Parallel BLAS).
Enfin, des librairies scientifiques généralistes sont aussi disponibles, comme la GSL (GNU
Scientific Library) :
[Link]
 1
RÉSOLUTION DES
SYSTÈMES LINÉAIRES

On est souvent, voire tout le temps, amené à résoudre des systèmes linéaires de la forme
Ax = b pour la simulation sur ordinateur ; par exemple en éléments finis [6], une analyse
statique conduit à résoudre un système que l’on notera Ku = f où K est la matrice de rigidité,
u le déplacement en certains points de la structure et f, les forces généralisées.
Une analyse non-linéaire implicite, la recherche d’une base modale. . . sont souvent ba-
sées sur la brique de base qu’est la résolution de systèmes linéaires. Dans ces cas de figure,
comme dans d’autres, les matrices A et b correspondantes ont souvent des particularités
dont il sera judicieux de tirer parti pour les avantages, ou qu’il faudra avoir en tête pour les
inconvénients. Par exemple :

1. une grande dimension dim A = n, classiquement de 103 à 106 , donc :

– on n’inverse pas A (ce serait résoudre n systèmes linéaires) ;
– il faut prendre en compte des considérations de stockage adapté, à la fois pour ré-
duire les ressources en mémoire et le nombre d’opérations à effectuer.

2. la nécessité de traiter plusieurs second membres :

– simultanément ;
– successivement.

3. des propriétés souvent rencontrées :

– des coefficients réels, une matrice A carrée,
– A inversible ⇔ A régulière ⇔ det A 6= 0
– A symétrique ⇔ A = AT avec ∀x, y, x · y = xT y, xT Ay = yT AT x
– A positive ⇔ ∀x, xT Ax ¾ 0
– A définie ⇔ xT Ax = 0 ⇒ x = 0
– A mal conditionnée, typiquement cond K = O ( h12 ) où h est la taille de la maille en
éléments finis massifs.

Pour plus d’informations concernant ces particularités et les techniques mises en œuvre, on
pourra se référer à [7, 8].
 6 RÉSOLUTION DES SYSTÈMES LINÉAIRES

1.1 Conditionnement
1.1.1 Définition
Pour A symétrique et inversible,

cond A = kAk · kA−1 k (1.1)

Avec une norme matricielle cohérente avec celle définie pour les vecteurs :

kAxk
kAk = sup (1.2)
kxk6=0 kxk
P
où kxk2 = xT x = i x i2 est la norme euclidienne, on a :

|λ|max
cond A = (1.3)
|λ|min

où les λ sont les valeurs propres de A.

En général, il est très coûteux d’atteindre le conditionnement de A, par contre un certains
nombre de techniques peuvent permettre de l’estimer. On pourra en particulier s’inspirer
des techniques développées dans le chapitre concernant la recherche de valeurs propres.
Un conditionnement grand conduit à des problèmes de précision (cumul d’arrondis numé-
riques). Pour les illustrer, on est amené à considérer un système perturbé (par des perturba-
tions généralement supposées petites) Ãx̃ = b̃.

Exemple 1
Considérons le système Ax = b suivant :
  ! ! !
1 1 x1 2 2
  = ⇒ x= (1.4)
1 1+ǫ x2 2 0

Le polynôme caractéristique λ2 − (2 + ǫ)λ + ǫ = 0 permet de déterminer les valeurs propres,

à savoir :
p
2+ǫ± 4 + ǫ2
λ1,2 = (1.5)
2

Avec ǫ petit, le conditionnement est donc très grand puisque cond A ≈ 4/ǫ. Considérons
maintenant le système Ãx̃ = b̃ perturbé tel que :
  ! ! !
1 1 x1 2 2 − αǫ
  = ⇒ x= (1.6)
α
1 1+ǫ x2 2+α ǫ

Ce résultat entraîne la discussion suivante :

− si α ≪ ǫ ≪ 1 alors x ≈ [2 0]T ;

− si α ≈ ǫ ≪ 1, ce qui est tout à fait possible avec les arrondis numériques, la solution
est très différente puisque x ≈ [1 1]T .
 1.1 Conditionnement 7

Exemple 2
Examinons le système suivant :
  ! !
1,2969 0,8648 0,8642 0,8642
A=  b= x= (1.7)
0,2161 0,1441 0,1440 0,1440

Le conditionnement de A est cond A = 3,3·108 ; le résidu est petit puisque b − Ax = [−10−8 10−8 ]T
alors que la solution exacte est très différente de x puisque c’est [2 − 2]T . Il semble donc que
le conditionnement du système soit très fortement lié à la qualité de la solution (en terme
d’influence des arrondis numériques) que l’on puisse espérer avoir.
Dans le cas des éléments finis massifs (autre que plaques en flexion, poutres. . . ), le condi-
tionnement varie en O ( h12 ). La figure 1.1 illustre cela dans un cas 2D très simple. Le condi-
tionnement a été obtenu en calculant explicitement la plus grande et la plus petite valeur
propre de K avec les méthodes qui seront présentées ultérieurement.

105

L
h
104
cond K

103

102

10
10−2 10−1 1
h/L
(a) maillage (b) conditionnement

fig. 1.1 - Conditionnement en éléments finis

1.1.2 Influence des perturbations

Considérons le système Ax̃ = b̃. Si b̃ = b + δb avec kδbk ≪ kbk, a-t-on kδxk ≪ kxk pour
x̃ = x + δx ? Calculons simplement, en utilisant l"inégalité de Milkowski :

Ax = b ⇒ kAx = bk ⇒ kAk · kxk ¾ kbk (1.8)

autrement dit :

1 kAk
¶ (1.9)
kxk kbk
et maintenant :

A(x + δx) = b + δb ⇒ Aδx = δb ⇒ δx = A−1 δb ⇒ kδxk ¶ kA−1 k · kδbk (1.10)

 8 RÉSOLUTION DES SYSTÈMES LINÉAIRES

soit :

kδxk ¶ kA−1 k · kδbk (1.11)

par conséquent :
kδxk kδbk
¶ cond A (1.12)
kxk kbk
Considérons maintenant le système perturbé Ãx̃ = b. Si à = A + δA et x̃ = x + δx, alors :

(A + δA)(x + δx) = b ⇒ Aδx + δA(x + δx) = 0 ⇒ δx = −A−1 δA(x + δx)

(1.13)
⇒ kδxk ¶ kA−1 k · kδAk · kx + δxk
autrement dit :
kδxk kδAk
¶ cond A (1.14)
kx + δxk kAk
Il n’est donc pas étonnant qu’on se trouve face à une difficulté lorsque le système à résoudre
est très mal conditionné, autant avec les méthodes directes qui fournissent une solution
erronée qu’avec les méthodes itératives qui ne convergeront pas.

1.2 Techniques de stockage et coût

1.2.1 Paramètres influents
Dans les applications envisagées, en particulier pour les systèmes de grande taille, il est né-
cessaire d’envisager des techniques de stockage adaptées. Par exemple, si la taille du sys-
tème est n = 103 , le stockage de A pleine nécessite 8 Mo de mémoire ; avec n = 106 , 8 000 Go
sont nécessaires.
Pire encore, le nombre d’opérations pour résoudre le système Ax = b est lié au nombre
de termes non nuls (en fait au nombre de termes stockés, donc supposés non nuls) de la
matrice, et aussi de la numérotation des inconnues, comme on le verra. Par exemple, une
méthode de Cramer sur une matrice pleine demande un nombre d’opérations en virgule
flottante de n op ≈ (n +1)!, ce qui devient rapidement inextricable, même avec les ordinateurs
les plus puissants [9, 10].
Dans le cas des systèmes issus d’une discrétisation par éléments finis, le système est dit
creux 1 . Il faut donc en tenir compte au niveau du stockage, en utilisant par exemple un
stockage :
– bande (largeur de bande, bandwidth, l b ) ;
– ligne de ciel, skyline ;
– morse ;
– semi-morse.
Quant à l’influence de la renumérotation des degrés de liberté (ddl), des méthodes heuris-
tiques sont employées car la détermination de la numérotation optimale est un problème
extrêmement coûteux (en fait, plus coûteux que le système que l’on cherche à résoudre). Les
techniques les plus employées se nomment :
1. en anglais, un système creux est dit sparse
 1.2 Techniques de stockage et coût 9

– Cuthill-MacKee inverse [11] ;

– nested dissection [12] ;
– degré minimum [13, 14]. . .

1.2.2 Coût
Le coût d’un algorithme s’estime de différentes manières. Une première indication est la
taille de la mémoire requise pour traiter un problème donné. Une autre est le temps mis
pour effectuer la résolution du problème. Deux types de prise de temps sont généralement
employés : le temps CPU (Central Processing Unit), lié aux opérations effectuées, et le temps
horloge, lié au temps extérieur écoulé pendant cette résolution. Il est certain que ce dernier
est une mesure très globale puisqu’elle fait intervenir la machine utilisée 2 , les compétences
du programmeur, et peuvent dépendre fortement du problème traité. Elle ne caractérise
donc pas uniquement l’algorithme mais plutôt globalement le logiciel réalisé pour un pro-
blème donné sur une machine donnée.
D’autres indicateurs sont également disponibles, comme le nombre d’opérations en vir-
gule flottante nécessaire à la résolution. C’est ce que l’on appelle la complexité de l’algo-
rithme. Elle devient plus liée à l’algorithme même en faisant des estimations sur les types
de problèmes à traiter, mais elle ne fait pas intervenir les éventuels transferts en mémoire,
opérations sur les entiers. . .
La vraie qualification d’une technique doit donc faire appel à tous ces indicateurs. Ici,
seule la complexité sera présentée par la suite. L’estimation du coût effectif en temps CPU
des différentes opérations en virgule flottante peut ensuite être proposée par des essais nu-
mériques et des prises de temps sur des architectures de machines cibles.
À titre d’indication, un petit test a été réalisé dans ce sens, qui fait intervenir essentielle-
ment les performances du processeur. Les coûts indiqués dans le tableau 1.1 ont été obtenus
en janvier 2000. Entre parenthèses, on trouvera les coûts obtenus l’année ou les deux années
précédentes. Le code est écrit en Fortran 77 ; les différentes machines testées ainsi que les
options de compilation sont les suivantes :
– fronsac est une station de travail HP du LMT-Cachan (HP 9000/879, système B.10.20),
la compilation est faite avec f77 +O2 ;
– bacchus est une machine parallèle SGI 02000 du Pôle Parallélisme Île de France Sud
(SGI IP27 système 6.5), et la compilation, f90 -O2 ;
– artemis est une machine IBM SP de l’Antenne de Bretagne de l’ENS de Cachan (IBM
AIX système 3), et la compilation, f77 -O2.
Les opérations dénommées daxpy (série de 2 opérations en fait : addition et multiplication),
dscal (série de multiplications), dgemv (calcul d’un produit matrice A par vecteur x, pon-
déré par un coefficient β et ajouté à un précédent vecteur y, lui même multiplié par un co-
efficient α) sont appelées avec la librairie BLAS, déjà mentionnée. L’opération dénommée
αy + βAx est donc l’équivalent de dgemv, programmé directement en Fortran. À partir de

2. les performances du système, du compilateur, les accès aux disques

 10 RÉSOLUTION DES SYSTÈMES LINÉAIRES

type d’opération coût (10−8 s CPU)

(moyennée sur 108 termes) fronsac bacchus artemis
addition 5,2 (5,7 ; 7) 5,54 (6,2) 2,67
multiplication 5,2 (5,7 ; 7) 5,54 (6,2) 2,67
division 9 9 (11,7) 2,73
racine carrée 9 14,5 (18,3) 18,2
daxpy 8,9 (9,6 ; 10) 4,66 (5,5) 2,02
dscal 4,97 (5,3 ; 6) 2,56 (3,9) 1,35
type d’opération (moyennée sur coût (10−2 s CPU)
100 vecteurs de 1000 termes) fronsac bacchus artemis
αy + βAx 4,95 (5,2) 3,4 (3,9) 0,84
dgemv - 2 (2,3) 1,1

tab. 1.1 - estimation du coût des opérations élémentaires

ces quelques résultats, on peut tirer quelques conclusions : outre le fait qu’utiliser les op-
tions d’optimisation des compilateurs est intéressante, l’utilisation de librairies (ici la BLAS)
est plus efficace que reprogrammer des routines élémentaires soi-même. En fait, l’utilisation
de routines standards augmente la lisibilité du programme, permet de programmer à un ni-
veau plus élevé et donc réduit les temps de « debuggage », et permet souvent d’atteindre de
meilleures performances, puisque ces librairies sont programmées de façon plus efficace
(souvent directement en assembleur) et sont optimisées pour la machine sur laquelle elles
sont utilisées.
D’autre part, pour l’utilisateur normal, le coût effectif (en secondes CPU) dépend d’un
grand nombre de paramètres comme l’efficacité du compilateur, de l’environnement de dé-
veloppement, des temps de latence de la mémoire, des caches, de l’accès aux disques, etc.
En particulier, le changement de version du système d’exploitation d’une année sur l’autre
influe sur les résultats du test présenté précédemment.

1.3 Méthodes directes

Encore majoritairement employée dans les codes actuels pour leur robustesse, et d’une pré-
vision fine du coût (nombre d’opérations, stockage nécessaire), elles deviennent coûteuses
en temps et en encombrement lorsque le système à résoudre devient très grand. Pour plus
d’informations sur ces méthodes, on pourra avantageusement consulter [15].

1.3.1 Méthode d’élimination de Gauss

Principe
On simplifie le système pour le rendre triangulaire supérieur par utilisation de combinaison
linéaires de lignes (c’est-à-dire d’équations).

Préliminaire : élimination de Gauss par bloc

Écrivons le système à résoudre sous une forme de blocs :
  ! !
A11 A12 x1 b1
  = (1.15)
A21 A22 x2 b2
1.3 Méthodes directes 11

Si A11 est inversible, la première équation conduit à exprimer x1 en fonction de x2 :

x1 = A−1
11 (b1 − A12 x2 ) (1.16)

En utilisant cette expression dans la deuxième équation, on se ramène à un système avec x2

comme seule inconnue :

A22 − A21 A−1 A12 x2 = b2 − A21 A−1 b1 (1.17)

| {z 11 } | {z 11 }
S c
S est appelé complément de Schur ou matrice condensée sur les degrés de liberté x2 . Connais-
sant x2 , on remonte à x1 avec (1.16). Le système de départ (1.15) est équivalent à :
  ! !
A11 A12 x1 b1
  = (1.18)
0 S x2 c
et on peut le réécrire :
   !   !
I 0 A A12 x1 I 0 b1
   11  =  (1.19)
−1
A21 A11 I 0 S x2 A21 A−1
11 I c
| {z }
A
Utilisation pour la factorisation LU
Si on prend comme bloc 11 la première équation scalaire du système, il faut avoir A11 6= 0
(c’est le premier pivot). En continuant sous la forme de l’équation (1.19), récursivement sur
les sous-matrices qui restent, on obtient la factorisation A = LU (en n − 1 étapes) avec U
(pour upper), matrice triangulaire supérieure et L (pour lower), matrice triangulaire infé-
rieure avec diagonale unitaire.
Pour des raisons d’efficacité de stockage, L et U écrasent en mémoire A si on ne stocke
pas la diagonale de L qu’on sait unitaire. À chaque itération k , la matrice de départ A est
donc modifiée en une matrice A(k ) .

Conditions d’emploi
(k −1)
Si à chaque étape de la factorisation, le pivot Ak k est non nul, on peut progresser sans pro-
blème. La seule condition est que A soit inversible, mais quand on trouve un pivot nul, il faut
procéder à des échanges d’équations. En pratique, on a bien entendu aussi des problèmes
si le pivot est petit en valeur absolue ; il est parfois conseillé de procéder systématiquement
à un pivotage. Cette opération est cependant coûteuse en terme d’accès et de transfert de
données en mémoire, surtout si elle nécessite des accès aux disques. Si la matrice assem-
blée est initialement bande, les permutations détruisent généralement cette propriété, ce
qui peut alors devenir très pénalisant.

Stratégies de choix du pivot

Il peut s’agir du pivotage partiel ou du pivotage total :
pivotage partiel
on choisit dans la même colonne k que celle du pivot à changer, la ligne l telle que :

|Al ,k | = max |Am ,k | (1.20)

m =k ,...,n
 12 RÉSOLUTION DES SYSTÈMES LINÉAIRES

pour remplacer le pivot ;

pivotage total
on peut aussi permuter les colonnes ; on cherche alors le nouveau pivot dans toute la
sous-matrice restant à traiter, on choisit donc la ligne l et la colonne c telles que :

|Al ,c | = max |Am ,p | (1.21)

m ,p =k ,...,n

Particularités
Sans recours au pivotage, il y a unicité de la factorisation. Dans le cas où il est nécessaire
de traiter plusieurs second membres, on opère simultanément sur les b. Si les résolutions
doivent être réalisées successivement, et pour chaque autre second membre b′ , on résoud :

LUx′ = b′ (1.22)

– par résolution de Ly′ = b′ appelée descente ;

– puis par résolution de Ux′ = y′ appelée montée.
À chaque fois, il ne s’agit plus que de résoudre un système triangulaire (supérieur ou in-
férieur). On peut d’autre part calculer à faible coût le déterminant de A une fois que l’on
connaît sa décomposition LU :
Y
det A = det(LU) = det L det U = det U = Uk ,k (1.23)
k =1,n

Complexité
n3
Dans le cas d’une matrice stockée pleine, la factorisation LU nécessite ≈ 3
opérations ; une
n2
montée (ou 1 descente) nécessite quant à elle ≈ 2
opérations.

1.3.2 Factorisation de Cholesky

Cette technique a été décrite pour la première fois dans [16] en 1924.

Caractéristiques de la méthode
Même idée que la factorisation LU quand A est symétrique et inversible : A = BT B avec B à
diagonale positive. Dans le cas, le stockage et le calcul ne portent que sur un peu plus de la
moitié de la matrice, par exemple sur sa partie inférieure.

Propriété
Si A est symétrique, définie, positive (donc inversible), il n’est pas utile de pivoter.

Complexité
Pour A pleine, la factorisation BT B requière environ n 3 /6 opérations. Pour A bande, la facto-
risation BT B requière approximativement nl b2 /2 opérations, une montée (ou une descente)
coûte de l’ordre de nl b opérations. Le stockage nécessite aux alentours de nl b réels (géné-
ralement 8 octets). De façon courante, on a l b ≪ n 3 .
Un exemple d’utilisation de cette factorisation pour l’orthogonalisation des vecteurs d’un
sous-espace est présentée dans l’annexe A.
lb 1 1
3. n
= 10
à 100
 1.3 Méthodes directes 13

1.3.3 Factorisation de Crout

C’est la méthode la plus employée. La factori-
sation se met sous la forme A = LDLT avec L boucle sur j = 1, 2, . . . , n
boucle sur i = 1, 2, . . . , j − 1
triangulaire inférieure à diagonale unitaire et D
v i ← A j ,i Ai ,i
diagonale. Cette factorisation s’applique pour fin
v j ← A j ,j − Aj ,1:j −1 v1:j −1
A symétrique, définie positive, mais traite aussi
A j ,j ← v j
les cas où les valeurs propres sont négatives ou Aj +1:n ,j ← v1 (Aj +1:n ,j − Aj +1:n ,1:j −1 v1:j −1 )
j
fin
nulles.

Complexité alg. 1.2 - Factorisation LDLT

stockée pleine
Elle est identique à celle de la factorisation de
Cholesky, mais on n’a plus besoin d’extraire n racines carrées (pour la diagonale).
L’algorithme 1.3.3 présente cette factorisation dans le cas d’une matrice symétrique A
stockée pleine. Il requière l’utilisation d’un vecteur v supplémentaire. Sans pivotage, comme
pour Cholesky et pour Gauss, la factorisation conserve le profil de la matrice écrasée en
mémoire, comme on le verra dans la suite.

1.3.4 Lien avec la condensation de Schur

Super-éléments
Comme c’était le cas pour la factorisation
zone d’intérêt particulier
LU, la factorisation de Cholesky ou Crout est
très liée à la condensation de Schur. Dans le
cas des éléments finis, on peut bâtir ainsi la i
S
notion de super-élément. Considérons une
b
partie d’une structure qui représente une
zone d’intérêt particulier (par exemple, une
fig. 1.3 - Utilisation d’un super-élément
zone qui reste en élasticité alors que le reste
de la structure plastifie, où une zone dans la-
quelle des fonctions de bases particulières ont été ajoutées comme dans la méthode de par-
tition de l’unité. . . ) et qui a été maillée par éléments finis, figure 1.3. L’interface avec le reste
de la structure se fait par le bord de cette zone. Notons avec un indice i les degrés de libertés
des nœuds strictement internes à cette zone et b ceux des nœuds de l’interface, alors :
 
Ki i Ki b
K=  (1.24)
Kb i Kbb

Si on condense sur les nœuds bord on obtient :

S = Kbb − Kb i K−1
i i Ki b (1.25)

c’est la matrice de rigidité du super-élément considéré. Elle est similaire à une matrice de
rigidité élémentaire, pour un élément qui peut être lui-même une petite structure.
 14 RÉSOLUTION DES SYSTÈMES LINÉAIRES

Conservation de la bande
La figure 1.4 montre, à l’itération k , c’est-à-dire lors de la condensation de l’inconnue k sur
les suivantes, les termes de la matrice qui vont être touchés, autrement dit, ceux qui a priori
sont ou deviennent non-nuls. Clairement, la bande est conservée lors de la factorisation. On
peut de la même façon montrer que seuls les termes sous le profil sont modifiés, et donc que
la factorisation conserve le profil de la matrice.

1.3.5 Autres techniques

L’adaptation de la factorisation de Crout au cas où A est
non-symétrique est classique, c’est la factorisation LDMT . k

D’autres variantes des techniques précédentes per-

mettent de tirer mieux parti des conditions de stockage, en
particulier pour optimiser les échanges de données avec
k
la partie des informations stockées sur disque, comme la
coefficients
méthode frontale [17]. modifiés à
l’étape k
D’autres versions encore permettent de réaliser des
factorisations en parallèle de façon assez performante,
comme les méthodes multi-frontales [18]. fig. 1.4 - Étape k de la
factorisation

1.4 Quelques particularités en éléments finis

1.4.1 Conditions aux limites en déplacement
La prise en compte de déplacements imposés non-nuls a une influence sur l’utilisation des
solveurs dans le sens où la matrice de rigidité est modifiée par ces contraintes.

Séparation des degrés de liberté

Le problème correspondant de minimisation du potentiel sous contrainte est le suivant :

1
min J1 (u) = uT Ku − fT u (1.26)
Cu=ud 2
C est une matrice booléenne dont le but est d’extraire une partie des degrés de libertés, ceux
sur lesquels la condition d’égalité aux déplacements imposés ud est imposée.
Reprenons l’écriture par blocs déjà présentés dans l’équation (1.15) en séparant les de-
grés de liberté imposés u2 = Cu = ud . Le problème (1.26) est ainsi équivalent à :

1
min uT1 K11 u1 − fT1 u1 + uT1 K12 ud + constante (1.27)
u1 2

dont l’équation d’Euler correspondante est K11 u1 = f1 − K12 ud . On a ainsi un problème de

taille réduite où les charges extérieures sont modifiées par le terme −K12 ud 4 . Cette méthode
est effectivement employée dans des codes comme abaqusTM. L’inconvénient de cette mé-
thode réside dans le tri explicite nécessaire dans les degrés de liberté qui rempli donc le
4. ce sont les réactions à l’appui qui permettent d’avoir un déplacement bien égal à ud
1.4 Quelques particularités en éléments finis 15

terme K12 , et qui entraîne un surcoût important lorsque le nombre de degrés de liberté blo-
qués est grand.

Multiplicateurs de Lagrange
Une autre technique consiste à reformuler le problème (1.26) en relaxant la contrainte à
l’aide de multiplicateur de Lagrange l et donc de transformer le problème de minimum en
celui de la recherche du point selle du Lagrangien :

extr J1 (u) − lT (Cu − ud ) (1.28)

u,l

Il n’est plus nécessaire de séparer les degrés de libertés, mais la taille du problème a aug-
menté ; en effet les équations d’Euler sont cette fois :
  ! !
K −CT u f
  = (1.29)
−C 0 l −ud
| {z }
A

On peut remarquer que la première équation montre que les multiplicateurs de Lagrange
sont bien les réactions à l’appui.
L’inconvénient vient du fait que la taille du problème augmente, et du fait que A est une
matrice symétrique mais plus définie positive ; elle possède des valeurs propres négatives.
On a cependant le résultat suivant si K est symétrique, définie positive :

A inversible ⇔ C injective ⇔ kerr C = ∅ (1.30)

c’est-à-dire que l’on n’écrit que des liaisons indépendantes. Une factorisation de Gauss risque
de détruire la structure bande de la matrice K et la factorisation de Crout va demander des
pivotages.
Une modification de cette technique permet, avec une numérotation ad hoc qui conserve
la bande, de ne pas avoir besoin de pivoter. Il s’agit d’une technique de double multiplica-
teurs utilisée dans le code Cast3MTM . Elle peut être illustrée de la manière suivante :
 T   
l !T !
1 1   −1 1 1  l1  u l
     d 1
extr J1 (u) + l2   1 −1 1  l2  − (1.31)
u,l1 ,l2 2     ud l2
u 1 1 0 u2
| {z }
1
− (l 1 − l 2 )2 + (l 1 + l 2 )T (Cu − u d )
2

dont les équations d’Euler sont :

Ku = f − CT (l 1 + l 2 )
l1 = l2 (1.32)

Cu = ud

Bien entendu, en cas de nombreux blocages, la taille du problème croît en conséquence.

 16 RÉSOLUTION DES SYSTÈMES LINÉAIRES

Valeurs très grandes

On peut aussi considérer une condition aux limites en déplacement imposé comme une
modélisation de l’action de l’extérieur. Celui-ci pourrait aussi réagir comme un ressort ayant
une grande raideur k , avec un déplacement imposé à la base. Le potentiel élastique de l’en-
semble structure + ressort extérieur est alors :
1
J1 (u) + (Cu − ud )T k (Cu − ud ) (1.33)
2
Les nouvelles équations d’Euler de cette formulation sont alors :

(K + CT k C)u = f + CT k Cud (1.34)

On ne fait donc que rajouter des termes dans la matrice de raideur sur les degrés de libertés
incriminés et dans le vecteur des forces généralisées sur les efforts duaux.
Le problème est le choix de k : pour des valeurs trop petites, la condition est mal impo-
sée, pour des valeurs trop grandes, le système devient extrêmement mal conditionné, voire
numériquement singulier.

1.4.2 Relations linéaires entre degrés de liberté

Cette situation se présente par exemple fans le cas d’utilisation de symétrie par rapport à
un plan non parallèle à un des plan des axes globaux de l’analyse, dans le cas de conditions
aux limite périodiques sur une cellule de base. . . À chaque fois, il s’agit de contraintes de
la forme Cu = ud où C n’est cette fois-ci pas forcément booléenne. Les techniques précé-
dentes s’appliquent bien entendu sur ces cas de figure tout aussi bien. . . et avec les mêmes
limitations.

1.4.3 Système singulier : semi défini positif

Exemple
Les exemples sont nombreux : il s’agit de structure avec des modes rigides, de problèmes de
sous-intégration. . . il y a alors des déplacements non nuls à énergie de déformation nulle,
autrement dit des modes rigides :

1
uT Ku = 0 (1.35)
2
Dans la suite, on s’intéressera au cas où cas est A symétrique, sachant que la généralisation
est tout à fait possible.

Noyau
L’ensemble des solutions indépendantes, stockées dans les colonne de R, telles que AR = 0,
forme le noyau de A. S’il y en a p , le rang de A est n − p .

Existence de solutions
Pour qu’il existe des solutions, une condition nécessaire et suffisante est que b soit orthogo-
nal au noyau de A, c’est-à-dire RT b = 0. Exemple : les efforts ne sollicitent pas les modes de
solide rigide.
 1.5 Méthodes itératives stationnaires 17

– Si A est non symétrique, la condition est que b appartienne au noyau adjoint de A.

– Si x0 est une solution de ce système, l’ensemble des solutions est de la forme x0 + Rα
où α est quelconque.

Inverses généralisées (ou pseudo-inverses)

A+ telle que AA+ A = A. Si A est inversible, il n’y a qu’une pseudo-inverse A+ = A−1 . Pour
A symétrique, toutes les pseudo-inverses sont de la forme A+ + BT RT + RB, avec B quel-
conque [19].

Factorisation et construction d’une inverse généralisée

Supposons que le système puisse s’écrire sous la forme de blocs à l’image de (1.15), où A11
est la plus grande sous-matrice régulière. On a alors S = A22 − A21 A−1
11 A12 = 0. Une pseudo-
inverse de A est :
 
A−1 0
 11  (1.36)
0 0

et la solution est :
! !  
x1 A−1
11 b1 −A12
= +  x2 (1.37)
x2 0 Id
| {z } | {z } |{z}
x0 R α

La factorisation fait progresser l’élimination tant qu’un pivot non nul existe. Dans notre cas,
avec pivotage, on tombe sur une sous-matrice nulle. Pour les matrices semi-définies, il n’est
pas nécessaire de pivoter car un pivot nul apparaît lorsque la ligne du pivot est une com-
binaison linéaire des précédentes. L’inconnue correspondante est repérée et fixée à 0, voir
x0 dans (1.37), c’est-à-dire qu’on ajoute des liaisons complémentaires pour fixer les modes
de solide rigide — à condition que le terme correspondant au second membre condensé
soit nul aussi. La colonne correspondante de la matrice stocke le vecteur associé, voir R
dans (1.37).

1.5 Méthodes itératives stationnaires

De manière générale avec les techniques itératives, il faut se poser des questions telles que :
– quelles sont les conditions d’emploi de ces techniques afin d’assurer la convergence
vers la solution cherchée ?
– quel critère d’arrêt des itérations employer ?
– quels sont les coûts de ces méthodes ?
Dans le cas des méthodes itératives stationnaires, on rencontre les particularités suivantes :
– on accède à A uniquement par le produit Aw (phase la plus coûteuse). Éventuellement,
comme en éléments finis, on peut le faire sous forme élémentaire :
X
Aw = Ae we (1.38)
e
 18 RÉSOLUTION DES SYSTÈMES LINÉAIRES

– on construit une suite de solutions approchées :

x(k ) → A−1 b (1.39)

k →∞

à partir d’une récurrence linéaire x(k ) = Bx(k −1) + c où B et c sont indépendantes de k 5 .

Ce sont des méthodes de type point fixe : si on converge, c’est vers :

x = Bx + c ⇒ (x(k ) − x) = B(x(k −1) − x) (1.40)

Pour avoir kx(k ) − xk →k →∞ 0, il faut et il suffit que ρ(B) < 1 où ρ(B) est le rayon spectral de B :
c’est le sup des valeurs absolues des valeurs propres de B. Le taux de convergence asympto-
tique est défini par :

τ = − ln ρ(B) (1.41)

Il peut servir à estimer le nombre d’itérations nécessaire pour atteindre le niveau de conver-
gence souhaité ; en effet on a :

kx(k ) − xk kx(k ) − xk
(k ) k (0) k
(x − x) = B (x − x) ⇒ ¶ kB k ⇒ − ln ¶ kτ (1.42)
kx(0) − xk kx(0) − xk

L’intérêt est de pouvoir estimer le nombre d’itérations nécessaires pour diviser le niveau
d’erreur par un coefficient donné (utilisation d’un critère d’arrêt), connaissant le taux de
convergence asymptotique de la méthode. Pour plus d’informations sur les méthodes itéra-
tives stationnaires, le lecteur pourra consulter [20].

1.5.1 Méthode de Jacobi

Pour les systèmes linéaires, cette méthode, d’algo-
rithme 1.5, est décrite pour mémoire seulement, Vecteur de départ x(0)
car elle possède un faible taux de convergence et tant que R > ǫ
itération sur k
dans nos applications, est remplacée par d’autre boucle sur i = 1, n
(k −1) (k −1)
x̄i ← Ai ,1:i −1 x1:i −1 + Ai ,i +1:n xi +1:n
techniques plus efficaces. Elle illustre cependant
x̄i ← (bi − x̄i )/Ai ,i
bien les techniques de point fixe : si on suppose fin
x(k ) ← x̄
connues toutes les inconnues sauf l’inconnue i , fin
cette dernière peut être déterminée à l’aide de
l’équation i . Procédant ainsi sur toutes les incon- alg. 1.5 - Méthode de Jacobi
nues, on obtient un nouvel itéré xk à partir de l’an-
cien xk −1 . Cependant, comme la méthode de Gauss-Seidel, elle peut présenter un intérêt
dans le cas de systèmes non-linéaires. Sous forme matricielle, cela se traduit par :

xk = (I − D−1 A)xk −1 + D−1 b (1.43)

où D est la matrice diagonale, contenant la diagonale de A. La matrice d’itération est donc

ici B = I − D−1 A. On peut noter l’utilisation d’un vecteur de stockage x̄ supplémentaire.
5. Cette caractéristique n’existe pas pour les méthodes non stationnaires.
 1.5 Méthodes itératives stationnaires 19

Convergence
Si A est symétrique, définie positive, alors ρJ = ρ(B) < 1, mais en éléments finis, on a une
évaluation de ce rayon spectral qui varie comme :

h2
ρJ = 1 − + O (h 4 ) (1.44)
2
d’où un taux de convergence faible lorsque le maillage est raffiné (donc lorsque h est faible) :

h2
τ≈ (1.45)
2

1.5.2 Méthode de Gauss-Seidel

L’idée est similaire à celle de la méthode de Jacobi, à ceci
Vecteur de départ x(0) près que l’on utilise les nouvelles valeurs estimées dès
tant que R > ǫ
qu’elles sont disponibles, et non pas à la fin de l’itération.
itération sur k
boucle sur i = 1, 2, . . . , n Sous forme matricielle, cela se traduit par :
(k −1)
xi ←0
(k −1)
β ← Ai ,1:n x1:n xk = (D − E)−1 (F xk −1 + b) (1.46)
(k )
xi ← (bi − β)/Ai ,i
fin où −E est la partie triangulaire inférieure de A, et −F sa
fin
partie triangulaire supérieure. La matrice d’itération est
alg. 1.6 - Méthode de ici :
Gauss-Seidel
B = (D − E)−1 F = (I − D−1 E)−1 D−1 F (1.47)

Quand A est symétrique, définie positive, on est assuré de la convergence puisque ρGS < 1.
Mais l’estimation de la valeur du rayon spectral est difficile cas il dépend en particulier de la
numérotation des degrés de liberté employée.

1.5.3 Méthode de relaxation

Cette technique, SOR pour successive over relaxation, est attribuée à Southwell [21] et date de
1946. Elle utilise comme base la méthode de Gauss-Seidel ainsi qu’une technique mainte-
nant souvent utilisée avec d’autre approches : la tech-
Vecteur de départ x(0) nique de relaxation.
tant que R > ǫ
Le principe, expliqué par l’algorithme 1.7, consiste
itération sur k
boucle sur i = 1, 2, . . . , n à modifier (par relaxation) l’itéré xk fourni par l’algo-
(k −1)
α ← xi
(k −1) rithme de Gauss-Seidel à chaque itération. Pour cela, il
xi ←0
(k −1)
β ← Ai ,1:n x1:n faut conserver l’itéré précédent xk −1 et modifier l’itéré
(k )
xi ← (1 − ω)α + ωβ
fin
courant comme suit :
fin
xk ← ωxk + (1 − ω)xk −1 (1.48)
alg. 1.7 - Méthode SOR
où ω est appelé le paramètre de relaxation. Si ω < 1, on
parle de sous-relaxation, et si ω > 1, de sur-relaxation. On peut montrer que la matrice d’ité-
ration correspondante est :

B = (I − ωD−1 E)−1 [(1 − ω)I + ωD−1 F] (1.49)

 20 RÉSOLUTION DES SYSTÈMES LINÉAIRES

Concernant la convergence, le théorème de Kahan donne une condition nécessaire : 0 <

ω < 2. Si ω = 0, l’algorithme ne progresse plus. Le théorème d’Ostrowski-Reich permet de
conclure à la convergence ∀ω ∈]0, 2[ quand A est symétrique, définie positive.
Le choix de ω est fait à partir d’heuristiques,
comme ω = 2 − O (h) pour une discrétisation Vecteur de départ x(0)
tant que R > ǫ
éléments finis dont la taille de maille est h. Il
itération sur k
existe aussi des techniques d’optimisation de ω boucle sur i = 1, 2, . . . , n
(k −1) (k −1)
x̄i ← Ai ,1:i −1 x1:i −1 + Ai ,i +1:n xi +1:n
en cours de route, comme la technique d’accélé-
x̄i ← (bi − x̄i )/Ai ,i
ration d’Aitken, voir annexe C. fin
x(k ) ← x(k −1) + ω(x̄ − x(k −1) )
On peut aussi s’inspirer de résultats dans des fin
cas particuliers. Par exemple dans un cas tridia-
gonal bloc, on montre que la valeur optimale de alg. 1.8 - Méthode de Jacobi relaxée
ω est :
2
ωopt = p (1.50)
1 + 1 − ρJ2
La valeur de ρJ est ainsi à estimer sur le problème en question. On a alors ρSOR = ωopt − 1, et,
dans le cas éléments finis, τ ≈ 2h (un ordre est gagné par rapport à la méthode de Jacobi).
La technique de relaxation peut être appliquée à d’autres algorithmes. Par exemple, un
algorithme de Jacobi relaxé s’écrit sous la forme donnée dans l’algorithme 1.8.

1.5.4 Méthode d’Uzawa

Dans certains cas de figure comme les problèmes de point selle ou les liaisons traitées par
multiplicateur de Lagrange par exemple, on peut être amené à résoudre des problèmes de
la forme :
  ! !
K C x f
  = (1.51)
CT 0 l g
Pour que le problème continue d’être inver-
sible, il suffit que C soit injective. Vecteurs de départ x(0) , l(0)
tant que R > ǫ, itérer sur k
La méthode d’Uzawa développée en 1958
résoudre Kx(k ) = b − Cl(k −1)
consiste à résoudre itérativement de façon dé- réactualiser l(k ) ← l(k −1) +α(CT x(k ) −c)
fin
couplée, comme le montre l’algorithme 1.9 où α
est un paramètre de la méthode. Pour analyser alg. 1.9 - Méthode d’Uzawa
cet algorithme, on peut condenser la variable x(k )
sur la variable l(k ) pour trouver :

l(k ) = (I − αCT K−1 C) l(k −1) + α(CT K−1 f − g) (1.52)

| {z }
B

Si on appelle λ1 , . . . , λp , les valeurs propres de la matrice CT K−1 C, on en déduit l’équivalence

suivante :

̺(B) < 1 ⇔ ∀i , |1 − αλi | < 1 (1.53)

 1.6 Méthodes itératives non-stationnaires 21

Avec CT K−1 C symétrique, définie positive, λi > 0 donc il faut prendre :

2
0<α< (1.54)
maxi λi

1.6 Méthodes itératives non-stationnaires

Ces méthodes se décrivent de la même façon que les précédentes, à ceci près que les opéra-
teurs B et c dépendent des itérations.

1.6.1 Méthode générique de projection

Reformulons le problème de départ en définissant le résidu :

r = b − Ax (1.55)

Il s’agit maintenant d’avoir r = 0. L’idée est ici de se trouver

avec un problème de plus petite taille, de structure plus b − Ax

simple (par exemple, tridiagonal) par plusieurs étapes de

projection. x
U
Considérons un sous-espace U de l’espace dans lequel
la solution est cherchée. Le problème approché est défini fig. 1.10 - Réduction dans
un sous-espace
par x = Qx̂ ∈ U avec b − Ax ⊥ U , comme décrit sur la fi-
gure 1.10. Il est clair que si U est l’espace complet, ce problème est le problème de départ.
Dans tous les cas, il conduit à :

QT (b − AQx̂) = 0 ⇒ (QT AQ) x̂ = QT b (1.56)

| {z } |{z}
 b̂

Une méthode de projection peut donc se schématiser de la façon suivante :

1. r ← b − Ax

2. x̂ = Â−1 QT r

3. x ← x + Qx̂

4. itération suivante

Le choix de Q à chaque itération dépend bien entendu de la méthode employée.

1.6.2 Méthode du gradient

Pour A symétrique, définie positive, la solution approchée est modifiée de la manière sui-
vante :

x(k +1) ← x(k ) + µ(k ) d(k ) (1.57)

où d(k ) est la direction de recherche dans la direction du gradient :

d(k ) = Ax(k ) − b (1.58)

 22 RÉSOLUTION DES SYSTÈMES LINÉAIRES

µ(k ) est le paramètre d’optimalité de la nouvelle solution dans la direction d(k ) . Il se calcule
comme :
T T
d(k ) (b − Ax(k ) ) d(k ) d(k )
(k )
µ = T
=− T
(1.59)
d(k ) Ad(k ) d(k ) Ad(k )
Il s’agit d’une méthode de projection pour laquelle la taille du sous-espace est 1 car on choi-
sit pour Q le sous-espace généré par le seul vecteur résidu Q = r.
La convergence de cette méthode est assez lente, d’où l’intervention d’une technique de
conjugaison pour la méthode suivante.

1.6.3 Méthode du gradient conjugué

Cette méthode est généralement attribuée à Hestenes et Stiefel [22] en 1952. Elle est utilisée
dans le cas où A est symétrique, définie positive.

Description de la méthode
L’initialisation consiste à choisir un vecteur de départ x(0) , à calculer le résidu associé r(0) =
b − Ax(0) et de prendre comme première direction de recherche d(0) = r(0) .
Les itérations sont ensuite de la forme :
1. itéré x(k +1) = x(k ) + µ(k ) d(k ) avec le pas optimal
T T
r(k ) d(k ) r(k ) r(k )
(k )
µ = T
= T
(1.60)
d(k ) Ad(k ) d(k ) Ad(k )
2. résidu r(k +1) = b − Ax(k +1)
3. direction de recherche d(k +1) = r(k +1) + λ(k ) d(k ) qu’on orthogonalise par rapport à la
précédente :
T T
r(k +1) Ad(k ) r(k +1) r(k +1)
(k +1) T (k ) (k )
d d =0 ⇒ λ =− T
= T
(1.61)
d(k ) Ad(k ) r(k ) r(k )
4. itération suivante

Propriétés
T T T T
Si i < j , alors r(i ) r(j ) = 0 = d(i ) r(j ) = d(i ) Ad(j ) = r(i ) Ad(j ) . On construit donc des vecteurs
orthogonaux au sens de A, d(0) , d(1) , . . . qui forment donc une base. On minimise à chaque
fois sur un sous-espace de taille de plus en plus grande.
En théorie, cette méthode est donc une méthode directe en n étapes. En pratique, on a
progressivement perte d’orthogonalité due aux arrondis, c’est donc en ce sens une méthode
itérative.
On a alors la propriété de convergence suivante :
p (k )
 cond A − 1
kx(k ) − xkA ¶ 2  p (0)
 kx − xkA (1.62)
cond A + 1

d’où, en éléments finis :

 
1
τ=O p = O (h) (1.63)
cond A
 1.7 Critères d’arrêt 23

1.6.4 Préconditionnement
Si on peut avoir facilement une approximation
Vecteur de départ x(0) spectrale M de A, l’idée est de chercher à ré-
Résidu d’équilibre r(0) ← b − Ax(0)
soudre :
tant que R > ǫ, itérer sur k
résoudre Mz(k −1) = r(k −1)
T M−1 Ax = M−1 b (1.64)
ρk −1 ← r(k −1) z(k −1)
si k = 1 alors
M symétrique, définie positive. Ce problème de-
d(1) ← z(0)
sinon vient non symétrique, mais on peut lui substi-
ρ
βk −1 ← ρk −1
k −2 tuer :
d(k ) ← z(k −1) + βk −1 d(k −1)
fin
€ Š
M−1/2 AM−1/2 M1/2 x = M−1/2 b (1.65)
q(k ) ← Ad(k )
ρ
αk ← (k )kT−1(k )
d q Le choix du préconditionneur dépend du pro-
x(k ) ← x(k −1) + αk d(k )
r(k ) ← b(k ) − Ax(k ) = r(k −1) − αk q(k ) blème à traiter. Dans le meilleur des cas, il est fa-
fin cile à inverser et représente bien la matrice A. En
fait, si M = I il est facile à inverser, mais est très
alg. 1.11 - Gradient conjugué
préconditionné différent de A, et si M = A il y a convergence en
une itération, mais dans laquelle il faut résoudre
le problème de départ !
La convergence est cette fois-ci en :
 
 1 
τ = O p  (1.66)
cond M−1 A

Comme exemples de préconditionneurs classiques, on a la factorisation de Cholesky incom-

plète, SSOR, etc. L’algorithme 1.11 présente cette méthode. Le choix et le développement de
préconditionneurs adaptés font encore aujourd’hui l’objet de recherches.

1.6.5 Autres méthodes itératives

Pour information, on peut citer :
– les méthodes dites « rapides » dédiées à des applications particulières : Fourier (O (n log n)),
multigrilles (O (n)) ;
– pour A symétrique et définie : Chebyshev [23] ;
– pour A symétrique mais non définie : MinRES (résidu minimum, minimal residual),
SymmLQ (symmetric LQ) [24] ;
– pour A non symétrique connaissant AT : QMR (quasi minimal residual) [25] ;
– pour A non symétrique : CGS (conjugate gradient square) [26], BiCGStab (bi-conjugate
gradient stabilized) [27], GMRES (generalized minimal residual) [28], GMRES restar-
ted. . .

1.7 Critères d’arrêt

Le critère d’arrêt recherché est celui de l’écart en solution tel que e(k ) = x(k ) − x soit suffi-
samment petit. Malheureusement, x est la solution à convergence et on ne peut atteindre
24 RÉSOLUTION DES SYSTÈMES LINÉAIRES

l’erreur en solution au cours des itérations. La quantité que l’on peut espérer atteindre est
liée au résidu :

r(k ) = b − Ax(k ) = −Ae(k ) (1.67)

Un critère souvent employé est le suivant :

kr(k ) k
<ǫ (1.68)
kAx(k ) k + kbk

où ǫ est un seuil fixé a priori, que l’on peut quantifier par rapport à la précision souhaitée et
la précision machine. Parfois, on remplace au dénominateur kAx(k ) k par kbk. Le problème
vient du fait, comme on a déjà pu le voir dans un exemple, que l’erreur en solution et le
résidu sont liés de la façon suivante :

ke(k ) k ¶ kA−1 k (1.69)

donc :

r(k ) e(k )
<ǫ ⇒ < ǫ cond A (1.70)
b x
et, en tout état de cause, le paramètre ǫ devrait être choisi aussi fonction du conditionne-
ment de A, qu’il faut donc estimer pour chaque problème traité.
Des critères d’arrêt en cas de non-convergence sont aussi mis en place. Par exemple,
kr(k ) k décroît trop lentement ou ne décroît plus (mais il faut faire attention aux plateaux
possibles avec certaines méthodes), ou la résolution est trop longue (critère en nombre maxi
d’itérations). Il s’agit alors de sortie en cas d’erreur.
 2
SYSTÈMES AUX
VALEURS PROPRES

On considère ici les problèmes, dits aux valeurs propres, de la forme : trouver les caractéris-
tiques propres, c’est-à-dire les couples (λ, x) où x 6= 0 tels que Ax = λx et l’on se restreindra
aux cas où A est symétrique.
Un problème aux valeurs propres généralisé est de la forme Kx = ω2 Mx, avec x 6= 0, K
symétrique, définie positive et M symétrique positive. Dans ce cas, on peut utiliser la factori-
sation de Cholesky de M, soit M = LLT et avec y = LT x, on est ramené au problème précédent
L−1 KL−T y = ω2 y.
En général, dans les problèmes qui nous concernent, on ne veut qu’environ 30 valeurs
et/ou vecteurs propres, de préférence à basse fréquence (dynamique lente, vibrations [29,
30, 31], flambage. . . ). Dans le cas où il y a de l’amortissement, le système à résoudre est :

Mẍ + Cẋ + Kx = f(t ) (2.1)

La matrice d’amortissement, C, est souvent mal maîtrisée et si l’amortissement est peu in-
fluent à basse fréquence, ce qui n’est pas évident pour tous les cas, on est ramené à :

Mẍ + Kx = f(t ) (2.2)

Avec une discrétisation éléments finis standard, les hautes fréquences obtenues numérique-
ment sont éloignées des vraies.
Enfin, on prévoit parfois d’avoir des problèmes si les valeurs propres sont mal séparées,
proches ou même confondues en cas de symétries, par exemple. Une référence de base
concernant ces problèmes est [32, 33].

2.1 Préliminaires
2.1.1 Rappels
Matrice quelconque
λ est valeur propre de A ssi det(A − λI) = 0. Dans C, A admet n valeurs propres distinctes ou
non, réelles ou complexes. A a n valeurs propres distinctes ssi A est diagonalisable.
 26 SYSTÈMES AUX VALEURS PROPRES

Matrice réelle
A est symétrique ssi A est diagonalisable dans R. A est symétrique, définie positive ssi :
– toutes les valeurs propres sont strictement positives ;
– ∃B symétrique, définie positive telle que A = B2 .
Dans C, A admet n valeurs propres distinctes ou non, réelles ou complexes. A a n valeurs
propres distinctes implique que A est diagonalisable.

Quotient de Rayleigh pour une matrice réelle, symétrique

Le quotient de Rayleigh est défini comme étant, pour x 6= 0 :

xT Ax xT Kx
R(x) = ou R(x) = (2.3)
xT x xT Mx

Il vient alors l’équivalence suivante :

(
R est extremum en x 6= 0
(λ, x) est caractéristique propre de A ⇔ (2.4)
R(x) = λ

La i e valeur propre λi est le minimum de R(x), x 6= 0 sur le sous espace orthogonal aux précé-
dents vecteurs propres. De plus, en dimension finie, la valeur propre maxi est λn = maxx6=0 R(x)
dans l’espace complet de dimension n.

2.1.2 Domaine d’emploi des différentes méthodes

En comparaison de la résolution d’un système linéaire, la résolution d’un système aux va-
leurs propres est beaucoup plus coûteuse. En fait, il s’agit d’un problème non-linéaire, et
beaucoup de méthodes s’appuient sur la résolution d’une succession de systèmes linéaires,
en interne. Elles souffriront alors des problèmes affectant ces solveurs, en particulier le mau-
vais conditionnement. L’ordre de grandeur de la taille des systèmes aux valeurs propres que
l’on va considérer dans cette partie est de l’ordre de n = 250 à n = 25 000.
La première méthode est celle de l’équation caractéristique. On peut en effet recher-
cher les zéros du polynôme caractéristique de la matrice, mais si les coefficients de ce po-
lynôme peuvent être assez facilement atteints, une petite perturbation sur ceux-ci conduit
à d’énormes variations sur ses racines. Cette méthode est donc réservée à des tailles de sys-
tème réduites, jusqu’à n ≈ 10. On peut montrer mathématiquement que, dès que n ¾ 5, la
méthode employée ne peut être qu’itérative : le problème ne peut plus être résolu numéri-
quement en un nombre fini d’opérations.
Les méthodes de Jacobi, des puissances et de Householder vont ensuite pouvoir être em-
ployées pour des systèmes dont la taille va jusqu’à n ≈ 250.
Au delà, elle doivent être adaptées à un stockage bande, et le remplissage des matrices
au cours des itérations doit être étudié. Avec les méthodes de réduction, ces techniques per-
mettent d’atteindre des systèmes dont la taille va jusqu’à n ≈ 2 500.
Jusqu’à n ≈ 25 000, les méthodes des puissances, plus ou moins couplées aux méthodes
de sous-espaces, la méthode de Lanczos avec stratégie de réorthogonalisation permettent
encore de traiter le problème.
 2.2 Méthode de Jacobi 27

Pour dépasser ces tailles de problème, il est ensuite nécessaire de faire appel à des tech-
niques particulières. Par exemple, les techniques de parallélisation, vectorisation, optimisa-
tion des entrées-sorties, les méthodes de condensation et sous-structuration (Guyan-Irons,
Craig-Bampton. . . ) Ces dernières méthodes ne seront pas étudiées ici.

2.2 Méthode de Jacobi

Cette méthode, datant de 1846 [34], est décrite ici, dans le cas particulier où A est symé-
trique, à des fins de rappel uniquement puisqu’elle est peu efficace pour les problèmes qui
nous concernent.
Dans cette méthode, on cherche à ramener la
tant que R > ǫ, itérer sur k matrice A à une forme diagonale pour trouver tous
choix du terme à annuler Ap,q
(k −1) (k −1) ses valeurs propres et vecteurs propres, et ce par
Aq,q −Ap ,p
ν← (k −1)
2Ap ,q une séquence de transformations orthogonales :
calcul de la tangente t p
calcul du cosinus c ← 1/ 1 + t 2 T
A(k +1) = R(k ) A(k ) R(k ) (2.5)
calcul du sinus s ← c t
A(k ) ← A(k −1)
T
(k ) (k −1) (k −1)
Ap,1:n ← c Ap,1:n − s Aq,1:n avec les rotations R(k ) vérifiant R(k ) R(k ) = I, trans-
(k ) (k −1) (k −1)
Aq,1:n ← s Ap,1:n + c Aq,1:n formations qui sont des changements de base
(k ) (k −1) (k −1)
Ap,p ← Ap,p − t Ap,q
(k )
Aq,q
(k −1) (k −1)
← Aq,q + t Ap,q
de type isométries ne modifiant pas les valeurs
(k )
Ap,q ←0 propres.
(k )
Aq,p ←0
On commence par chercher l’élément hors dia-
fin
gonale de module maximal :
alg. 2.1 - Méthode de Jacobi
(k )
|A(k )
p,q | = max |Ai ,j | (2.6)
i 6= j

On construit alors la matrice de rotation d’un angle θ ad hoc dans le plan (ep , eq ) de façon à
(k +1) (k +1)
annuler dans l’itéré suivant les termes Aq,p et Ap,q :

R(k ) = 1 + (cos θ − 1)(ep eTp + eq eqT ) + (sin θ − 1)ep eTp − (sin θ + 1)eq eqT (2.7)

avec :
1 h i 1
A(k +1)
p,q = sin 2θ A(k
p,p
)
− A(k ) (k )
q,q + cos 2θAp,q = 0 (2.8)
2 2
ce qui entraîne :
(k ) (k )
Ap,p − Aq,q
cotan 2θ = (k )
(2.9)
Ap,q

Seuls changent les termes sur les lignes ou colonnes p ou q . Dans l’algorithme 2.1, la tan-
gente t est calculée comme suit :

 1 si ν = 0
t= p (2.10)
 −ν + sgn(ν) 1 + ν2 sinon
 28 SYSTÈMES AUX VALEURS PROPRES

Au cours des itérations, un terme nul peut redevenir non nul mais :
X
(k ) 2
a (k ) = Ai ,j → =0 (2.11)
k →∞
i 6= j

En effet :
X h i 2 X h i2 h i2 ‹ X h i2 h i2 ‹
(k +1) (k +1) (k +1) (k +1) (k +1)
a (k +1) = Ai ,j + Ai ,p + Ai ,q + Ap,j + Aq,j
i 6= j ,p | {z } i 6=p,q | {z } j 6=p,q | {z }
q,j 6=p,q
h i2 h i2 h i2 h i2 h i2
(k )
Ai ,j (k ) (k )
Ai ,p + Ai ,q
(k )
Ap,j + Aq,j
(k ) (2.12)
h i2
= a (k ) − 2 A(k
p,q
)

donc a (k +1) < a (k ) .

Cet algorithme est stable et simple, mais la convergence est lente et il s’avère coûteux en
accès mémoire si n est grand. On a les estimations suivantes :
 
2 2
a (k +1) ¶ 1 − a (k ) et τ ≈ 2 (2.13)
n(n − 1) n

2.3 Méthode des puissances

On se place ici dans le cas où A est symétrique, définie positive et dont les valeurs propres
sont toutes simples, rangées par ordre de module décroissant λi .
On construit un itéré de départ x(0) puis les itérés
x (0) , kx (0) k = 1 successifs par récurrence de la manière suivante :
tant que R > ǫ, itérer sur k
v ← Kx (k −1) x(k +1) = Ax(k ) (2.14)
résoudre Mx(k ) = v
α ← kx(k ) k
λ(k ) ← R(x(k ) ) =
T
x(k ) Kx(k ) Si (vi )i est la base propre, on peut écrire les itérés
T
x(k ) Mx(k )
normalisation x(k ) ← xα
(k )
dans cette base. En particulier :
fin X
x(0) = α i vi (2.15)
alg. 2.2 - Méthode des puissances i

alors :
  k 
X n
X λi
x (k )
= λki αi vi = λk1 α1 v1 + αi vi  (2.16)
i i =2
λ1

|λi |
et si α1 6= 0 puisque |λ1 |
< 1 alors :

x(k ) → λk1 α1 v1 (2.17)

k →∞

En pratique, il est utile de normer x(k ) à chaque itération pour éviter l’explosion et dans ce
cas obtenir x(k ) → v 1 .
P
Pour comprendre le comportement de la convergence, prenons x(0) = i vi , alors :
n 
X 
x(k ) λ k i
= v1 + vi (2.18)
λk1 i =2
λ1
2.3 Méthode des puissances 29

donc :
 
 (k +1)  
 x  Xn  2(k +1) 1/2
  λ
 λk +1 − v1  
i

 1   λ 1

λ2
 
  =  i =2n    → (2.19)
 (k )   X λi 2k  k →∞ λ1
 x   
  
 λk − v1  i =2
λ 1
 1 

Le taux de convergence est alors estimé comme suit :

λ2 λ1
τ ≈ − ln = ln (2.20)
λ1 λ2

Bien entendu, on peut avoir des ennuis si au départ x(0) est orthogonal à v1 . Par rapport à la
méthode de Jacobi, la convergence est maintenant indépendante de la taille du problème,
n.

Cas de valeurs propres multiples ou proches

Dans le cas de valeurs propres multiples λ1 , la convergence se fait vers λ1 et vers une com-
binaison linéaire des vecteurs propres associés.
δλ
Dans le cas de valeurs propres proches, λ2 = λ1 + δλ, on a τ ≈ λ1
et la convergence est
donc très lente.

Contrôle de la convergence
Comme critère de convergence sur les valeurs propres, on peut utiliser le fait que kx(k +1) k/kx(k ) k
se stabilise, et pour la convergence sur les vecteurs propres, kx(k +1) − x(k ) k < ǫ ou :

kAx(k +1) − λ(k +1) x(k +1) k

< ǫ′ (2.21)
kAx(k +1) + λ(k +1) x(k +1) k

On conseille aussi d’utiliser la norme infinie pour ces tests.

Obtention des modes suivants

La technique utilisée est la déflation orthogonale : on
v1 prend x(0) orthogonal à v1 . Pour éviter les pertes d’orthogo-
x
nalité par cumul d’arrondis, on réorthogonalise à chaque
itération. Pour cela, on utilise un projecteur comme décrit
sur la figure 2.3 :
fig. 2.3 - Projection ortho-
v1 vT1
gonale à v1 P1 = I − (2.22)
vT1 v1

d’où :
vT1 x
P1 x = x − v1 (2.23)
vT1 v1

On a alors convergence vers (λ2 , v2 ). On peut procéder de façon similaire pour accéder à la
troisième caractéristique propre, et ainsi de suite.
 30 SYSTÈMES AUX VALEURS PROPRES

2.4 Méthode des itérations inverses

Comme décrit dans l’algorithme 2.4, on choisit en-
x(0) , kx(0) k = 1 core un x(0) arbitraire : les solutions classiques sont
tant que R > ǫ, itérer sur k
de prendre un vecteur de composantes aléatoires,
v ← Mx(k −1)
résoudre Kx(k ) = v normé, quitte à essayer avec plusieurs vecteurs de
α ← kx(k ) k
T
x(k ) Kx(k )
départ pour vérifier la convergence vers la même
λ(k ) ← R(x(k ) ) = T
x(k ) Mx(k )
(k ) valeur.
normalisation x(k ) ← xα
fin On résoud ensuite Ax(k +1) = x(k ) et on norme
x(k +1) . Il s’agit bien de la méthode des puissances,
alg. 2.4 - Itérations inverses appliquée à A−1 , qui commence donc par converger
vers la valeur propre de A de module le plus élevé.
Bien évidemment, il est simple de démontrer que le taux de convergence est tel que :

|λn−1 |
τ ≈ ln (2.24)
|λn |
Comme pour la méthode des puissances, on converge d’autant mieux que les valeurs propres
sont séparées. . . d’où la méthode du paragraphe suivant.

2.5 Méthode des itérations inverses avec décalage spectral

Appelée encore méthode de Wielandt, cette méthode suppose que l’on connaît une approxi-
mation µ de la valeur propre cherchée λ.

2.5.1 Principe et convergence

On utilise la matrice décalée A⋆ = A − µI, on s’intéresse donc à A⋆ x = λ⋆ x et λ⋆ = λ − µ.
Comme on converge vers la valeur propre |λ⋆ | la plus petite, on va donc trouver λi la plus
proche de µ. Le taux de convergence est :
  
|λi − µ| |λi − µ|
τ = − ln inf , (2.25)
|λi −1 − µ| |λi +1 − µ|

et il est par conséquent très élevé si µ est très proche de λi . Voici quelques remarques :
– µ est proche de λi ; le système A⋆ x(k +1) = x(k ) est très mal conditionné. En fait, la mé-
thode y est peu sensible : l’erreur en solution tend à s’aligner dans la direction du vec-
teur propre et est éliminée par la normalisation ;
– Cette méthode permet aussi d’accéder aux valeurs propres d’une région donnée du
spectre ;
– On peut construire de cette manière la méthode des itérations inverses de Rayleigh :
on change l’estimation µ par λ(k ) à chaque itération. On obtient alors une convergence
cubique :
€ Š
kx(k +1) − vk = O kx(k ) − vk3 (2.26)

Cette méthode est coûteuse puisque l’on factorise à chaque itération.

 2.6 Transformation de Householder 31

2.5.2 Exemples d’utilisation dans un code éléments finis

On considère le système Kx = ω2 Mx.

Problème 1
On veut trouver le mode proche de la pulsation ω0 . On procède alors de la façon suivante :

1. décalage : K⋆ ← K − ω20 M

2. factorisation LDLT de K⋆ ; soit N⋆ le nombre de termes strictement négatifs de D

3. x(0) aléatoire puis itérations inverses avec décalage de ω0 sur K⋆ x = λ⋆ Mx pour obtenir
(λ⋆ , x)

4. ω2 ← ω20 + λ⋆

5. le numéro du mode trouvé est N⋆ si ω2 < ω20 ou N⋆ + 1 si ω2 ¾ ω20 .

Problème 2
On veut trouver tous les modes dans l’intervalle [ω1 , ω2 ]. On procède alors de la façon sui-
vante :

1. pour trouver le nombre de modes dans [ω1 , ω2 ], on compte le nombre de termes stric-
tement négatifs de D dans les factorisations de K − ω21 M et K − ω22 M ;

2. on procède par itérations inverses avec décalage et déflations orthogonales succes-

sives.

2.6 Transformation de Householder

Cette méthode date de 1958 [35].

2.6.1 Factorisation QR
On cherche à mettre A sous la forme A = QR où R est triangulaire supérieure et Q orthogo-
nale, par applications successives de transformations de Householder à A → HA [36].
H est de la forme H = I − 2uuT avec kuk = 1. Il s’agit cette fois d’une symétrie par rapport
à un plan de normale u comme illustré sur la figure 2.5. En conséquence, H est symétrique,
orthogonale, et conserve les valeurs propres.
Ces transformations de Householder sont appli-
h i
u quées sur les colonnes de A = a1 a2 · · · an suc-
a1
cessivement. Par exemple, sur la première, on veut
Ha1 = αe1 d’où le choix :

α = ka1 k = kHa1 k > 0 (2.27)

a1
e1
alors :
fig. 2.5 - Transformation de
Householder de a1 a1 − 2u(uT a1 ) = αe1
(2.28)
⇒ a1 − αe1 = 2u(uT a1 )
 32 SYSTÈMES AUX VALEURS PROPRES

En prenant le produit scalaire par a1 , on obtient :

α2 − αaT1 e1 = 2(uT a1 )2 (2.29)

d’où :
Ç
T 1
u a1 = α(α − aT1 e1 ) (2.30)
2
On peut noter que le terme (α − aT1 e1 ) est non nul si a1 n’est pas déjà parallèle à e1 . Finale-
ment, on en déduit le vecteur cherché :
a1 − αe1
u= (2.31)
2uT a1
Le procédé se poursuit en considérant ensuite la deuxième colonne de A et H telle que :
 
I 0
H=  (2.32)
0 H2

Le vecteur suivant u2 est déterminé à partir de e2 et a2 . C’est une méthode directe en n − 1

étapes.
On peut regarder les propriétés de remplis-
sage de A lors des transformations mais c’est as- boucle sur k = 1, 2, . . . , n − 1
(k ) T (k )
δ ← Ak :n ,k Ak :n ,k
sez délicat. En général, on ne stocke pas la matrice p
α← δ
Q ; l’application de Q à un vecteur est coûteuse β← 1
(k )
δ−αAk ,k
(4n opérations), cette méthode n’est donc pas (k ) (k )
vk :n ← Ak :n ,k
(k ) (k )
appliquée à la résolution de systèmes linéaires. vk ← vk − α
(k ) T (k )
zTk +1:n ← βvk :n Ak :n ,k +1:n
Cet algorithme se trouve aussi sous le nom de (k +1) (k ) (k )
Ak :n ,k +1:n ← Ak :n ,k +1:n − vk :n zTk +1:n
Householder-Businger-Golub. Rk ,k ← α
(k +1)
Rk ,k +1:n ← Ak ,k +1:n
La méthode est présentée par l’algorithme 2.6. fin
En ce qui concerne le stockage, vk écrase A1:k ,k
en mémoire, Rk +1,k +1:n écrase Ak +1,k +1:n . Rk ,k doit alg. 2.6 - Transformation de
Householder pour
être stocké ailleurs pour chaque itération, d’où un la factorisation QR
vecteur de dimension n supplémentaire. Q n’est
en général pas stockée, elle pourrait être reconstruite à partir des vk . Enfin, z est un vecteur
de stockage supplémentaire.
La complexité de cet algorithme pour une matrice pleine de taille n est de l’ordre de 23 n 3 .

2.6.2 Détermination des valeurs propres par la méthode QR

Cette méthode est attribuée à Francis et Vera N. Kublanovskaya en 1961. On considère le
problème Ax = λx avec des valeurs propres simples. On procède alors de la façon suivante :

1. A(1) = A ;

2. factorisation QR de A(k ) = Q(k ) R(k ) ;

3. calcul de A(k +1) = R(k ) Q(k ) ;

4. itération suivante.
 2.6 Transformation de Householder 33

Propriétés
– si les valeurs propres sont toutes simples, A(k ) tend vers la matrice diagonale des va-
leurs propres. Q(k ) tend vers la matrice des vecteurs propres (propriété peu utilisée) ;
– il y a conservation des bandes de A. C’est une méthode performante quand A est tri-
diagonale de petite taille. . . le coût d’une itération est alors d’environ 10n. Quand A est
pleine, une factorisation coûte ≈ 23 n 3 opérations et un produit QR, ≈ 12 n 3 opérations.
On peut aussi appliquer des translations (décalages) :

A(k ) − α(k ) I = Q(k ) R(k )

(2.33)
A(k +1) = R(k ) Q(k ) + α(k ) I

Les matrices restent toujours semblables, donc conservent les valeurs propres.
(k )
Concernant la convergence de la méthode, avec des décalages, et en prenant α(k ) = An,n :

|λn − A(k
n,n
)
| → 0 (2.34)
k →∞

alors :

k(A − λk +1 I)q(k ) k = O (k(A − λk I)q(k −1) k3 ) (2.35)

Sans translation, la convergence reste linéaire.

2.6.3 Forme de Hessenberg

Comme son nom l’indique, cette méthode est due à Gerhard Hessenberg (1874-1925) pour
qui on pourra consulter [37] et est dédiée aux matrices A réelles symétriques.
On cherche cette fois à rendre A tri-diagonale. Il
s’agit de la même technique colonne par colonne que
0 précédemment. On prend ici une transformation qui
ne s’applique qu’à la deuxième colonne pour la pre-
mière étape :
0
ak  
I 0
H=  (2.36)
0 H̃

fig. 2.7 - Matrice en cours de mais comme A est symétrique, il faut appliquer les
tri-diagonalisation
transformations à gauche et à droite :

H−1 AH = HT AH = HAH (2.37)

A(k ) prend donc la forme de la figure 2.7 et le principe consiste à travailler sur le vecteur
h iT
(k ) (k )
Ak ,k +1 · · · Ak ,n pour le rendre colinéaire à ek +1 .
Cette méthode a la propriété intéressante de ne pas modifier les valeurs propres. C’est
une méthode directe en n − 2 étapes et A reste symétrique. Travaillons sur le premier itéré :
 34 SYSTÈMES AUX VALEURS PROPRES

ploplo
boucle sur k = 1, 2, . . . , n − 1  
(k ) T (k )
δ ← Ak +1:n ,k Ak +1:n ,k A1,1 ãT1
p
α ← δ€ Š
A (1)
=  (2.38)
β ← 1/ δ − αAk +1,k
(k ) ã1 Ã(1)
(k ) (k )
vk +1:n ← Ak +1:n ,k on prend :
(k ) (k )
v k +1← v k +1 − α
(k ) (k )
zk +1:n ← βAk +1:n ,k +1:n vk +1:n H̃(1) = 1 − 2ũũT (2.39)
(k ) T
γ ← β/2vk +1:n zk +1:n
zk +1:n ← zk +1:n − γvk +1:n
(k ) en souhaitant avoir :
(k +1) (k )
Ak +1:n ,k +1:n ← Ak +1:n ,k +1:n
(k ) (k ) T
H̃(1) ã1 = αẽ1 (2.40)
−vk +1:n zTk +1:n − zk +1:n vk +1:n
fin En fait, ici ẽ1 = e2 et on procède comme précédem-
ment pour avoir ũ. On cherche donc H̃(1) = I − βỹỹT
alg. 2.8 - Transformation
1
de Householder avec α = ãT1 ã1 , β = α(α−A2,1 )
et ỹ = ã1 −αẽ1 . On calcule
ensuite la transformée par :
H̃(1) Ã(1) H̃(1) = (1 − βỹỹT )Ã(1) (1 − βỹỹT )
(2.41)
= Ã(1) − βỹ(ỹT Ã(1) ) − β(Ã(1) ỹ)ỹT + β2 ỹ(ỹT Ã(1) ỹ)ỹT
β
et en utilisant les intermédiaires z̃ = βÃ(1) ỹ et t̃ = z̃ − 2 (ỹT z̃)ỹ :

H̃(1) Ã(1) H̃(1) = Ã(1) − ỹz̃T − z̃ỹT + β(z̃T ỹ)ỹỹT

(2.42)
= Ã(1) − ỹt̃T − t̃ỹT
C’est la technique mise en œuvre dans l’algorithme 2.8, dans lequel vk n’est en général pas
conservé et où z est un vecteur de stockage supplémentaire. La réduction s’effectue en n − 1
étapes.

2.7 Exemple d’utilisation dans un code éléments finis

On considère Kx = λMx avec K symétrique, définie positive, M symétrique positive pour
l’obtention de m valeurs propres, en basse fréquence.

2.7.1 Réduction dans un sous-espace

On cherche à trouver toutes les m valeurs propres dans un sous-espace de dimension p ¾
h i
m , défini par p vecteurs rangés dans la matrice de base V = v1 · · · vp . Éventuellement,
les vi sont orthonormés au sens de K, soit VT KV = K̂ = I.
On cherche alors les vecteurs propres x sous la forme x = Vx̂, x̂ est de dimension p . Le
problème est alors KVx̂ = λ̂MVx̂, d’où :

(VT KV) x̂ = λ̂ (VT MV) x̂ (2.43)

| {z } | {z }
K̂ M̂
Le problème K̂x̂ = λ̂M̂x̂ a été construit par une réduction de Ritz. Sa résolution complète
permet de trouver p caractéristiques propres (λ̂, x̂).
On retourne alors dans l’espace de départ où x = Vx̂ est l’approximation d’un vecteur
propre, λ̂ est l’approximation de la valeur propre associée, approximations dont la qualité
dépend du sous-espace choisi. . . Ce choix va être discuté dans les paragraphes suivants.
 2.8 Méthode de Lanczos 35

2.7.2 Technique itérative

Rapidement, il s’agit de :

1. choisir p vecteurs de départ ;

2. réduire à K̂x̂ = λ̂M̂x̂ dont on cherche les valeurs propres (si on a une approximation
des valeurs propres, on fait un décalage) par une méthode QR ;

3. itérer en sens inverse avec décalage (souvent une seule itération suffit) pour avoir les
approximations des vecteurs propres associés dans l’espace de départ ;

4. passer à l’itération suivante.

Comme la convergence, dépend du rapport des modules des valeurs propres, on a intérêt à
conserver p > m approximations. La règle conseillée est choisir une taille du sous-espace :

p = min(2m , m + 8) (2.44)

|λi |
La convergence pour une valeur propre λi est proportionnelle à |λp +1 |
.

2.8 Méthode de Lanczos

Cette méthode est due à Lanczos [38] en 1950. On travaille sur un sous-espace généré par
h i
la suite des itérés inverses d’un vecteur v(0) de départ : v(0) A−1 v(0) · · · et sur une base
orthonormée.
La convergence est très rapide : 3 à 4 itérations inverses par valeur propre convergée, et
elle est indépendante de n ; cependant, la taille du sous-espace croît au cours des itérations.
Un inconvénient est la perte rapide d’ortho-
gonalité sur les vecteurs successifs (elle se pro- v(0) ← 0, β(0) ← 0
duit à chaque fois qu’une valeur propre converge), Vecteur de départ v(1) , kv(1) k = 1
tant que R > ǫ, itérer sur k
d’où des instabilités numériques. En pratique, résoudre Az = v(k )
α(k ) ← v(k ) T z
pour pouvoir traiter de grands problèmes, il
d ← z − αv(k ) − β(k −1) v(k −1)
est nécessaire d’utiliser des stratégies de ré- β(k ) ← p1
dT d

orthogonalisation. Ces dernières ne seront pour- v(k +1) ← β(k ) d

recherche de (x̂, λ̂) tel que Tx̂ = λ̂x̂
tant pas décrites ici. fin
En termes simples, l’algorithme se résume à
alg. 2.9 - Méthode de Lanczos
choisir un vecteur de départ v(0) puis de résoudre le
système Az(k ) = v(k ) , d’orthogonaliser z(k ) par rap-
z(k )
port aux v(j ) , j ¶ k puis de normer v(k +1) ← kz(k ) k
.
Pour démontrer le taux de convergence, il suffit de chercher v(k +1) sous la forme :

β(k ) v(k +1) = A−1 v(k ) − α(k ) v(k ) − β(k −1) v(k −1) (2.45)

ce que l’on peut montrer par récurrence, avec v(k ) T v(j ) = 0 pour j = k − 1 et j = k − 2,
et v(k ) T v(k ) = 1, pour avoir orthogonalité par rapport à tous les v(j ) , j < k . Ces conditions
36 SYSTÈMES AUX VALEURS PROPRES

donnent les valeurs de α(k ) = v(k ) A−1 v(k ) et β(k −1) = v(k ) A−1 v(k −1) . Matriciellement, en no-
tant V de taille (n, k ) la matrice qui contient tous les v(j ) disponibles à l’itération k , la rela-
tion (2.45) est équivalente à :
 
α(0) β(0) 0
 
 .. ..  h i
 β(0) . . 
A V − V
−1

 = 0 · · · 0 β(k )
 v (k +1) (2.46)
.. ..
 . . β(k −1) 
 
0 β(k −1) α(k )
| {z }
T

ce qui implique :

VT A−1 V − T = 0 (2.47)

Par conséquent, T a les mêmes valeurs propres que VT A−1 V (voir le principe des sous-espaces).
On cherche ensuite les caractéristiques propres de T, Tx̂ = λ̂x̂, λ̂ sont des approximations
des valeurs propres de A−1 et Vx̂, celles des vecteurs propres associés.
Afin d’illustrer le comportement de cette méthode, la figure 2.10 présente l’estimation
des valeurs propres aux cours des itérations de la méthode de Lanczos sur le problème test
de Wilkinson (de taille 21) comportant plusieurs valeurs propres multiples et plusieurs va-
leurs propres très proches (les cercles sont les valeurs propres exactes). On voit que la multi-
plicité et/ou proximité des valeurs propres échappe à l’algorithme, et que les valeurs propres
se mettent à converger quasiment l’une après l’autre.

12

10

8
valeurs propres

−2
0 2 4 6 8 10 12 14
itérations de Lanczos

fig. 2.10 - Comportement de la méthode de Lanczos sur le test de Wilkinson

Critère d’arrêt
Classiquement, on borne la norme du résidu r = A−1 Vx̂ − λ̂Vx̂ avec :
h i
A−1 Vx̂ = VTx̂ + 0 · · · 0 β(k ) v (k +1) x̂ (2.48)

avec Tx̂ = λ̂x̂, d’où r = β(k ) v(k +1) x̂k .

L’algorithme 2.9 décrit cette méthode, pour lequel les α(k ) et les β(k ) sont les composantes
de T.
 3
RÉSOLUTION DE SYSTÈMES
NON-LINÉAIRES PARTICULIERS

Dans cette partie, des préoccupations directement issues du calcul de structure par élé-
ments finis sont abordées. Il s’agit ici d’appliquer les outils précédents à certains problèmes
particuliers, et non pas de disposer d’un cours entier consacré aux éléments finis. Il en sera
de même dans la partie 4.4 pour les relations de comportement du matériau. Pour une pré-
sentation dédiée aux éléments finis, le lecteur est invité à se reporter, par exemple, aux nom-
breux cours de qualité disponibles sur le réseau. Par exemple, ceux de C. Felippa : Nonlinear
Finite Element Methods, essentiellement pour les non-linéarités d’origine géométriques, et
Advanced Finite Element Methods, pour la technologie des éléments finis, donnés au Aeros-
pace Engineering Sciences Department de l’University of Colorado at Boulder.

3.1 Non-linéarités matérielles de type plasticité

Le problème à résoudre est celui du calcul de structures dont le matériau possède un com-
portement non linéaire. Il peut par exemple s’agir d’un matériau élastique non linéaire (sans
dissipation, donc) ou d’un problème d’évolution comme la plasticité ou la visco-plasticité.
Dans ce dernier cas, une approche incrémentale transforme le problème d’évolution en une
série de problèmes non-linéaires stationnaires à chaque incrément de temps. Ceux-ci sont
alors résolus par exemple par des méthodes classiques de type Newton [39], décrites ici.
Pour plus d’informations, on pourra consulter [40].

3.1.1 Formulation du problème

Le problème de référence consiste à vérifier :
– les équations d’admissibilité cinématique : les déformations ǫ doivent être compa-
tibles avec un champ de déplacement admissible U (régulier et satisfaisant les condi-
tions aux limites en déplacement imposé)

ǫ = ǫ(U) (3.1)

– les équations d’admissibilité statique : le champ de contraintes σ doit être en équilibre

avec les charges extérieures ;
 38 RÉSOLUTION DE SYSTÈMES NON-LINÉAIRES PARTICULIERS

– les relations de comportement [41]. On considère un matériau dont l’état est décrit
par (outre la déformation ou la contrainte) un jeu de variables internes X, avec leurs
variables thermodynamiques associées Y.
La partition de la déformation en partie élastique et non élastique s’écrit :

ǫ = ǫ e + ǫp (3.2)

Les lois d’état reliant les variables cinématiques et statiques sont associées à l’énergie
libre (qui dépend de l’état) ψ(ǫe , X) :

∂ψ
σ=
∂ ǫe (3.3)
∂ψ
Y=
∂X
L’évolution de l’état, pour un modèle standard associé, se déduit du potentiel de dissi-
pation Φ(σ, Y) par la loi de normalité :

∂Φ
ǫ̇p = λ̇
∂σ (3.4)
∂Φ
−Ẋ = λ̇
∂Y
avec, dans le cas non visqueux, Φ ¶ 0, λ̇ ¾ 0 et Φλ̇ = 0. λ est le multiplicateur plastique.

3.1.2 Résolution pour un incrément de temps

Pour les problèmes qui nous préoccupent, à chaque incrément de temps, on est amené à
résoudre, après une discrétisation en déplacement, un problème de la forme :

r(u) = K(u)u − f = 0 (3.5)

r(u) est le gradient à annuler. La matrice de rigidité dépend de la solution u, par l’intermé-
diaire de la relation de comportement du matériau. Les méthodes de type Newton sont habi-
tuellement utilisées : si u0 est une approximation de la solution, cette dernière est u = u0 +δu
et un développement de Taylor de r donne au voisinage de u :

∂r
(u0 )δu = −r(u0 ) (3.6)
∂u
ce qui suggère un algorithme itératif de la forme :

ui +1 − ui = −s i Hi ri (3.7)
∂r
où ri = r(ui ), Hi est une approximation de l’inverse de la matrice jacobienne J(ui ) = ∂u
(ui ).
s i est une amplitude destinée à réduire ri +1 1 .
En fait, les méthodes présentées ci-après sont issues d’une classe de problèmes plus
vaste, celle de la minimisation de fonctionnelles non-linéaires f (U) ; dans ce cadre, f est
parfois appelée la fonction coût.
1. Dans un premier temps, on peut choisir s i = 1 mais une procédure de type line search peut être utilisée comme une
procédure itérative additionnelle pour améliorer sa valeur
 3.1 Non-linéarités matérielles de type plasticité 39

3.1.3 Méthode de la plus grande pente

Dans ce cas :

H=I (3.8)

3.1.4 Méthodes de gradient conjugué

Ces méthodes n’imposent pas le stockage des matrices Hi car seuls quelques vecteurs sont
nécessaires au calcul de Hi ri . Cela est réalisé avec une formule de récurrence sur une direc-
tion de recherche d :
(
ri pour i = 1
di = (3.9)
ri + ζi di −1 pour i > 1

et :

ui +1 = ui − s i di (3.10)

ζi est une scalaire utilisé pour conjuguer les directions de recherche, i.e. dTi di −1 = 0, pour
lequel deux expressions classiques sont :

1. Fletcher-Reeves [42] :

rTi ri
ζi = (3.11)
rTi −1 ri −1

2. Polak-Ribière [43] :

rTi (ri − ri −1 )
ζi = (3.12)
rTi −1 ri −1

Si la fonctionnelle à minimiser est quadratique, donc si K ne dépend pas de u, alors toutes

les directions de recherche successives sont conjuguées entre elles.
Pour mémoire, il exite d’autres méthodes de grandient conjugué, comme la méthode du
shortest residual [44].

3.1.5 Méthode de Newton-Raphson

Ici, s i = 1 et :
 −1  −1
∂r   ∂K 
Hi = J(ui )−1 =  (ui ) = K(ui ) + (ui ) (3.13)
∂u ∂u

à chaque itération, il faut évaluer la matrice jacobienne et résoudre le précédent système

linéaire.

3.1.6 Méthode de Newton-Raphson modifiée

Elle utilise un Jacobien figé, évalué à une précédente estimation de la solution u⋆ :

Hi = J(u⋆ )−1 (3.14)

 40 RÉSOLUTION DE SYSTÈMES NON-LINÉAIRES PARTICULIERS

3.1.7 Mises à jour de type Quasi-Newton

L’idée est de modifier simplement Hi à chaque itération plutôt que de le laisser fixé comme
MNR, et plutôt que de le recalculer entièrement, comme pour NR.

Hi +1 = Hi + δHi (3.15)

en s’assurant que la condition sécante :

Hi +1 (ui +1 − ui ) = ri +1 − ri (3.16)

est satisfaite. Elle correspond à une approximation au premier ordre de la matrice jaco-
bienne au voisinage de ui .
Si Hi est mis à jour, on parle de mise à jour inverse, alors qu’une mise à jour directe
modifie H−1
i .
Le rang de la mise à jour est le rang de la matrice de correction qui est utilisée. Pour
une mise à jour de rang 1, δH = abT où a et b sont deux vecteurs, la formule de Sherman-
Morrison-Woodbury [45, 46, 47] permet d’exprimer directement l’inverse, tant que H n’est
pas singulier et que 1 + bTH−1 a 6= 0 :
H−1 abT H−1
(H + abT )−1 = H−1 − (3.17)
1 + bTH−1 a
L’algorithme générique consiste à itérer sur :
– direction de recherche di = −Hi ri ;
– mins i >0 kri +1 k par exemple par une procédure de line search ;
– mise à jour de Hi .
Le but est de construire une approximation de J−1 (u), qui s’améliore au cours des itérations,
en utilisant l’information des gradients successifs. Par exemple, avec la précédente correc-
tion de rang 1, en vérifiant la symétrie et la condition sécante, on obtient :

(δ i − Hi γi )(δ i − Hi γi )T
Hi +1 = Hi + (3.18)
γTi (δ i − Hi γi )
où δi = ui +1 − ui et γi = ri +1 − ri .
Une propriété intéressante pour l’efficacité est d’avoir H symétrique, au moins quand
le problème à résoudre l’est, et définie positive pour avoir effectivement une direction de
descente telle que dTi ri < 0. Avec la précédente mise à jour de rang 1, il n’est pas garanti
d’avoir Hi défini positif.

Méthode de Davidson-Fletcher-Powell
Cette méthode de rang 2, proposée dans [48] préserve le caractère défini positif. En démar-
rant avec H0 symétrique, défini positif, la mise à jour est :

δi δTi (Hi γi )(Hi γi )T

HDFP
i +1
= Hi + − (3.19)
δTi δi γTi Hi γi
Dans le cas où la fonction coût est quadratique, les directions di sont K-orthogonales entre
elles ; après n itérations, Hn = J−1 . Si de plus H0 = I, il s’agit de la méthode du gradient
conjugué.
3.1 Non-linéarités matérielles de type plasticité 41

Méthode de Broyden-Fletcher-Goldfarb-Shanno
Cette méthode de rang 2, proposée dans [49], préserve aussi le caractère défini positif. En
démarrant avec H0 symétrique, défini positif, la mise à jour est :
  T
 δ i H i γ i  δ i H i γ i 
HBFGS DFP T
i +1 = Hi +1 + (γi Hi γi )  T − T  T − T 
δ i γi γi H γ
i i δ i γi γi H γ
i i
  (3.20)
T
 γT
i H i γi  δi δi δi γTi Hi + Hi γi δTi
= Hi + 1 + T  T −
δi γi δi γi δTi γi

Cette méthode est réputée pour être l’une des plus robustes [50].
La famille de méthodes dite de Broyden est définie par :

Hi +1 = αHBFGS DFP
i +1 + (1 − α)Hi +1 (3.21)

α étant tel que le caractère défini positif est assuré (0 ¶ α ¶ 1 le fait).

Versions avec mémoire limitée

Même si elles sont symétriques, les corrections présentées précédemment souffrent d’une
grande demande en mémoire, pour les problèmes éléments finis de grande taille (les ma-
trices Hi ne sont en effet pas creuses). Une solution consiste à ne stocker Hi uniquement
de façon implicite, c’est-à-dire à être capable de construire Hi ri avec H0 (qui lui est creux)
et des vecteurs successifs. De plus, le calcul de la direction de recherche est alors souvent
moins coûteux que dans le cas de l’utilisation explicite d’un produit matrice-vecteur.
Par exemple, le gradient conjugué préconditionné (PCG) par H0 , correspond à :

Hi +1 = Vi H0 (3.22)

avec :

γ j δTj
Vj = I − (3.23)
δTj γ j

et la méthode BFGS, à :
 T    T  
Y i Y i X i Y i T Y i
      δ k δ k 
Hi +1 =  V j  H0  Vj  +  Vj  T  Vj  (3.24)
j =0 j =0 j =0 k =j +1
δ γ
k k k =j +1

Une troncature peut alors être effectuée pour limiter encore les ressources en mémoire, en
éliminant les contributions des plus anciens termes. Elle correspond alors à la version limi-
ted memory de BFGS (L-BFGS) [51, 52].
 4
RÉSOLUTION DE
SYSTÈMES DIFFÉRENTIELS

Il est possible de consulter sur le sujet les références [53, 54]. Dans les applications qui nous
intéressent, on est conduit à résoudre de tels systèmes, par exemple dans les cas suivants :
intégration de relation de comportement, problème de thermique transitoire, problème de
dynamique rapide. . .

4.1 Équation scalaire du premier ordre

4.1.1 Problème de Cauchy
Il s’agit de trouver la fonction x (t ), t ∈ [0, T], telle que :
dx
= f (x , t )
dt (4.1)
x (0) = u 0

Théorème de Cauchy-Lipschitz
Il y a existence et unicité de la solution si f est assez régulière 1 en x :

∃k > 0, ∀t ∈ [0, T], ∀(y , z ) ∈ R2 , | f (z , t ) − f (y , t )| ¶ k |z − y | (4.2)

4.1.2 Ré-écriture du problème

Pour se ramener aux techniques d’intégration, il est plus intéressant d’écrire le problème
sous forme intégrale :
Zt
x = u0 + f (x , t ) dt (4.3)
0

Une discrétisation de l’intervalle de temps [0, T] est la séquence d’instants 0 = t 0 < t 1 < . . . <
t m = T. Le pas de temps est h = t i +1 − t i . Le principe est de chercher à approcher la valeur
de la fonction aux piquets de temps, x (t i ), par la valeur x i , ce qui permet de définir l’erreur
e i = |x (t i ) − x i |. Dans la suite, on notera f i = f (x i , t ). Le problème incrémental est, ayant la
solution approchée jusqu’à t i , de la déterminer à t i +1 avec :
Zt
x = xi + f (x , t ) dt (4.4)
ti

1. Plus précisément en termes mathématiques, f est continue et lipschitzienne.

 44 RÉSOLUTION DE SYSTÈMES DIFFÉRENTIELS

Exemple : schéma d’Euler explicite

La figure 4.3(a) illustre l’estimation de l’intégrale cherchée. Le schéma s’écrit :

x i +1 = x i + h f i
(4.5)
x0 = u 0

4.1.3 Schémas explicites à un pas

Ce sont les schémas de la forme :

x i +1 = x i + hΦ(x i , t i , h)
(4.6)
x0 = u 0

convergence
Si on a la propriété supi |e i | →h→0 0 (en l’absence d’erreur d’arrondi !)

stabilité
On considère le schéma perturbé :

x̃ i +1 = x̃ i + hΦ(x̃ i , t i , h) + ǫi +1
(4.7)
x̃ 0 = u 0 + ǫ0

Le schéma considéré est stable ssi :

m
X
∃S > 0, sup |x̃ i − x i | ¶ S |ǫi | (4.8)
i i =0

consistance
Pour définir cette notion, on a besoin de l’erreur de consistance (local truncation error) :


c i +1 = x (t i +1 ) − x (t i ) + hΦ(x (t i ), t i , h) (4.9)

Le schéma est consistant ssi :

X
|c i | → 0 (4.10)
h→0
i

ordre
On dit d’un schéma qu’il est d’ordre p quand ∃M > 0, supi |c i | ¶ Mh p +1 .
Les trois premières notions précédentes sont illustrées sur la figure 4.1.

h →0

h →0
t t t

(a) convergence (b) consistance (c) stabilité

fig. 4.1 - Propriétés de base d’un schéma d’intégration

 4.1 Équation scalaire du premier ordre 45

Propriétés
Il existe des liens entre les différentes caractéristiques précédentes, en particulier :
– stabilité et consistance ⇔ convergence
– Φ lipschitzienne en x ⇒ stabilité
– consistance ⇔ Φ(x , t , 0) = f (x , t ) ⇔ d’ordre p ¶ 1
– d’ordre p ⇒ |e i | ¶ Nh p

Exemple : schéma d’Euler explicite

Dans ce cas, Φ(x , t , h) = f (x , t ) et par conséquent, le schéma est consistant et stable, donc
convergent :


c i +1 = x (t i +1 ) − x (t i ) + h f (x (t i ), t i ) = O (h 2 ) (4.11)

or :

x (t i +1 ) = x (t i ) + h ẋ (t i ) + O (h 2 )
(4.12)
ẋ (t i ) = f (x (t i ), t i )

donc c’est un schéma d’ordre 1.

4.1.4 Quelques problèmes pratiques

Influence des erreurs d’arrondi
Perturbons le schéma :

x̃ i +1 = x̃ i + hΦ(x̃ i , t i , h) + ǫi +1
(4.13)
x̃ 0 = u 0 + ǫ0

Si le schéma est convergent, e i = |x (t i ) − x i | → 0. Qu’en est-il de ẽ i = |x (t i ) − x̃ i | ?

h→0

ẽ i = |x (t i ) − x i + x i − x̃ i | ¶ |x (t i ) − x i | + |x i − x̃ i | (4.14)
| {z } | {z }
① ②
avec, dans l’équation (4.14) :

− le schéma est d’ordre p donc le terme ① ¶ Nh p ;

”Pm Pm —
− le schéma est stable donc le terme ② ¶ S i =0 |ǫ| + i =1 h|ǫ| .

Au bilan :

ẽ i ¶ S[(m + 1)|ǫ| + hm |ǫ|] + Nh p (4.15)

avec hm = T, d’où :
 
T
ẽ i ¶ S|ǫ| 1 + T + + Nh p (4.16)
h
Il n’est donc pas forcément avantageux de réduire trop fortement le pas de temps pour la
précision du schéma. En particulier, on peut dire que si le pas de temps n’est pas grand de-
vant les erreurs numériques d’évaluation de la fonction Φ, la solution est entachée d’erreurs
par cumul des erreurs numériques.
 46 RÉSOLUTION DE SYSTÈMES DIFFÉRENTIELS

Problème numériquement mal posé

Considérons le problème suivant :
dx
= 4x − 1
dt (4.17)
1
x (0) =
4
dont la solution est :
1
x (t ) ≡ (4.18)
4
Perturbons maintenant les conditions initiales :
1
x (0) = +ǫ (4.19)
4
La nouvelle solution est alors très différente de la précédente :
1
x (t ) = + ǫe 4t (4.20)
4
Ce problème est numériquement mal posé. On risque d’avoir des difficultés à le résoudre
car il est très sensible aux conditions initiales 2 .

4.1.5 A-stabilité
De façon pratique, la propriété de stabilité d’un schéma n’est pas suffisante. Prenons par
exemple le schéma d’Euler explicite déjà décrit comme stable, et résolvons le problème :
dx
= −ω2 x
dt (4.21)
x (0) = 1
dont la solution exacte est :
2t
x (t ) = e −ω (4.22)

La figure 4.2 présente les solutions obtenues numériquement en même temps que la solu-
tion exacte pour deux pas de temps. Ces solutions diffèrent fortement, et de plus en plus
au cours du temps, si le pas de temps n’est pas assez petit. On parlera alors d’un schéma
conditionnellement stable. Un schéma inconditionnellement stable 3 [55] se définit pour
les équations linéaires. On teste alors le schéma sur un problème du type :
dx
= cx (4.23)
dt
où c est une constante complexe et pour lequel le cas où le schéma est tel que x i +1 = r (c h)x i ,
où r est un polynôme ou une fraction rationnelle. On définit alors le domaine de A-stabilité
comme suit :

Dstab = {z ∈ C, |r (z )| ¶ 1} (4.24)

Un schéma est alors dit A-stable (ou inconditionnellement stable) si et seulement si :

{z ∈ C, Re (z ) ¶ 0} ⊂ Dstab (4.25)
2. on peut d’ailleurs se poser des questions quant à la notion même de solution pour un tel problème
3. ou A-stable pour absolute stability
 4.1 Équation scalaire du premier ordre 47

1
0,8
2
0,6
0,4 1
0,2

0,2 0,4 0,6 0,8 1 0,2 0,4 0,6 0,8 1

−0,2
−1
−0,4
−2

(a) h = 1/10 (b) h = 1/7

fig. 4.2 - Stabilité du schéma d’Euler explicite

Exemple : schéma d’Euler explicite

La relation est la suivante :

x i +1 = x i + hc x i = (1 + hc )x i ⇒ r (z ) = 1 + z (4.26)

La condition de stabilité est donc :

|1 + hc | ¶ 1 (4.27)

2
Pour c = −ω2 ∈ R− , elle s’écrit h ¶ ω2
. Il s’agit donc d’un schéma conditionnellement stable.
Pour l’exemple précédent, avec c = −15 (la solution tend vers 0), le critère de stabilité donne
1
0¶h ¶ 7,5
. Ceci est bien en accord avec les résultats de la figure 4.2.

La condition de stabilité dépend bien entendu du schéma, mais aussi du problème. Son
interprétation est que, pour une solution bornée, il faut que le schéma donne une ap-
proximation bornée.

4.1.6 Méthode des trapèzes généralisée

Cette méthode est décrite avec un paramètre α ∈ [0, 1] par :
 
x i +1 = x i + h (1 − α)f i + α f i +1
(4.28)
x0 = u 0

Le calcul de x i +1 fait intervenir α f i +1 = α f (x i +1 , t i +1 ). Il s’agit donc d’une méthode implicite

si α 6= 0.
Le cas α = 0 conduit à la méthode d’Euler explicite ; dans le cas où α = 1, il s’agit de
la méthode d’Euler implicite ; enfin si α = 12 , c’est la méthode des trapèzes, appelée aussi
méthode de Crank-Nicolson et qui date de 1947. La figure 4.3 illustre l’approximation de
l’intégration réalisée par chacun de ces cas particuliers.
 48 RÉSOLUTION DE SYSTÈMES DIFFÉRENTIELS

f f f

t t t
ti t i +1 ti t i +1 ti t i +1
(a) α = 0 (b) α = 1 (c) α = 0,5

fig. 4.3 - Différents types d’intégration

Stabilité
On considère le problème :
ℑ(hc )
dx
= cx (4.29)
α = 1/2

α > 1/2 dt
α < 1/2 α= le schéma donne :
1
α=0
ℜ(hc ) x i +1 = x i + h[(1 − α)c x i + αc x i +1 ] (4.30)
−1 1
donc :

x i +1 = (1 − αhc )−1 (1 + (1 − α)hc ) x i (4.31)

| {z }
fig. 4.4 - Domaines de stabilité pour la B
méthode des trapèzes généra-
lisés Il faut avoir ̺(B) < 1, ce qui est le cas si
toutes les valeurs propres de B sont simples,
c’est-à-dire si la condition suivante est vérifiée :

|1 + (1 − α)hc |
¶1 (4.32)
|1 − αhc |
1
− si α ¾ 2
le schéma est inconditionnellement stable ;
1 2
− si α < 2
avec c = −ω2 , la condition de stabilité s’écrit h ¶ (1−2α)ω2
.

Les différents domaines de stabilité correspondants sont tracés sur la figure 4.4.

Ordre
En utilisant la même méthode que celle déjà employée pour le schéma d’Euler explicite, on
peut facilement montrer que le schéma des trapèzes généralisés est d’ordre 1 si α 6= 12 , et
qu’il est d’ordre 2 pour α = 12 , c’est-à-dire pour Crank-Nicolson.

4.1.7 Exemple de méthode à pas multiples : Adams

L’idée des méthodes à pas multiple est d’utiliser non seulement l’évaluation f i mais aussi
un certain nombre de valeurs passées f i −1 . . .
À titre d’illustration, la méthode d’Adams, datant de 1883, utilisant un développement
de Taylor de la solution est maintenant décrite.
4.1 Équation scalaire du premier ordre 49

Adams explicite
On utilise pour cela un développement de Taylor (avant) de x (t i +1 ) :
h2
x (t i +1 ) = x (t i ) + h ẋ (t i ) + ẍ (t i ) + . . . (4.33)
2!
df
avec ẋ (t i ) = f (x (t i ), t i ) et ẍ (t i ) = dt
(x (t i ), t i ). On garde ensuite les termes d’un ordre suffi-
sant pour obtenir l’ordre voulu du schéma. Par exemple, pour Adams explicite d’ordre 1, on
conserve les termes en O (h) ; on obtient x i +1 = x i + h f i + O (h 2 ), c’est Euler explicite !
Pour Adams explicite d’ordre 2, on conserve les termes en O (h 2 ) qui impliquent la dé-
rivée première de f ; pour l’estimer on utilise encore le développement d’une fonction g :

g (t i +1 ) = g (t i ) − h ġ (t i ) + O (h 2 ) (4.34)

qui permet d’avoir l’expression de la dérivée ġ :

1
ġ (t i ) = [g (t i ) − g (t i +1 )] + O (h) (4.35)
h
df
Avec ce procédé appliqué à dt
, on obtient la méthode cherchée :

h
x i +1 = x i + [3 f i − f i −1 ] + O (h 3 ) (4.36)
2
La démarche se prolonge pour les ordres de schéma supérieurs. On obtient ainsi les for-
mules d’Adams-Bashforth :

x i = x i −1 + h f i −1
h
x i = x i −1 + (3 f i −1 − f i −2 )
2
h
xi = x i −1 + (23 f i −1 − 16 f i −2 + 5 f i −3 )
12
h (4.37)
xi = x i −1 + (55 f i −1 − 59 f i −2 + 37 f i −3 − 9 f i −4 )
24
h
xi = x i −1 + (1901f i −1 − 2774 f i −2 + 2616 f i −3 − 1274 f i −4 + 251 f i −5 )
720
h
xi = x i −1 + (4277f i −1 − 7923 f i −2 + 9982 f i −3 − 7298 f i −4 + 2877 f i −5 − 475 f i −6 )
1440
On ne peut cependant pas augmenter l’ordre du schéma impunément : l’inconvénient
réside dans un rétrécissement du domaine de stabilité. Celui-ci est illustré sur la figure 4.5.

Adams implicite
L’idée est la même que précédemment, en utilisant un développement de Taylor arrière de :
h2
x (t i ) = x (t i +1 − h) = x (t i +1 ) − h ẋ (t i +1 ) + ẍ (t i +1 ) + . . . (4.38)
2!
La méthode d’Adams implicite d’ordre 1 est ainsi :

x i = x i +1 − h f i +1 + O (h 2 ) (4.39)

d’où :

x i +1 = x i + h f i +1 + O (h 2 ) (4.40)
 50 RÉSOLUTION DE SYSTÈMES DIFFÉRENTIELS

3
1, 5
2
1 k =4

0, 5 1
k =3
k =1 k =2 k =2 k =3 k =4
ℑ(z )

ℑ(z )
0

−0, 5 −1

−1 −2

−1, 5 −3

−2, 5 2 1, 5 −1 −0, 5 0 0, 5 1 −7 −6 −5 −4 −3 −2 −1 0 1
ℜ(z ) ℜ(z )

(a) Adams explicite : Adams-Bashforth (b) Adams implicite : Adams-Moulton

fig. 4.5 - Domaine de stabilité dans le plan complexe

c’est Euler implicite !

L’utilisation de la même technique que pour Adams explicite permet de construire les
schémas implicites d’ordres plus élevés, et on obtient alors les formules d’Adams-Moulton :

x i = x i −1 + h f i
h
x i = x i −1 + (f i + f i −1 )
2
h
xi = x i −1 + (f i + 8 f i −1 − f i −2 )
12
h (4.41)
xi = x i −1 + (9 f i + 19 f i −1 − 5 f i −2 + f i −3 )
24
h
xi = x i −1 + (251 f i + 646 f i −1 − 264 f i −2 + 106 f i −3 − 19 f i −4 )
720
h
xi = x i −1 + (475 f i + 1427 f i −1 − 798 f i −2 + 482 f i −3 − 173 f i −4 + 27 f i −5 )
1440
Le fait d’utiliser des schémas implicites augmente la stabilité, mais la même tendance est
observée lorsque l’ordre croît, voir figure 4.5.

Amorçage du schéma
Aux premiers pas de temps, on ne dispose pas forcément de tous les f i −q . Dans ce cas, on
peut utiliser une méthode à 1 pas pour démarrer l’intégration jusqu’à avoir suffisamment
progressé pour pouvoir appliquer la méthode à pas multiples voulue.

4.1.8 Exemple de méthode de prédiction correction basée sur Adams

Le problème des méthodes implicites réside dans le fait, qu’en non-linéaire au moins, on
ne peut exprimer explicitement la valeur au pas courant. On est alors obligé de mettre en
œuvre une méthode itérative, généralement de type point fixe.
(0)
On part ainsi de la version implicite ; par exemple Adams d’ordre 2 : connaissant x i +1 ,
(0)
on évalue f (x i +1 , t i +1 ) qui permet, en utilisant la formule d’intégration d’Adams implicite,
 4.1 Équation scalaire du premier ordre 51

(1)
d’atteindre une valeur x i +1 . Cette dernière étape est appelée correction. On itère ainsi jusqu’à
(0)
une convergence souhaitée. Pour avoir la première valeur x i +1 , la prédiction, on utilise la
version implicite de la méthode, par exemple ici, Adams explicite d’ordre 2. En pratique,
une à trois itérations suffisent.
On peut aussi apporter une amélioration à faible surcoût avec une correction du prédic-
teur avec l’erreur de troncature. Pour Adams d’ordre 2, on procède comme suit :
1. prédicteur (AE2) :
(0) h
x i +1 ← x i + (3 f i − f i −1 ) (4.42)
2
2. modificateur :
(0) (0) 1  (0)

x̂ i +1 ← x i +1 − xi − xi (4.43)
12
3. correcteur (AI2) :
h€ˆ Š
x i +1 ← x i + f i +1 + f i −1 (4.44)
2
 
(0) 1
où fˆi +1 = f x̂ i +1 , t i +1 . Le coefficient 12
est celui issu de la formule de (AI3).

4.1.9 Méthode de Runge-Kutta

Historiquement présentée dans [56] en 1901, cette méthode est de la forme :

x i +1 = x i + hΦ(x , t , h)
(4.45)
x0 = u 0
avec :
q
X
Φ(x , t , h) = a l k l (x , t , h)
l =1

k 1 (x , t , h) = f (x , t ) (4.46)
 
X l −1
 (l ) 
k l (x , t , h) = f x + β j k j (x , t , h), t + αl h 
j =1

L’idée est d’intégrer en utilisant des valeurs en des points particuliers, avec des poids d’inté-
(l )
gration particuliers. Les constantes a l , αl et β j sont ajustées pour obtenir l’ordre souhaité.

Exemple de schéma à l’ordre 2

On prend Φ = a f (x , t ) + b f (t + αh,x + βh f (x , t )) et on détermine a , b , α, β pour avoir :
1   
x (t + h) − x (t ) − Φ x (t ), t , h = O (h 2 ) (4.47)
|h {z } | {z }
① ②
avec, dans l’équation (4.47) :
 
h ∂f ∂f
①= f + + f + O (h 2 )
2 ∂t ∂x
 
∂f ∂f 2
② = a f + b f + αh + βh f + O (h )
∂t ∂x
 52 RÉSOLUTION DE SYSTÈMES DIFFÉRENTIELS

Il faut donc avoir :

h h
0 = 1−a −b 0= − b αh 0= − b βh (4.48)
2 2
1
d’où les coefficients cherchés a = 1 − b , α = β = 2b et le schéma :
 
h h
Φ(x , t , h) = (1 − b )f (x , t ) + b f x + f (x , t ), t + (4.49)
2b 2b
où b est une constante arbitraire non nulle. Des cas particuliers existent :
– b = 1 : il s’agit de la méthode d’Euler modifiée, d’ordre 2 ;
1
– b= 2
: il s’agit de la méthode de Heun, d’ordre 2.

Exemple de schéma à l’ordre 4

Bien entendu, cette technique se poursuit pour
les ordres supérieurs. Un ordre souvent utilisé
4
avec cette méthode est l’ordre 4. Le schéma se pré-
sente alors comme suit :
2 1
Φ(x , t , h) = (k 1 + 2k 2 + 2k 3 + k 4 )
6
k 1 = f (x , t )
−4 −3 −2 −1 0 1  
h h
k 2 = f x + k 1, t + (4.50)
2 2
 
−2 h h
k 3 = f x + k 2, t +
2 2
k 4 = f (x + hk 3 , t + h)
−4
Bien entendu, des méthodes de prédiction-
fig. 4.6 - Domaine de stabilité pour correction basées sur la méthode de Runge-Kutta
Runge-Kutta (ordre 1 à 8)
peuvent aussi être mises en place.
Cette méthode requière beaucoup d’évalua-
tions de f , ce qui peut la pénaliser au niveau du coût. Concernant les domaines de stabilité,
on peut observer la même tendance que pour les méthodes d’Adams, voir figure 4.6.

4.1.10 Autres méthodes

Citons pour mémoire les méthodes de Hanning, de Milne-Simpson, les méthodes de tir
(shooting), les méthodes de collocation optimale, et une méthode bien adaptée à l’intégra-
tion des relations de comportements décrits par variables internes à 1 seul seuil, la méthode
de retour radial [57].
Il existe aussi une famille de méthodes appelée Galerkin discontinu en temps, plus ins-
pirée des formulations variationnelles en espace. Pour l’illustrer, considérons à nouveau le
problème de Cauchy (4.1). Il peut être réécrit en utilisant une fonction test x ⋆ à choisir dans
un espace adéquat :
Z T 
dx
− f (x , t ) x ⋆ dt = 0 (4.51)
0
dt
 4.2 Systèmes différentiels du premier ordre 53

L’idée est de pouvoir utiliser des discrétisations en temps qui tolèrent a priori des discon-
tinuités aux pas de temps : à t = t i , on peut définir une valeur à gauche x i− et une valeur
à droite x i+ différentes. Le saut est noté ¹x ºi = x i+ − x i− . L’intégrale précédente est donc à
prendre en un sens généralisé, qui fait « travailler » les sauts 4 :
Z T  Z t   
dx X i +1
dx
− f (x , t ) x ⋆ dt =  − f (x , t ) x ⋆ dt + ¹x ºi x i⋆  = 0 (4.52)
0
dt i t
dt
i

La fonction test est choisie dans le même espace que x , et on définit sa valeur en t i comme
x i⋆ = αx i⋆+ +(1−α)x i⋆− . Dans le cas standard, on prend α = 1 5 . Le problème peut alors s’écrire
sur chaque intervalle de temps :
Z t i +1  
dx
− f (x , t ) x ⋆ dt + ¹x ºi x i⋆ = 0 (4.53)
t
dt
i

Cette méthode est très ressemblante à une méthode de Runge-Kutta à un pas implicite (A-
stable) d’ordre 2l + 1 si on prend les espaces polynômiaux de degré inférieur ou égal à l par
incrément de temps (Lesaint et Raviart 1974). Il faut aussi prévoir une méthode numérique
de calcul d’intégrale en temps qui intègre exactement des polynômes de degré 2l (Cockburn
et Shu 1989).

4.2 Systèmes différentiels du premier ordre

Cette situation est rencontrée en éléments finis, lors de l’intégration de relations d’un com-
portement décrit par variables internes et avec un système de taille bien plus grande, lors de
calculs de diffusion, par exemple, un calcul de conduction thermique transitoire, qui amène
à résoudre un problème de la forme :

Mẋ + Kx = g(t )
(4.54)
x(0) = u0

où K, généralement symétrique positive, est la matrice de conduction, et M, généralement

symétrique, définie positive, est la matrice de capacité. x est le vecteur des températures
inconnues aux nœuds de la discrétisation.

4.2.1 Emploi des méthodes précédentes

Les méthodes précédentes se généralisent sans grande difficulté à ce cas de figure.

Exemple
Prenons le schéma d’Euler explicite. Le problème précédent peut s’écrire :

ẋ = −M−1 Kx + M−1 g (4.55)

On applique alors le schéma voulu :

xi +1 = xi + h(−M−1 Kxi + M−1 gi ) (4.56)

4. On suppose cependant f (x , t ) suffisamment régulière.
5. one-sided upwinding
 54 RÉSOLUTION DE SYSTÈMES DIFFÉRENTIELS

Pour éviter de faire intervenir directement M−1 , on préfère le mettre sous la forme :

Mxi +1 = Mxi + h(Kxi + gi )

(4.57)
x0 = u0

Stabilité
On emploie la même démarche que dans le cas des équations scalaires, dans la base propre
(pour les équations linéaires), avec la valeur propre la plus contraignante, c’est-à-dire la pul-
sation ω la plus grande.
2
Dans le cas d’Euler explicite, on a donc le critère de stabilité h ¶ ω2M
et on peut borner la
e
pulsation propre maxi du système par la pulsation propre maxi élémentaire ωM ¶ ωmax , voir
annexe B.
Pour avoir un vrai schéma explicite, il faut ne pas avoir à résoudre de système linéaire
couplé dans (4.57). Il faut donc avoir une matrice de masse diagonale. Une modification du
problème (une approximation supplémentaire dans la modélisation) consiste à remplacer
la matrice de masse par une matrice de masse « lumpée »diagonale. Pour cela, plusieurs
techniques existent, comme une matrice proportionnelle à la masse diagonale, avec respect
de la masse totale par élément, ou même la somme des lignes de la matrice de départ.
Outre une approximation supplémentaire [58], on peut se poser la question de la conser-
vation des propriétés du schéma. . .

4.2.2 Choix d’un schéma : influence de l’ordre d’intégration

Si on considère que le coût est directement lié
106 au nombre d’évaluations de la fonction f , on
Euler explicite
peut tracer le coût en fonction de l’erreur ob-
tenue pour différents schémas.
Heun
104 Par exemple, si on choisit le cas test suivant :
! !
ẋ y
Runge-Kutta 4 =
ẏ (1 − x 2 )y − x
! ! (4.58)
102
10−8 10−6 10−4 10−2 x (0) u0
=
y (0) 0
fig. 4.7 - Influence de l’ordre d’un
schéma sur le coût Il s’agit d’un système différentiel du premier
ordre, non linéaire, pour t ∈ [0, T]. Si α ≈
2,008619861986087484313650940188 et avec T ≈ 6,6632868593231301896996820305, la so-
lution est périodique, de période T.
Une mesure de l’erreur peut donc être de prendre à la fin de l’intervalle d’étude :
! !
xT x0
e =k − k (4.59)
yT y0

La figure 4.7 reporte en échelles logarithmiques le nombre d’évaluations de f en fonction du

niveau d’erreur e pour plusieurs schémas d’ordres différents : Euler explicite d’ordre 1, Heun
 4.3 Systèmes différentiels du deuxième ordre 55

d’ordre 2 et Runge-Kutta d’ordre 4. On voit donc l’intérêt de mettre en œuvre des schémas
d’ordre élevé, en gardant une stabilité raisonnable.

4.3 Systèmes différentiels du deuxième ordre

On trouve ce genre de problème par exemple en dynamique transitoire des structures :

Mẍ + Cẋ + Kx = g(t )

ẋ(0) = v0 (4.60)

x(0) = u0

M est la matrice de masse, C est la matrice d’amortissement et K la matrice de rigidité.

De façon générale, l’ordre d’un système linéaire peut être abaissé, ici à l’ordre 1, mais
en doublant la taille du système. Ici, on s’intéressera aux traitement direct du système au
second ordre, dans le cas où on n’a pas d’amortissement, c’est-à-dire où C = 0. L’équilibre
est donc :

Mẍ + Kx = g(t ) (4.61)

4.3.1 Méthode de Newmark

Il s’agit d’une famille de méthodes implicites à un pas [59] et définie en 1959. Elle se présente
sous la forme suivante :
1

xi +1 = xi + h ẋi + h 2 (1 − 2β)ẍi + 2βẍi +1
2 (4.62)


ẋi +1 = ẋi + h (1 − γ)ẍi + γẍi +1

On retrouve bien le caractère implicite par l’apparition de quantités indicées i +1 au second

membre. On peut aussi remarquer que si β = γ = 0, on retrouve des développements de
Taylor puisque le schéma de Newmark ressemble à un développement au troisième ordre :
h2 h 3 ẍi +1 − ẍi
xi +1 = xi + h ẋi + ẍi + 6β
2 6 h (4.63)
h 2 ẍi +1 − ẍi
ẋi +1 = ẋi + h ẍi + 2γ
2 h
Il reste à écrire l’équilibre (4.61) à l’instant t i +1 avec le schéma précédent. On peut alors
prendre l’accélération comme inconnue :
   
2
2 1
βh K + M ẍi +1 = gi +1 − K xi + h ẋi + h − β ẍi (4.64)
2
et on obtient le système à résoudre pour atteindre ẍi +1 ; mais on peut aussi l’écrire avec le
déplacement comme inconnue :
   
2
2 2 1
βh K + M xi +1 = βh gi +1 + M xi + h ẋi + h − β ẍi (4.65)
2
qui est le système à résoudre pour atteindre xi +1 . Dans ce dernier cas, cependant, un pro-
blème subsiste si β = 0. En effet, soit on utilise l’équilibre pour revenir à l’accélération :

ẍi +1 = M−1 (gi +1 − Kxi +1 ) (4.66)

56 RÉSOLUTION DE SYSTÈMES DIFFÉRENTIELS

ce qui peut être coûteux ; soit on utilise la relation inverse de (4.62) :

1 1 1 − 2β
ẍi +1 = 2
(xi +1 − xi ) − ẋi − ẍi (4.67)
βh βh 2β

mais on s’aperçoit bien du problème quand β = 0.

Cas particuliers
− γ = 12 , β = 1
4
: méthode des accélérations moyennes ;

− γ = 12 , β = 1
6
: méthode des accélérations linéaires ;

− γ = 12 , β = 0 : méthode des différences centrées, qui peut être explicite 6 ;

− γ = 12 , β = 1
12
: méthode de Fox-Goodwin.

Ordres
1
On peut montrer qu’il s’agit d’une méthode d’ordre 1 pour γ 6= 2
et d’ordre 2 pour γ = 12 .

Stabilité
Dans le cas d’une seule équation scalaire, ẍ = −ω2 x , le schéma donne :
  !   !
1 + (hω)2 β 0 x i +1 1 − (hω)2 ( 21 − β) h xi
  =  (4.68)
hω2 α 1 ẋ i +1 −hω2 (1 − α) 1 ẋ i
| {z } | {z } | {z } | {z }
A1 yi +1 A2 yi

d’où la relation yi +1 = A−1 −1

1 A2 yi et la matrice d’itération B = A1 A2 . La condition de stabilité
est donc ̺(B) < 1 ou bien ̺(B) ¶ 1 avec des valeurs propres simples. L’équation caractéris-
tique permettant de déterminer les valeurs propres de B est λ2 − 2b λ + c = 0 avec :
   
1 (hω)2 1 (hω)2
2b = 2 − α + c =1− α− (4.69)
2 1 + β(hω)2 2 1 + β(hω)2

les conditions de stabilité sont donc :

c ¶1 1 − 2b + c ¾ 0 1 + 2b + c ¾ 0 (4.70)

sauf en ce qui concerne les points isolés (b, c ) = (−1, 1) et (1, 1). D’où les conditions en terme
de coefficients de la méthode :

1 hω(2β − γ)
γ¾ 2+ ¾0 (4.71)
2 1 + β(hω)2

En conclusion, si 2β ¾ γ ¾ 12 , le schéma est inconditionnellement stable (IS) ; si γ ¾ 1

2
et
γ > 2β, le schéma est conditionnellement stable (CS) et la condition de stabilité est :

1
h¶ p (4.72)
ω γ/2 − β

La figure 4.8 illustre ces résultats.

6. Cette méthode est aussi appelée Velocity Verlet algorithm ou Leap Frog en dynamique moléculaire !
4.3 Systèmes différentiels du deuxième ordre 57

c β

1 IS
racines
complexes
1
conjuguées
β γ

−1 1 0 1/2
racines
réelles

−1

fig. 4.8 - Stabilité des schémas de Newmark

Influence de l’amortissement et comportement hautes fréquences

Considérons le cas où l’amortissement n’est pas négligé dans l’équation scalaire 7 :

ẍ + 2ζωẋ + ω2 x = 0 (4.73)

On peut montrer dans ce cas que la condition de stabilité s’écrit :

€ Š € Š1/2
ζ γ − 12 + 12 γ − β + ζ2 (γ − 12 )2 1
h̄ ¶ 1
(4.74)
γ−β ω
2

Dans le cas où γ > 12 , l’amortissement a un effet stabilisant puisque sur les conditions de
stabilité, on a h̄ < h.
Concernant les modes à haute fréquence, ils ne sont pas représentatifs du modèle continu
mais sont néanmoins sollicités lors de l’intégration et conduisent à la présence de parasites
non désirés. On utilise classiquement une dissipation numérique pour amortir ces hautes
fréquences.
En choisissant γ > 12 , on prend β pour maximiser l’amortissement (région des racines
complexes conjuguées sur la figure 4.8), ce qui conduit à :
 
1 1 2
β¶ γ+ (4.75)
4 2

Le point de dissipation Haute Fréquence maximal est déterminé par le couple (β = 1, γ =

3
2
), mais il présente le gros inconvénient d’être peu précis en Basse Fréquence. Un choix
généralement recommandé pour les méthodes de Newmark est alors :
 
1 1 1 2
γ> β= γ+ (4.76)
2 4 2

La figure 4.9 présente l’utilisation des schémas de Newmark sans amortissement pour le cas
de réflexions d’ondes dans une poutre unidimensionnelle, ce qui met en évidence la pré-
sence de perturbations à hautes fréquences. La figure 4.10 présente quant à elle l’influence
de l’amortissement numérique pour ces mêmes schémas.
7. Cette écriture dans la base propre d’un système entier correspond à un modèle d’amortissement de Rayleigh.
58 RÉSOLUTION DE SYSTÈMES DIFFÉRENTIELS

contraintes (MPa) contraintes (MPa)

200 200

100 100

0 0

0 200 400 600 800 0 200 400 600 800

temps temps
(a) différences centrées (b) Fox-Goodwin

contraintes (MPa) contraintes (MPa)

200 200

100 100

0 0

0 200 400 600 800 0 200 400 600 800

temps temps
(c) accélération linéaire (d) accélération moyenne

fig. 4.9 - Parasites hautes fréquences (d’après [60])

Toujours pour illustrer l’amortissement numérique, considérons un bilan d’énergie entre

les instants t i et t i +1 [61]. L’écriture du schéma de Newmark (4.62) en différence entre ces
deux instants, en utilisant la notation suivante :
1
〈xi 〉 = (xi +1 + xi )
2 (4.77)
[xi ] = xi +1 − xi

conduit à :
h2
[xi ] = h 〈ẋi 〉 + (2β − γ)[ẍi ]
2  (4.78)
1
[ẋi ] = h 〈ẍi 〉 + h γ − [ẍi ]
2
De même, la différence des deux équilibres (4.61) devient :

M[ẍi ] + K[ẋi ] − [gi ] = 0 (4.79)

En multipliant cette dernière relation par [ẋi ]T et en utilisant la propriété 〈xi 〉T [xi ] = [ 12 x2i ],
avec des matrices M et K positives, on obtient le bilan :
     
1 T ⋆ 1 T 1 T 1
ẍ M ẍi + ẋ Kẋi = [ẋi ] [gi ] − γ − [ẍi ]T M⋆ [ẍi ] (4.80)
2 i 2 i h 2

où M⋆ = M+ 12 h 2 (2β−γ)K. On trouve donc, dans l’ordre, un terme lié à la variation de l’éner-

gie cinétique, un terme lié à la variation de l’énergie interne, un terme lié au travail des
 4.3 Systèmes différentiels du deuxième ordre 59

contraintes (MPa) contraintes (MPa)

200
200

100 100

0
0
0 200 400 600 800 0 200 400 600 800
temps temps
(a) g = 0,55 (b) g = 0,6

contraintes (MPa) contraintes (MPa)

200 200

100 100

0 0
0 200 400 600 800 0 200 400 600 800
temps temps
(c) g = 0,75 (d) g = 1,5

fig. 4.10 - Influence de l’amortissement numérique (d’après [60])

pertes énergétiques (%)

0
g = 0,55

g = 0,6

g = 0,75
−0,05

g = 1,5

−0,1
0 200 400 600 800
temps

fig. 4.11 - Influence de l’amortissement numérique (suite) (d’après [60])

efforts extérieurs et un dernier terme lié à l’amortissement numérique (dissipation). Cela

permet de conclure que si γ ¾ 0,5, alors la solution est d’amplitude bornée et que si γ = 0,5,
le schéma numérique ne dissipe pas d’énergie.

4.3.2 Méthode de Gear implicite à deux pas

Cette méthode [62] proposé en 1969, appelée aussi méthode BDF pour Backward Differen-
tiation Formulas, est réputée être très stable. Pour information, elle utilise le développement
 60 RÉSOLUTION DE SYSTÈMES DIFFÉRENTIELS

15
3

2 10

1 5
k =1 k =2 k =3 k =4 k =5 k =6
ℑ(z )

ℑ(z )
0 0

−1 −5

−2 −10
−3
−15

−1 0 1 2 3 4 5 6 7 −10 −5 0 5 10 15 20 25 30
ℜ(z ) ℜ(z )

fig. 4.12 - A-stabilité dans le plan complexe pour la méthode de Gear : k est le nombre de pas

suivant :
1
ẋi +1 = (3xi +1 − 4xi + xi −1 )
2h (4.81)
1
ẍi +1 = (3ẋi +1 − 4ẋi + ẋi −1 )
2h
De la même façon, l’équilibre (4.61) à l’instant t i +1 donne le système à résoudre pour at-
teindre ẍi +1 :
   
4 2 4 1 8 2
M + h K ẍi +1 = gi +1 − K xi − xi −1 + h ẋi − h ẋi −1 (4.82)
9 3 3 9 9

La figure 4.12 présente le domaine de stabilité de cette méthode.

4.3.3 θ -méthode
Historiquement attribuée à Wilson en 1968 [63], cette méthode estime une accélération li-
néaire entre les instants t i et t i +θ = t i + θh, θ ¾ 0 étant un paramètre de la méthode. On a
donc :
ẍi +θ − ẍi
ẍ(t i + τ) = τ + ẍi (4.83)
θh
et par intégration, on a les expressions :

1
ẋi +θ = ẋi + θh (ẍi + ẍi +θ )
2   (4.84)
1 1
xi +θ = xi + θh ẋi + θ 2 h 2 ẍi + ẍi +θ
3 6

Cette fois-ci, l’équilibre (4.61) à l’instant t i +θ conduit au système en xi +θ :

   
6 6 6
K + 2 2 M xi +θ = gi +θ + M xi + ẋi + 2ẍi (4.85)
θ h θ2h 2 θh

− θ ¾ 1,37 : la méthode est inconditionnellement stable ;

− θ = 1 : méthode des accélérations linéaires.
 4.4 Intégration de modèles de comportement 61

4.3.4 α-HHT
Cette méthode [64] proposée en 1977, ressemble a priori à la méthode de Newmark, à ceci
près qu’elle introduit une modification de l’équilibre (et donc en particulier du second membre)
en utilisant un paramètre additionnel α.
L’écriture du déplacement, comme celle de la vitesse sont identique à la méthode de
Newmark, et l’équilibre est écrit sous la forme :




Mẍi +1 + (1 + α)Cẋi +1 − αCẋi + (1 + α)Kxi +1 − αKxi = (1 + α)fi +1 − αfi (4.86)

Le démarrage est réalisé avec x0 et ẋ0 données ainsi que ẍ0 = M−1 (f0 − Cẋ0 − Kx0 ).
Pour le choix de paramètres de la forme β = 14 (1 − α)2 et γ = 1
2
− α, le schéma est du
second ordre pour les basses fréquences, il est inconditionnellement stable pour − 13 ¶ α ¶ 0,
la dissipation en hautes fréquences est maximum pour α = − 13 , enfin, si on impose de plus
α = 0 on retrouve la méthode des trapèzes.

4.3.5 Autres méthodes

Il existe un grand nombre de méthodes d’intégration plus ou moins adaptées à un type de
problème. En dynamique, citons encore la méthode de Park à trois pas [65].

4.4 Intégration de modèles de comportement

Si on considère le problème du calcul de structure non-linéaire de la section 3.1.1, on était
conduit dans les méthodes de résolution à estimer le résidu r(u) = K(u)u − f qui dépend de
l’état courant au travers du modèle de comportement du matériau :
Z
K(u)u = Trσ ǫ(u) dΩ (4.87)
Ω

puisque la contrainte σ dépend de l’état courant.

4.4.1 Exemple de modèle de plasticité

Utilisons les notations de la section 3.1.1. Les lois d’état sont de la forme :

σ = D(ǫ − ǫp )
(4.88)
Y = G(X)

où D est l’opérateur de Hooke et G la fonction d’écrouissage, caractéristiques du matériau.

Dans le cadre d’un modèle associé, le potentiel de dissipation est pris égal au seuil de plas-
ticité, sur lequel s’annule la fonction caractéristique f dépendant de l’état f (σ, Y).
Les lois d’évolution s’écrivent avec l’utilisation d’un multiplicateur plastique λ̇ à cause
du caractère non différentiable du potentiel de dissipation :
∂f
ǫ̇p = λ̇ (4.89)
∂σ
∂f
−Ẋ = λ̇ (4.90)
∂Y
 62 RÉSOLUTION DE SYSTÈMES DIFFÉRENTIELS

On a alors les conditions de complémentarité de Khun-Tucker :

f ¶ 0, λ̇ ¾ 0, f λ̇ = 0 (4.91)

L’expression du multiplicateur est donnée dans le cas de l’écoulement plastique par la condi-
tion de cohérence :

∂f ∂f
f˙ = 0 = σ̇ + Ẏ (4.92)
∂σ ∂Y

(dans le cas contraire λ̇ = 0 et le comportement est localement élastique). En remplaçant

σ̇ et Ẏ par leurs expressions issues de la dérivation des lois d’état, et en utilisant les lois
d’évolution, on obtient l’expression du multiplicateur :

∂f
∂σ
Dǫ̇
λ̇ = ∂f ∂f ∂f
(4.93)
G ∂f
D + ∂ Y ∂∂ X
∂σ ∂σ ∂Y

On peut donc formellement écrire (à partir des lois d’évolution, de l’expression précédente
du multiplicateur et en remplaçant chaque fois que nécessaire σ et Y par leurs expressions
en ǫ, ǫp et X issues des lois d’état) l’évolution du comportement sous la forme d’un système
d’équations non-linéaires différentielles ordinaires (ODE) :
! ! !
ǫ̇p ǫp
=g , ǫ, ǫ̇ (4.94)
−Ẋ X

Il est à noter que le temps n’intervient pas directement dans l’évolution, contrairement au
cas des comportements visqueux. Remarquons aussi que ce système ne peut être résolu :
il comporte trop d’inconnues par rapport au nombre d’équations disponibles. En effet, il
faut une sollicitation pour que le comportement évolue. En utilisation dans un code aux
éléments finis en déplacement, il est agréable de conserver une déformation admissible tout
le temps. Le modèle de comportement est donc généralement piloté en intégration par le
fait que l’évolution de la déformation totale est imposée. Dans le système précédent, ǫ et ǫ̇
sont alors des données.
La résolution de ce système différentiel peut être effectuée avec les techniques décrites
dans le chapitre 4. On peut ainsi utiliser des méthodes explicites (avec des sous-pas plus
nombreux) ou implicites (avec parfois un seul sous-pas par pas de chargement de la struc-
ture, c’est-à-dire par résolution de système différentiel). Classiquement, Runge-Kutta d’ordre
4 avec contrôle de la taille du sous-pas et Euler implicite sont utilisés. En fait, ces méthodes
sont utilisées sur une partie du problème à résoudre seulement, car souvent, des méthodes
de type retour radial sont employées.

4.4.2 Méthodes de type retour radial

Ces méthodes sont employées car elles sont mieux adaptées au type de problème à traiter.
On souhaite en effet retourner assez précisément sur le seuil après intégration du système
différentiel. De plus, il faut pouvoir prendre en compte le cas de la décharge élastique. Pour
4.4 Intégration de modèles de comportement 63

cela, on prend en compte la propriété de séparation en partie élastique et plastique des mo-
dèles de comportement, pour séparer en deux temps l’intégration. Considérons un incré-
ment de charge sur la structure, et un pilotage en incrément de déformation totale, noté ∆ǫ.
Cette méthode est aussi considérée comme une méthode de partage d’opérateur (operator
split, cf. prédiction élastique et correction plastique par un retour sur le seuil ou return map).
La notation δ représente une variation sur un éventuel sous-pas d’intégration. La technique
employée consiste en :
prédiction élastique : elle est telle que ∆ǫe + ∆ǫp = ∆ǫ, avec ∆ǫp = 0 et ∆X = 0, les lois
d’état (4.88) donnent sur cet incrément :
!   !
∆σ D 0 ∆ǫ
=  (4.95)
∆Y 0 0 0
Si, au bout de cet incrément, le nouvel état (σ + ∆σ, ǫp , X, Y + ∆Y) vérifie f ¶ 0, on est
à l’intérieur du seuil et le comportement est resté élastique. Sinon, il faut corriger la
prédiction de la façon suivante.
correction : à partir de l’état déterminé précédemment, on procède à une nouvelle correc-
tion (toujours notée avec ∆) ∆ǫe +∆ǫp = 0 et on cherche à vérifier les équations d’état,
les équations d’évolution et, généralement, directement le retour sur le seuil f = 0 à la
fin de la correction.
Cette phase de correction peut être réalisée par utilisation des méthodes d’intégration
déjà citées. On peut aussi utiliser des méthodes plus dédiées comme celle du retour
radial d’Ortiz et Simo [66, 57].

Retour radial d’Ortiz et Simo

Pour résoudre la phase de correction, l’idée est d’utiliser une méthode de Newton sur le
seuil. On cherche à avoir :

f (σ + ∆σ, Y + ∆Y) = 0 (4.96)

on procède par sous-pas (notés avec un δ) en utilisant la matrice jacobienne des dérivées
partielles de f lors des développements successifs :
∂f ∂f
f (σ + δσ, Y + δY) ≈ f (σ, Y) + δσ + δY (4.97)
∂σ ∂Y
pour lesquels on cherche à rester sur le seuil :
∂f ∂f
f+ δσ + δY = 0 (4.98)
∂σ ∂Y
Comme δǫ = 0 lors de la correction, les équations (4.89) et (4.90) donnent alors la valeur
d’un « multiplicateur plastique »pour le sous-pas :
f
δλ = ∂f ∂f ∂ f ∂G ∂ f
(4.99)
D
∂σ ∂σ
+ ∂Y ∂X ∂Y

Le potentiel de dissipation étant convexe, le dénominateur de l’expression précédente est

toujours positif. Après avoir sommé ce sous-incrément à l’état courant, on itère sur la cor-
rection jusqu’à annuler le résidu f (σ, Y).
 64 RÉSOLUTION DE SYSTÈMES DIFFÉRENTIELS

Rigidité tangente cohérente

Il a été constaté que si l’on utilise la matrice de rigidité tangente continue précédente, pour
itérer avec une méthode de Newton-Raphson (voir le chapitre 3.1.5), on perd la convergence
quadratique de la méthode, et que pour la retrouver, il faut utiliser la matrice de rigidité
tangente cohérente avec le schéma d’intégration employé.
Ce dernier peut être vu comme un outil pour produire, à partir d’un état donné :

s = (σ, Y, ǫp , X) (4.100)

un nouvel état satisfaisant au comportement du matériau, piloté par exemple par un incré-
ment de déformation totale ∆ǫ. Ce nouvel état est noté :

s + ∆s = (σ + ∆σ, Y + ∆Y, ǫp + ∆ǫp , X + ∆X) (4.101)

Si la rigidité tangente continue à l’état s est DT telle que σ̇ = DT ǫ̇, la matrice tangente
cohérente DC est définie par δ∆σ = DC δ∆ǫ. Elle relie la variation de correction de contrainte
obtenue, à la variation de pilotage en incrément de déformation totale : piloter avec ∆ǫ+δ∆ǫ
à partir du même état s conduira à un nouvel état suivant s + δs + δ∆s .
On peut obtenir DC par une méthode de perturbation (on procède à plusieurs intégra-
tions avec des perturbations différentes δ∆ǫ sur le pilotage) mais cela est assez coûteux nu-
mériquement, et nécessite le choix de la valeur de la perturbation (trop petite, elle conduit
à des problèmes de précision numériques, trop grande, à une mauvaise approximation de
la dérivée). Une autre technique consiste à dériver analytiquement le schéma d’intégration,
ce qui peut s’avérer assez pénible suivant le schéma et le modèle de comportement.

4.4.3 Exemple de modèle de viscoplasticité

Par rapport au modèle précédente, un grand nombre de modèles de viscoplasticité s’écrivent
en n’imposant plus de signe sur la valeur de la fonction seuil f (σ, Y), mais avec un terme de
rappel sur le multiplicateur viscoplastique λ̇ : les conditions de complémentarité de Kuhn-
Tucker (4.91) sont alors remplacées par :
‚  Œn
f (σ, Y)
λ̇ = (4.102)
m +

où les parties positive et négative sont respectivement notées :

1 1
〈x 〉+ = max(0,x ) = (x + |x |) et 〈x 〉− = min(0,x ) = (x − |x |) = − 〈−x 〉+ (4.103)
2 2

m est un module de viscoplasticité et n un exposant de viscosité.

Il est possible de définir l’équivalent du seuil de plasticité de la façon suivante : si on


considère f¯ = f − f , alors f¯ = 0 correspond à f ¾ 0 c’est-à-dire à un écoulement visco-
+
plastique, et f¯ < 0 correspond à f ¶ 0 c’est-à-dire à une évolution élastique, et dans ce cas
f¯ = f .
4.4 Intégration de modèles de comportement 65

On peut aussi se ramener à une re-

Ayant les incréments donnés ∆ǫ et ∆t ,
lation de complémentarité : si on note prédiction élastique avec (4.104)
prédiction du seuil f (σ, Y)
µ = λ̇, on doit en effet avoir µ ¾ 0, f˜ =
si f > 0 alors
f − m µ1/n ¶ 0 et µ f˜ = 0. La principale correction viscoplastique
µ = 0, f˜ = f
différence par rapport au cas de la plasti- tant que f˜ > 0
cité est la dépendance avec le temps au estimer ǫ̇p et Ẋ avec (4.108)
calculer la correction δµ avec (4.109)
travers des vitesses des quantités cher- mise à jour de σ, Y, ǫp , X avec (4.107)
chées. L’intégration de la loi de compor- mise à jour du multiplicateur µ ← µ + δµ
mise à jour du seuil corrigé f˜
tement ne sera plus pilotée par le seul in- fin
crément de déformation ∆ǫ, mais aussi fin

par l’incrément de temps ∆t . Le principe

alg. 4.13 - Intégration du comportement
consiste encore à prédire une évolution viscoplastique
élastique :

∆σ = D∆ǫ, ∆Y = 0
(4.104)
∆X = 0, ∆ǫp = 0

et à estimer la valeur du seuil f (σ, Y). Si celle-ci est négative, la solution est effectivement
élastique. Sinon, il faut procéder à une correction viscoplastique ; cette étape est encore ité-
rative, mais il faut bien réévaluer toutes les quantités (en particulier les dérivées par l’inter-
médiaire du schéma d’intégration en temps choisi, ici la θ-méthode) par rapport au précé-
dent pas de temps t i (elle sont ici indicées par i ). La correction d’effectue toujours à défor-
mation totale nulle, et l’estimation du seuil corrigé est :

∂ f˜ ∂ f˜ ∂ f˜
f˜(σ + δσ, Y + δY, µ + δµ) ≈ f˜(σ, Y, µ) + δσ + δY + δµ (4.105)
∂σ ∂Y ∂µ
pour lesquels on cherche à rester sur le seuil corrigé :
∂f ∂f m
f (σ, Y) − m µ1/n + δσ + δY − µ1/n−1 δµ = 0 (4.106)
∂σ ∂Y n
∂G
avec δσ = −Dδǫp , δY = ∂X
δX et les incréments :

δǫp ∂f
− ǫ̇p
= (µ + δµ)
θ∆t ∂σ (4.107)
δX ∂f
= −(µ + δµ) − Ẋ
θ∆t ∂Y
où les dérivées sont estimées par rapport au précédent pas de temps (et doivent être rééva-
luée à chaque itéré de la correction) :
1−θ
ǫp − ǫp i
ǫ̇p = ǫ̇p i
−
θ∆t θ
(4.108)
X − Xi 1 − θ
Ẋ = − Ẋi
θ∆t θ
L’équation (4.105) permet alors d’obtenir l’incrément de multiplicateur viscoplastique :
f −m µ1/n ∂f ∂f G
θ∆t
+ ∂ σ Dǫ̇p − ∂ Y ∂∂ X Ẋ
δµ = ∂f ∂f ∂ f ∂G ∂ f m 1/n−1 1
(4.109)
D + ∂Y ∂X ∂Y + n µ
∂σ ∂σ θ∆t
 66 RÉSOLUTION DE SYSTÈMES DIFFÉRENTIELS

La valeur courante du multiplicateur doit alors aussi être mise à jour par µ + δµ, à chaque
itéré de la correction, qui est poursuivie jusqu’à retourner sur le seuil corrigé f˜.
Connaissant la solution à t i , on procède à une succession de prédictions élastiques jus-
qu’à rééquilibrer les forces extérieures (c’est-à-dire à annuler le résidu des équations d’équi-
libre) et d’intégration du comportement selon l’algorithme décrit dans 4.13. La solution à
t i +1 est celle qui doit vérifier à la fois l’équilibre et le comportement.

4.5 Cas de la dynamique non régulière

Dans les approches précédentes, on a souvent considéré que les évolutions temporelles
étaient suffisamment régulières pour pouvoir les approcher avec un développement limité.
Ce n’est pas forcément toujours le cas ; par exemple dans le cas de la fissuration dynamique,
lorsque les structures sont soumises à des chocs ou impacts. Dans ce dernier cas, en parti-
culier quand les corps sont des solides rigides, les vitesses ne sont plus dérivables (lors d’un
choc, un corps rigide change peu en position, mais beaucoup en vitesse, par exemple dans
le cas d’un rebond). Une technique (dite event driven) consiste à traiter de façon standard
les phases de vol libre c’est-à-dire sans impact, à repérer les instants de choc, et à traiter
ceux-ci de façon particulière, avec un saut de vitesse de part et d’autre de l’instant de choc.
Avec une modélisation en solides rigides, il manque une équation de comportement de
l’impact pour fermer le problème. Un tel comportement peut être modélisé par la loi de
Newton 8 . Elle postule que la vitesse normale du point de contact — d’où des difficultés
pour modéliser un impact sur une surface entière — après impact (notée v n+ ) est opposée
en direction à celle avant impact (notée v n− ) :

v n+ = −e v n− (4.110)

le saut de vitesse lors de l’impact étant v n+ − v n− . Le coefficient matériau e , appelé coefficient

de restitution d’énergie, traduit le comportement de l’impact ; 0 ¶ e ¶ 1. Le problème est
qu’il dépend de beaucoup de paramètres du problème, puisqu’il représente de façon phé-
noménologique le comportement complexe des propagations / réflections d’ondes lors de
l’impact, ou des impacts multiples. . .
Sans présence de frottement, il n’est pas nécessaire d’écrire autre choses sur les glisse-
ments tangentiels lors de l’impact. On se contentera de décrire ce type de problème ici.
Une autre situation particulière est celle des chocs nombreux dans les matériaux gra-
nulaires ou soumis à la fragmentation. . . Dans ce cas, une technique event driven est in-
applicable à cause des impacts qui peuvent être extrêmement nombreux. Une alternative
consiste à employer une technique de time-stepping. Le pas de temps est fixé et on s’em-
ploie à déterminer les vitesses immédiatement après les piquets de temps, en considérant
l’impact moyen pendant le pas de temps courant. On ne s’intéresse alors pas vraiment à
l’évolution du système continûment sur le pas de temps, mais plutôt à une configuration
8. Elle masque le fait que le solide est quand même légèrement déformable et que le rebond est dû à un échange entre
l’énergie cinétique et l’énergie de déformation, lors de la propagation des ondes de chocs dans le solide
4.5 Cas de la dynamique non régulière 67

possible en chaque piquet de temps.

Considérons tout d’abord un système dynamique régulier (sans choc), rigide, évoluant
sous l’action d’une force extérieure f . Son équation du mouvement est de la forme :

m ẍ = f (4.111)

où l’on suppose m constante. On se propose d’intégrer cette équation différentielle avec un

schéma particulier en diminuant l’ordre du système 9 . Pour cela, commençons par intégrer
l’équation du mouvement sur un pas de temps [t i , t i +1 ] ; on obtient :
Z t i +1
m (ẋ i +1 − ẋ i ) = f (t ) dt (4.112)
ti

L’intégrale de la force extérieure (régulière) peut être approximée par un schéma d’intégra-
tion, par exemple la méthode des trapèzes généralisés, voir section 4.1.6. Le problème s’écrit
alors :

m ẋ i +1 = m ẋ i + (1 − α)h f i + αh f i +1 (4.113)

La position est récupérée par une intégration similaire de la vitesse :

Z t i +1

x i +1 = x i + ẋ dt ≈ x i + (1 − α)h ẋ i + αh ẋ i +1 (4.114)
ti

Le lecteur pourra vérifier qu’il s’agit exactement de l’application de la méthode des trapèzes
généralisés au système du premier ordre :

m v̇ = f
(4.115)
ẋ = v

Cette façon de procéder va permettre de traiter les systèmes non réguliers comprenant des
chocs. En effet, l’équation de la dynamique entre deux instants immédiatement suivant les
piquets de temps t i et t i +1 va s’écrire :
Z t i +1

m (ẋ i +1 − ẋ i ) = f (t ) dt + R (4.116)
ti

où R est l’impulsion totale, c’est-à-dire, la somme des impulsions lors des impacts éventuels
sur le pas de temps traité.
Dans le cas d’une évolution statique, le contact unilatéral du système contre un mur
rigide placé en x = 0 par exemple ferait intervenir un effort F et non pas une impulsion R,
avec une loi de comportement du contact parfait du type 0 ¶ x ⊥ F ¾ 0, qui est une version
concise de x ¾ 0, F ¾ 0, x · F = 0. Le pendant de ce comportement en dynamique est la
condition de Signorini-vitesse issue du Lemme de Moreau :
– si x > 0, R = 0
– si x = 0, 0 ¶ ẋ ⊥ R ¾ 0
9. On travaillera alors à la fois sur les positions x i et les vitesses ẋ i aux piquets de temps t i .
68 RÉSOLUTION DE SYSTÈMES DIFFÉRENTIELS

Pour dicrétiser en temps ce comportement, on utilise en général une version explicite sur la
configuration — c’est-à-dire qu’on va utiliser un prédicteur pour la position x qui permettra
de sélectionner le cas correspondant, et qu’on notera x p :
– si x p > 0, Ri +1 = 0
– si x p ¶ 0, 0 ¶ ẋ i +1 ⊥ Ri +1 ¾ 0
Un prédicteur classique (saute-mouton) est x p = x i + h2 ẋ i . L’inconvénient est la possibilité
d’avoir une pénétration résiduelle, d’où l’inégalité utilisée sur les valeurs possiblement né-
gatives de x p , qui tend cependant vers 0 avec le pas de temps. Ce modèle de comportement
correspond à un impact plastique (sans rebond), qui dissipe le plus l’énergie. Dans le cas de
choc unique, on peut retrouver l’équivalent du modèle de Newton en substituant ẋ i +1 par
une vitesse dite formelle :

ẋ i +1 − ẋ i
ṽ = (4.117)
1+e
 RÉFÉRENCES

[1] B. Mohammadi and J.-H. Saiac. Pratique de la simulation numérique. Dunod /

Industries et Technologies, 2003. ISBN 210006407X. 1
[2] G.H. Golub and C.F. Van Loan. Matrix Computations. The Johns Hopkins University
Press, Baltimore, MD, USA, 2e edition, 1989. ISBN 0-8018-3772-3, 0-8018-3739-1. 2
[3] K. J. Bathe and E. L. Wilson. Numerical Analysis in Finite Element Analysis.
Prentice-Hall, 1982. 2
[4] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz,
A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK Users’ Guide.
Society for Industrial and Applied Mathematics, Philadelphia, PA, 3e edition, 1999.
ISBN 0-89871-447-8 (paperback). 2
[5] L.S. Blackford, J. Choi, A. Cleary, E. d’Azevedo, J. Demmel, I. Dhillon, J. Dongarra,
S. Hammarling, G. Henry, A. Petitet, K. Stanley, D. Walker, and R.C. Whaley.
ScaLAPACK Users’ Guide. Society for Industrial and Applied Mathematics,
Philadelphia, PA, 1997. ISBN 0-89871-397-8 (paperback). 3
[6] H. Oudin. Méthode des éléments finis. École Centrale de Nantes, France, 2008.
[Link] 5
[7] Y. Saad. Iterative Methods for Sparse Linear Sytems. 2000. 5
[8] R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
C. Romine, and H. Van der Vorst. Templates for the Solution of Linear Systems :
Building Blocks for Iterative Methods. SIAM, 1994.
[Link] 5
[9] J. J. Dongarra, H. W. Meuer, and E. Strohmeier. TOP 500 supercomputer sites. In
Proceedings of Supercomputer ’99 Conference, number 17, Mannheim, 1999.
Disponible sur [Link] et
[Link] 8
[10] H. W. Meuer, E. Strohmaier, J. J. Dongarra, and H. Simon. Top500 supercomputer
sites. In 16th International Supercomputer Conference : Applications — Architectures
— Trends, 2001. [Link]
[Link] 8
[11] E. Cuthill and J. McKee. Reducing the bandwidth of sparce symmetric matrices. In
Proc. of the 24th National Conference of the ACM, volume P-69, pages 157–172, 1969. 9
[12] A. George. Nested dissection of a regular finite-element mesh. SIAM Journal on
Numerical Analysis, 10 :345–363, 1973. doi:10.1137/0710032. 9
[13] W. F. Tinney and J. W. Walker. Direct solutions of sparse network equations by
optimally ordered triangular factorization. In Proc. of the IEEE, volume 55, pages
1801–1809, 1967. 9
[14] A. George and J. Liu. The evolution of the minimum degree ordering algorithm. SIAM
Review, 31(1) :1–19, 1989. 9
[15] I. S. Duff, A. M. Erisman, and J. K. Reid. Direct Methods for Sparse Matrices. Oxford
University Press, Oxford, 1990. 10
70 RÉFÉRENCES

[16] Comm. Benoît. Note sur une méthode de résolution des équations normales
provenant de l’application de la méthode des moindres carrés à un système
d’équations linéaires en nombre inférieur à celui des inconnues — application de la
méthode à la résolution d’un système d’équations linéaires (procédé du commandant
cholesky). Bulletin Géodésique, 2 :5–77, 1924. 12
[17] B. M. Irons. A frontal solution program for finite element analysis. International
Journal for Numerical Methods in Engineering, 2(1) :5–32, 1970.
doi:10.1002/nme.1620020104. 14
[18] I. S. Duff. Parallel implementation of multifrontal schemes. Parallel Computing, 3 :
192–204, 1986. 14
[19] C. Farhat and D. Rixen. Encyclopedia of Vibration, chapter Linear Algebra, pages
710–720. Academic Press, 2002. 17
[20] L. Hagemana and D. Young. Applied Iterative Methods. Academic Press, New York,
1981. 18
[21] R. Southwell. Relaxation Methods in Theoretical Physics. Clarendaon Press, Oxford,
1946. 19
[22] M. R. Hestenes and E. Stiefel. Methods of conjugate gradients for solving linear
systems. Journal of research of the National Bureau of Standards, 49 :409–436, 1952. 22
[23] G. Shortley. Use of Tschebyscheff-polynomial operators in the numerical solution of
boundary-value problems. Journal of Applied Physics, 24 :392–396, 1953.
doi:10.1063/1.1721292. 23
[24] C. C. Paige and M. A. Saunders. Solution of sparse indefinite systems of linear
equations. SIAM Journal on Numerical Analysis, 12 :617–629, 1975.
doi:10.1137/0712047. 23
[25] R. W. Freund and N. M. Nachtigal. QMR : A quasi-minimal residual method for
nonhermitian matrices. Numerische Mathematik, 60 :315–339, 1991. 23
[26] P. Sonneveld. CGS, a fast Lanczos-type solver for nonsymmetric linear systems. SIAM
Journal on Scientific and Statistical Computing, 10 :36–52, 1989. 23
[27] H. A. Van der Vorst. BiCGSTAB : A fast and smoothly converging variant of Bi-CG for
the solution of nonsymmetric linear systems. SIAM Journal on Scientific and
Statistical Computing, 13 :631–644, 1992. 23
[28] Y. Saad and M. H. Schultz. GMRES : A generalized minimal residual algorithm for
solving nonsymmetric linear systems. SIAM Journal on Scientific and Statistical
Computing, 7 :856–869, 1986. 23
[29] J. S. Rao. Advanced Theory of Vibration. John Wiley & Sons, 1989. 25
[30] M. Géradin and D. Rixen. Théorie des Vibrations. Application à la dynamique des
structures. Physique Fondamentale et Appliquée. Masson, 1993. 25
[31] M. Géradin and D. Rixen. Mechanical Vibrations. Theory and Application to Structural
Dynamics. Wiley, 1994. 25
[32] J. H. Wilkinson. The Algebraic Eigenvalue Problem. Oxford University Press, 1965.
ISBN 0-19-853403-5, 0-19-853418-3. 25
[33] Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst, editors. Templates for
the Solution of Algebraic Eigenvalue Problems : A Practical Guide. 1999. Draft. 25
[34] C. G. J. Jacobi. Ueber ein leichtes verfahren, die in der theorie der säcularstörungen
vorkommenden gleichungen numerish aufsulösen. Journal für die reine und
angewandle Mathematik, pages 51–94, 1846. 27
[35] R. A. Willoughby. Collection of articles honoring Alston S. Householder. ACM, 18 :
3–58, 1975. 31
 71

[36] A. S. Householder. Unitary triangularization of a nonsymmetric matrix. Journal of the

Association of Computing Machinery, 5 :339–342, 1958. 31
[37] G. Hessenberg. Transzendenz von e und ; ein Beitrag zur hoheren Mathematik vom
elementaren Standpunkte aus. von 1912. Johnson Reprint Corp., New York, 1965. 33
[38] C. Lanczos. An iteration method dor the solution of the eigenvalue problem of linear
differential and integral operators. Journal of research of the National Bureau of
Standards, 45 :225–280, 1950. 35
[39] C. T. Kelley. Solving Nonlinear Equations with Newton’s Method, volume 1 of
Fundamentals of Algorithms. Society for Industrial and Applied Mathematics,
Philadelphia, PA, 2003. ISBN 0-89871-546-6. 37
[40] M. A. Crisfield. Non-Linear Finite Element Analysis of Solids and Structures. John
Wiley, 1991. 37
[41] J. Lemaitre and J.-L. Chaboche. Mechanics of Solid Materials. Cambridge University
Press, 1994. 38
[42] R. Fletcher and C. Reeves. Function minimization by conjugate gradients. Comput. J.,
7 :149–154, 1964. 39
[43] E. Polak and G. Ribière. Note sur la convergence de directions conjuguées. Rev.
Française Informat. Recherche Opérationnelle, (16) :35–43, 1969. 39
[44] Y. Dai and Y. Yuan. Global convergence of the method of shortest residuals.
Numerische Mathematik, 83(4) :581–598, 1999. doi:10.1007/s002119900080. 39
[45] W. J. Duncan. Some devices for the solution of large sets of simultaneous linear
equations. Philos. Mag. Ser., 35(7) :660–670, 1944. 40
[46] W. W. Hager. Updating the inverse of a matrix. SIAM Review, 31(2) :221–239, 1989. 40
[47] C. Huang and G. Verchery. An exact structural static reanalysis method.
Communications in Numerical Methods in Engineering, 13 :103–112, 1997. 40
[48] R. Fletcher and M. J. D. Powell. A rapidly convergent descent method for
minimization. Computer Journal, 6 :163–168, 1963. 40
[49] C. G. Broyden. The convergence of a class of double-rank minimization algorithms.
Journal of the Institute of Mathematics and its Applications, 6 :76–90, 1970. 41
[50] L. Nazareth. A relationship between BFGS and conjugate gradient algorithms and its
implications for new algorithms. SIAM Journal on Numerical Analysis, 16 :794–800,
1979. doi:10.1137/0716059. 41
[51] T. G. Kolda, D. P. O’Leary, and L. Nazareth. BFGS with update skipping and varying
memory. SIAM Journal on Optimization, 8(4) :1060–1083, 1998. 41
[52] D.C. Liu and J. Nocedal. On the limited memory BFGS method for large scale
optimization. Numerische Mathematik, 45 :503–528, 1989. 41
[53] T. J. R. Hughes. The Finite Element Method, Linear Static and Dynamic Finite Element
Analysis. Prentice-Hall, Englewood Cliffs, 1987. 43
[54] D. W. Jordan and P. Smith. Nonlinear Ordinary Differential Equations. Oxford applied
mathematics and computing science series, 1987. 43
[55] G. Dahlquist. A special stability problem for linear multistep methods. BIT. Numerical
Mathematics, 3(1) :27–43, 1963. doi:10.1007/BF01963532. 46
[56] W. Kutta. Beitrag zur naherungsweisen integration totaler differentialgleichunge. Zeit.
Math. Physik, 46 :435–453, 1901. 51
[57] M. Ortiz and J. C. Simo. An analysis of a new class of integration algorithms for
elasto-plastic constitutive relations. International Journal for Numerical Methods in
Engineering, 23(3) :353–366, 1986. doi:10.1002/nme.1620230303. 52, 63
72 RÉFÉRENCES

[58] J.-P. Combe. Sur le contrôle des calculs en dynamique rapide. Application aux
problèmes d’impact. PhD thesis, ENS de Cachan, LMT-Cachan, 2000. 54
[59] N. M. Newmark. A method for computation of structural dynamics. In EMS, editor,
Proc. of the American Society of Civil Engineers, volume 85, pages 67–94, 1959. 55
[60] H. Lemoussu. Approche non incrémentale pour le calcul de choc avec contact
unilatéral avec frottement. PhD thesis, ENS de Cachan, LMT-Cachan, 1999. 58, 59
[61] T. J. R. Hughes and T. Belytschko. Nonlinear Finite Element Analysis. Zace sevices Ltd
- ICE Division, 1995. 58
[62] C. W. Gear. Numerical Initial Value Problems in Ordinary Differential Equation.
Prentice-Hall, 1971. 59
[63] K. J. Bathe and E. L Wilson. Stability and accuracy analysis of direct integration
method. Earthquake Eng. Struct. Dyn., 1 :283–291, 1973. 60
[64] H. M. Hilber, T. J. R. Hughes, and R. L. Taylor. Improved numerical dissipation for
time integration algorithm in structural dynamics. Earthquake Eng. Struct. Dyn., 5 :
283–292, 1977. 61
[65] K. C. Park. An improved stiffy stable method for direct integration of non linear
structural dynamics equations. Journal of Applied Mechanics, pages 164–470, 1975. 61
[66] M. Ortiz and E. P. Popov. Accuracy and stability of integration algorithms for
elastoplastic constitutive relations. International Journal for Numerical Methods in
Engineering, 21(9) :1561–1576, 1985. doi:10.1002/nme.1620210902. 63
[67] A. C. Aitken. On the iterative solution of a system of linear equations. Proceedings of
the Royal Society of Edinburgh, 63 :52–60, 1950. 77
 A
FACTORISATION
DE CHOLESKY -
ORTHOGONALISATION
D’UN SOUS-ESPACE

On considère un ensemble de m vecteurs de taille n vi indépendants, stockés dans les co-

lonnes d’une matrice rectangulaire V de taille (n, m ). Ils génèrent un sous-espace de l’es-
pace de départ. Dans un certain nombre d’applications, il est utile de construire à partir
de ces vecteurs, une nouvelle base orthogonale U. Pour ce faire, on utilise généralement le
procédé d’orthogonalisation de Gram-Schmidt. On cherche donc U engendrant le même
sous-espace que V et tel que UT U = I.
Pour ce faire, on peut utiliser une technique basée sur la factorisation précédente : si on
construit la matrice M = VT V, pleine et de taille réduite (n, n) et sa factorisation M = LLT ,
alors :

VT V = LLT ⇒ L−1 VT VL−T = I ⇒ UT = L−1 VT (A.1)

Il faut donc procéder à n résolutions du système triangulaire de taille n :

LUT = VT (A.2)

On réalise ainsi une orthonormalisation des vecteurs de départ.

 B
BORNAGE PAR LES VALEURS
PROPRES ÉLÉMENTAIRES

On se place ici dans le cas de figure où les problèmes sont issus d’une analyse par éléments
finis. L’élément fini courant est noté e et les quantités élémentaires associées porteront aussi
l’indice e .

B.1 Système aux valeurs propres généralisé

Le problème est de la forme Kx = λMx, et on défini les matrices diagonales par bloc conte-
nant les quantités élémentaires de la façon suivante :
   
Ke Me
   
 K   M 
K̄ =  e ′  et M̄ =  e ′  (B.1)
 ..   .. 
. .

L’assemblage consiste à utiliser formellement une matrice d’assemblage P booléenne telle

que :

K = PK̄PT (B.2)

Comme M est issue de la même discrétisation, la même matrice d’assemblage est utilisée :

M = PM̄PT (B.3)

Considérons le système, de taille supérieure à n :

K̄x̄ = λ̄M̄x̄ (B.4)

Il englobe le système initial car il suffit de l’écrire dans le sous-espace où x̄ est issu du désas-
semblage d’un vecteur x de raille n : x̄ = PT x 1 . Avec le quotient de Rayleigh :

x̄T K̄x̄
R̄(x̄) = (B.5)
x̄T M̄x̄
on sait que :

λ̄max = max R̄(x̄) (B.6)

x̄6=0

1. Les forces généralisées, elles, s’assembleraient de la façon suivante : f = Pf̄.

 76 BORNAGE PAR LES VALEURS PROPRES ÉLÉMENTAIRES

donc :

x̄T K̄x̄
λ̄max ¾ (B.7)
x̄ 6=0 x̄T M̄x̄

Au vu du problème désassemblé, λ̄max est la valeur maximale des valeurs propres des élé-
ments pris séparément et la relation (B.7) est en particulier vraie pour x̄ = PT x donc :

x̄T PK̄PT x̄ xT Kx
λ̄max ¾ = = R(x) (B.8)
x̄ 6=0 x̄T PM̄PT x̄ xT Mx

donc :

λ̄max ¾ max R(x) = λmax (B.9)

x6=0

B.2 Système aux valeurs propres

Cette fois-ci, le problème est de la forme Ax = λx. Concernant les contributions élémen-
taire, ce problème n’est plus équivalent au problème précédent ; en particulier, si on écrit
A = M−1 K et comme M−1 est pleine, on n’a plus la connectivité éléments finis dans A, et
si on écrit M = I, M n’est pas l’assemblage des matrices élémentaires identité. Le résultat
précédent ne s’applique donc plus ici. Par contre, on a toujours :

λmax = max R(x) (B.10)

x 6=0

où :

xT Ax
R(x) = (B.11)
xT x
donc :

xT Ax
λmax ¾ (B.12)
x 6=0 xT x

Prenons alors comme vecteur x particulier un vecteur nul partout sauf sur les nœuds d’un
élément quelconque e et notons xe sa restriction sur les degrés de liberté de cet élément
(vecteur de petite taille, donc). On a xTe xe = xT x et les contributions élémentaires xT Ax =
xT PĀPT x = xTe Ae xe + a où, si toutes les matrices élémentaires Ae sont positives, a ¾ 0, alors :

xTe Ae xe + a xTe Ae xe
λmax ¾ ¾ (B.13)
xe 6=0 xTe xe xTe xe

donc :

λmax ¾ λ̄max (B.14)

C’est le résultat contraire au précédent, concernant les systèmes aux valeurs propres géné-
ralisés.
 C
TECHNIQUE D’ACCÉLÉRATION
D’AITKEN

Cette technique [67] se sert des itérés précédents pour construire une meilleure nouvelle
approximation. Si on considère une série scalaire ηn qui tend vers η̄, construite par un algo-
rithme du premier ordre avec un taux de convergence κ ; en supposant les itérés ηn , ηn−1 ,
ηn−2 calculés, l’approximation au premier ordre s’écrit :

ηn − η̄ ≈ κ(ηn−1 − η̄)
(C.1)
ηn−1 − η̄ ≈ κ(ηn−2 − η̄)

Ces deux équations permettent de calculer une approximation de η̄, qui va servir pour
construire le nouvel itéré η′n , et κ, qui permet donc d’estimer le taux de convergence en
cours de route. On détermine ainsi l’expression du nouvel itéré :

(ηn−1 − ηn )2
η′n = ηn − (C.2)
ηn − 2ηn−1 + ηn
 D
REPÈRES BIOGRAPHIQUES

Alexander Craig Aitken

Aitken left the Otago Boys’ High School in Dunedin in 1913 having won a
scholarship to Otago University. He began to study languages and mathe-
matics with the intention of becoming a school teacher but his university
career was interrupted by World War I. Back to New Zealand in 1917, Aitken
followed his original intention and became a school teacher. At his old school Otago Boys’
High School. Encouraged by the new professor of mathematics at Otago University, Aitken
came to Scotland in 1923 and studied for a Ph.D. at Edinburgh under Whittaker. In 1925 he
was appointed to Edinburgh where he spent the rest of his life.
• Born : 1 April 1895 in Dunedin, New Zealand.

• Died : 3 Nov 1967 in Edinburgh, Scotland.

André-Louis Cholesky
Cholesky was a French military officer involved in geodesy and surveying
in Crete and North Africa just before World War I. He entered l’École Poly-
technique at the age of 20 and was attached to the Geodesic Section of the
Geographic Service, in June 1905. That was the period when the revision
of the French triangulation had just been decided to be used as the base of a new cadastral
triangulation. Cholesky solved the problem of the adjustment of the grid with the method
now named after him to compute solutions to the normal equations for the least squares
data fitting problem.
• Born : 15 Oct 1875 in Montguyon (Charente-Inférieure).

• Died : 31 Aug 1918.

Johann Carl Friedrich Gauss

Gaussian elimination, which first appeared in the text Nine Chapters of the
Mathematical Art written in 200 BC, was used by Gauss in his work which
studied the orbit of the asteroid Pallas. Using observations of Pallas taken
between 1803 and 1809, Gauss obtained a system of six linear equations in
six unknowns. Gauss gave a systematic method for solving such equations which is precisely
Gaussian elimination on the coefficient matrix.
80 REPÈRES BIOGRAPHIQUES

• Born : 30 April 1777 in Brunswick, Duchy of Brunswick (now Germany).

• Died : 23 Feb 1855 in Göttingen, Hanover (now Germany).

Issai Schur
In 1894, Schur entered the University of Berlin to read mathematics and
physics. Schur made major steps forward in representation theory of groups,
in collaboration with Frobenius, one of his teacher.
In 1916, in Berlin, he built his famous school and spent most of the rest
of his life there. Schur’s own impressive contributions were extended by his students in a
number of different directions. They worked on topics such as soluble groups, combinato-
rics, and matrix theory.
• Born : 10 Jan 1875 in Mogilyov, Mogilyov province, Belarus.

• Died : 10 Jan 1941 in Jerusalem, Palestine.

Joseph-Louis Lagrange
Joseph-Louis Lagrange is usually considered to be a French mathemati-
cian, but the Italian Encyclopedia refers to him as an Italian mathemati-
cian. They certainly have some justification in this claim since Lagrange
was born in Turin and baptised in the name of Giuseppe Lodovico Lagran-
gia.
The papers by Lagrange cover a variety of topics : beautiful results on the calculus of va-
riations, work on the calculus of probabilities. . . In a work on the foundations of dynamics,
Lagrange based his development on the principle of least action and on kinetic energy. The
‘Mécanique analytique’ which Lagrange had written in Berlin, was published in 1788. It sum-
marized all the work done in the field of mechanics since the time of Newton and is notable
for its use of the theory of differential equations. With this work Lagrange transformed me-
chanics into a branch of mathematical analysis.
« Lagrange, in one of the later years of his life, imagined that he had overcome the diffi-
culty (of the parallel axiom). He went so far as to write a paper, which he took with him to
the Institute, and began to read it. But in the first paragraph something struck him which
he had not observed : he muttered : ‘Il faut que j’y songe encore’, and put the paper in his
pocket. » A De Morgan Budget of Paradoxes.
• Born : 25 Jan 1736 in Turin, Sardinia-Piedmont (now Italy).

• Died : 10 April 1813 in Paris, France.

Karl Gustav Jacob Jacobi

Jacobi published three treatises on determinants in 1841. These were im-
portant in that for the first time the definition of the determinant was made
in an algorithmic way and the entries in the determinant were not specified
so his results applied equally well to cases were the entries were numbers
or to where they were functions. These three papers by Jacobi made the idea of a determi-
nant widely known.
 81

• Born : 10 Dec 1804 in Potsdam, Prussia (now Germany).

• Died : 18 Feb 1851 in Berlin, Germany.

Ludwig Philipp von Seidel

Philipp Seidel entered the University of Berlin in 1840 and studied under
Dirichlet and Encke. He moved to Kœnigsberg where he studied under Bes-
sel, Jacobi and Franz Neumann. Seidel obtained his doctorate from Munich
in 1846 and he went on to become professor there.
Seidel wrote on dioptics and mathematical analysis. His work on lens identified mathe-
matically five coefficients describing the aberration of a lens, now called ’Seidel sums’. These
Seidel sums correspond to spherical aberration, coma, astigmatism, Petzval curvature and
distortion. He also introduced the concept of nonuniform convergence and applied proba-
bility to astronomy.

• Born : 24 Oct 1821 in Zweibrœcken, Germany.

• Died : 13 Aug 1896 in Munich, Germany.

James Hardy Wilkinson

Wilkinson won a Foundation Scholarship to Sir Joseph Williamson’s Ma-
thematical School, Rochester at the age of 11. In 1940, Wilkinson began
war work which involved mathematical and numerical work on ballistics.
Wilkinson continued work becoming more involved in writing many
high quality papers on numerical analysis, particularly numerical linear algebra where he
developed backward error analysis methods. He worked on numerical methods for solving
systems of linear equations and eigenvalue problems.
As well as the large numbers of papers on his theoretical work, Wilkinson developed
computer software, working on the production of libraries of numerical routines. The NAG
(Numerical Algorithms Group) began work in 1970 and much of the linear algebra routines
were due to Wilkinson.

• Born : 27 Sept 1919 in Strood, Kent, England.

• Died : 5 Oct 1986 in London, England.

John William Strutt Lord Rayleigh

His first paper in 1865 was on Maxwell’s electromagnetic theory. He wor-
ked on propagation of sound and, while on an excursion to Egypt taken for
health reasons, Strutt wrote Treatise on Sound (1870-1). In 1879, he wrote
a paper on travelling waves, this theory has now developed into the theory
of solitons. His theory of scattering (1871) was the first correct explanation of why the sky is
blue.

• Born : 12 Nov 1842 in Langford Grove (near Maldon), Essex, England.

• Died : 30 June 1919 in Terling Place, Witham, Essex, England.

82 REPÈRES BIOGRAPHIQUES

Helmut Wielandt
Wielandt entered the University of Berlin in 1929 and there he studied ma-
thematics, physics and philosophy. There he was greatly influenced by Schmidt
and Schur. It was on the topic of permutation groups that Wielandt wrote
his doctoral dissertation and he was awarded a doctorate in 1935. At the
end of World War II Wielandt was appointed associate professor at the University of Mainz,
and in 1951, Ordinary Professor at the University of Tœbingen. For 20 years beginning in
1952, Wielandt was managing editor of Mathematische Zeitschrift.

• Born : 19 Dec 1910 in Niedereggenen, Lœrrach, Germany.

• Died : 1984.

Alston Scott Householder

Householder then taught mathematics in a number of different places. He
was awarded a Ph.D. by the University of Chicago in 1947 for a thesis on the
calculus of variations. However his interests were moving towards applica-
tions of mathematics, particularly applications of mathematics to biology.
In 1946, after the end of the war, Householder joined the Mathematics Division of Oak
Ridge National Laboratory. Here he changed topic, moving into numerical analysis which
was increasing in importance due to the advances in computers. He started publishing on
this new topic with Some numerical methods for solving systems of linear equations which
appeared in 1950.
In 1964, Householder published one of his most important books : The theory of matrices
in numerical analysis.

• Born : 5 May 1904 in Rockford, Illinois, USA.

• Died : 4 July 1993 in Malibu, California, USA.

Cornelius Lanczos
In 1938 at Purdue, Lanczos published his first work in numerical analy-
sis. Two years later he published a matrix method of calculating Fourier
coefficients which, over 25 years later, was recognised as the Fast Fourier
Transform algorithm. In 1946, with Boeing, he worked on applications of
mathematics to aircraft design and was able to develop new numerical methods to solve the
problems. In 1949 he moved to the Institute for Numerical Analysis of the National Bureau
of Standards in Los Angeles. Here he worked on developing digital computers and was able
to produce versions of the numerical methods he had developed earlier to program on the
digital computers.

• Born : 2 Feb 1893 in Székesfehérvr, Hungary.

• Died : 25 June 1974 in Budapest, Hungary.

 83

Joseph Raphson
Joseph Raphson’s life can only be deduced from a number of pointers. It is
through the University of Cambridge records that we know that Raphson
attended Jesus College Cambridge and graduated with an M.A. in 1692.
Rather remarkably Raphson was made a member of the Royal Society in
1691, the year before he graduated. His election to that Society was on the strength of his
book Analysis aequationum universalis which was published in 1690 contained the Newton
method for approximating the roots of an equation.
Raphson’s ideas of space and philosophy were based on Cabalist ideas. The Cabala was
a Jewish mysticism which was influential from the 12th century on and for which several
basic doctrines were strong influences on Raphson’s philosophical thinking. The doctrines
included the withdrawal of the divine light, thereby creating primordial space, the sinking
of luminous particles into matter and a cosmic restoration.

• Born : 1648 in Middlesex, England.

• Died : 1715.

Carle David Tolmé Runge

As a German mathematician, physicist, and spectroscopist, he is known to
be the co-developer and co-eponym of the Runge-Kutta method in the field
of numerical analysis. After taking his secondary school teachers examina-
tions he returned to Berlin where he was influenced by Kronecker. Runge
then worked on a procedure for the numerical solution of algebraic equations in which the
roots were expressed as infinite series of rational functions of the coefficients.

• Born : 30 Aug 1856 in Bremen, Germany.

• Died : 3 Jan 1927 in Göttingen, Germany.

Sir Isaac Newton

Isaac Newton is perhaps the best known renaissance scientist today, living
between 1642 and 1727. We think of gravity, celestial mechanics, and cal-
culus when we think of him. He certainly did develop the calculus by buil-
ding upon the ideas of Fermat and Barrow (the person whose chair he took
when he went to Cambridge). But he was not alone in developing calculus. And his focus
was really one of mechanics - how do bodies move. His focus was always on motion and
is reflected in the terminology he chose for calculus - what we call derivatives, he called
fluxions. And his integrals were simply inverse fluxions. You can see how flux was his focus.
He also developed the dot notation for calculus - one dot meant a first derivative, two dots a
second, and so forth. Today this is used, but only with respect to time derivatives. A general
derivative still needs to specify the variable of differentiation and more often than not, time
derivatives are done in this more general manner.
As indicated earlier, Newton and his followers argued vehemently with Huygens and his
followers over the nature of light. Newton subscribed to a corpuscular theory, where he en-
84 REPÈRES BIOGRAPHIQUES

visioned light as small compact bodies of energy. Huygens focussed on the wave like nature
and developed that theory. The diffraction properties of light were so obvious, that Huygens
school eventually won out, and the wave theory of light ruled science for the next three cen-
turies.

• Born : 4 Jan 1643 in Woolsthorpe, Lincolnshire, England.

• Died : 31 March 1727 in London, England.

Leonhard Euler
The publication of many articles and his book Mechanica (1736-37), which
extensively presented Newtonian dynamics in the form of mathematical
analysis for the first time, started Euler on the way to major mathematical
work. He integrated Leibniz’s differential calculus and Newton’s method of
fluxions into mathematical analysis. In number theory he stated the prime number theo-
rem and the law of biquadratic reciprocity. Euler made large bounds in modern analytic
geometry and trigonometry. He was the most prolific writer of mathematics of all time. His
complete works contains 886 books and papers.

• Born : 15 April 1707 in Basel, Switzerland.

• Died : 18 Sept 1783 in St Petersburg, Russia.

Martin Wilheim Kutta

Kutta was professor at the RWTH Aachen from 1910 to 1911. He then be-
came professor at Stuttgart and remained there until he retired in 1935.
In 1901, he had co-developed the Runge-Kutta method, used to solve or-
dinary differential equations numerically. He is also remembered for the
Zhukovsky-Kutta aerofoil, the Kutta-Joukowski theorem and the Kutta condition in aerody-
namics.

• Born : 3 Nov 1867 in Pitschen, Upper Silesia (now Byczyna), Poland.

• Died : 25 Dec 1944 in Frstenfeldbruck, Germany.

 INDEX

A conditions aux limites . . . . . . . . . . . . 14

accélération convergence . . . . . . . . . . . . . . . . . . . . . . 19
d’Aitken . . . . . . . . . . . . . . . . . . . . . . . 20, 77 discrétisation . . . . . . . . . . . . . . . . . . . 8, 25
méthode de l’accélération linéaire 56, problème de mémoire . . . . . . . . . . . . 41
60 super-élément . . . . . . . . . . . . . . . . . . . . 13
méthode de l’accélération moyenne 56 système différentiel du premier ordre
schéma de Newmark . . . . . . . . . . . . . 55 53
Aitken . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 équation
technique d’accélération . . . . . . 20, 77 caractéristique . . . . . . . . . . . . . . . . 26, 56
d’équilibre . . . . . . . . . . . . . . . . . . . . . . . . 66
B d’état . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
Broyden d’admissibilité cinématique . . . . . . 37
Méthode BFGS . . . . . . . . . . . . . . . . . . . 41 d’admissibilité statique . . . . . . . . . . . 37
d’Euler . . . . . . . . . . . . . . . . . . . . . . . . 14–16
C de comportement . . . . . . . . . . . . . . . . 66
Cauchy de la dynamique . . . . . . . . . . . . . . . . . . 67
problème de . . . . . . . . . . . . . . . . . . 43, 52 différentielle . . . . . . . . . . . . . . . . . . . . . . 67
Cholesky . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 différentielle non-linéaire . . . . . . . . 62
factorisation de . . . . . 12, 13, 23, 25, 73 du mouvement . . . . . . . . . . . . . . . . . . . 67
coût Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
CPU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 équation d’ . . . . . . . . . . . . . . . . . . . . 14–16
d’un algorithme . 8–10, 12, 15, 33, 52, schéma d’
54 explicite, 44, 46, 48
d’une méthode . . . . . . . . . . . . . . . . . . . 17 implicite, 47, 50, 62
fonction . . . . . . . . . . . . . . . . . . . . . . . . . . 40
fonction (optimisation) . . . . . . . . . . . 38 F
complexité . . . . . . . . . . . . . . . . . . 9, 12, 13, 32 factorisation
conditionnement . . . . . . . . . . . . . 6, 7, 24, 26 de Cholesky . . . . . . . . 12, 13, 23, 25, 73
convergence . . . . . . . . . . . . . . . . . . . . . . . 17–20 de Crout . . . . . . . . . . . . . . . . . . . . . . 13–15
asymptotique . . . . . . . . . . . . . . . . . . . . . 18 de Gauss . . . . . . . . . . . . . . . . . . . . . . 13, 15
cubique . . . . . . . . . . . . . . . . . . . . . . . . . . 30 LDLT . . . . . . . . . . . . . . . . . . . . . . . . . . 13, 31
gradient conjugué préconditionné 23 LDMT . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
linéaire . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 LU . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11–13
méthode de Jacobi . . . . . . . . . . . . . . . . 28 QR . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31, 32
méthode de Lanczos . . . . . . . . . . . . . . 35
méthode des itérations inverses . . 30
G
Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
méthode des puissances . . . . . . . . . . 28
élimination de . . . . . . . . . . . . . . . . . . . . 10
méthode du gradient . . . . . . . . . . . . . 22
algorithme de Gauss-Seidel . . . . . . . 19
quadratique . . . . . . . . . . . . . . . . . . . . . . 64
factorisation de . . . . . . . . . . . . . . . 13, 15
schéma d’Euler explicite . . . . . . . . . . 44
méthode de Gauss-Seidel . . . . . 18, 19
taux de . . . . . . . . . . . . . . . . . . . . 18, 19, 77
valeurs propres multiples . . . . . . . . . 29 H
Cramer Hessenberg
méthode de . . . . . . . . . . . . . . . . . . . . . . . . 8 forme de . . . . . . . . . . . . . . . . . . . . . . . . . . 33
Crout Householder . . . . . . . . . . . . . . . . . . . . . . . . . . 82
factorisation . . . . . . . . . . . . . . . . . . . . . . 13 algorithme de Householder-Businger-
factorisation de . . . . . . . . . . . . . . . 13–15 Golub . . . . . . . . . . . . . . . . . . . . . . . . 32
Cuthill-MacKee . . . . . . . . . . . . . . . . . . . . . . . . 9 factorisation QR . . . . . . . . . . . . . . . . . . 32
méthode de . . . . . . . . . . . . . . . . . . . . . . 26
E transformation de . . . . . . . . . . . . . 31, 34
éléments finis . . . . . . . 1, 5, 7, 13, 37, 75, 76
code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 J
conditionnement . . . . . . . . . . . . . . . . . . 7 Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80, 81
86 INDEX

algorithme de . . . . . . . . . . . . . . . . . . . . 27 schéma de . . . . . . . . . . . . . . . . . 55, 57, 58

méthode de . . . . . . . . 18–20, 26, 27, 29 stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
convergence, 28 Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80, 83
fluxions . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
K loi de . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
Kahan méthode de . . . . . . . . . . . . . . . 37, 38, 63
théorème de . . . . . . . . . . . . . . . . . . . . . . 20 modèle de . . . . . . . . . . . . . . . . . . . . . . . . 68
Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 Newton-Raphson . . . . . . . . . . . . . 39, 64
Kutta-Joukowski . . . . . . . . . . . . . . . . . . 84 quasi-newton . . . . . . . . . . . . . . . . . . . . . 40
méthode de Runge-Kutta . 51–53, 55, norme . . . . . . . . . . 6–8, 18, 22, 24, 29, 30, 36
62 infinie . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . 83 matricielle . . . . . . . . . . . . . . . . . . . . . . . . . 6
Zhukovsky-Kutta . . . . . . . . . . . . . . . . . 84 normer . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28, 35

L O
Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 orthogonalisation
multiplicateurs de . . . . . . . . . . . . 15, 20 d’un sous-espace . . . . . . . . . . . . . . . . . 12
Lanczos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 de Gram-Schmidt . . . . . . . . . . . . . . . . 73
itérations de . . . . . . . . . . . . . . . . . . . . . . 36 ré-orthogonalisation de grand problème
méthode de . . . . . . . . . . . . . . . 26, 35, 36 35
Leibniz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
linéaire(s) P
équation(s) . . . . . . . . . . . . . . . . . . . 46, 54 pivot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
combinaison(s) . . . . . . . . . . . . 10, 17, 29 Crout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
convergence(s) . . . . . . . . . . . . . . . . . . . 33 facrotisation LU . . . . . . . . . . . . . . . . . . 11
méthode des accélérations . . . . 56, 60 partiel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
récurrence(s) . . . . . . . . . . . . . . . . . . . . . 18 stratégie . . . . . . . . . . . . . . . . . . . . . . . . . . 11
relation(s) . . . . . . . . . . . . . . . . . . . . . . . . 16 total . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
système(s) . . . . 5, 18, 26, 32, 39, 54, 55 préconditionnement . . . . . . . . . . . . . . . . . . 23
gradient conjugué préconditionné 23,
M 41
matrice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 préconditionneur . . . . . . . . . . . . . . . . . 23
bande . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
booléenne . . . . . . . . . . . . . . . . . . . . . . . . 14 R
condensée . . . . . . . . . . . . . . . . . . . . . . . . 11 Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
d’amortissement . . . . . . . . . . . . . . 25, 55 Newton-Raphson . . . . . . . . . . . . . 39, 64
d’assemblage . . . . . . . . . . . . . . . . . . . . . 75 Rayleigh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
d’itération . . . . . . . . . . . . . . . . . . . . . 19, 56 amotissement de . . . . . . . . . . . . . . . . . 57
décalée . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 itérations inverses de . . . . . . . . . . . . . 30
de capacité . . . . . . . . . . . . . . . . . . . . . . . 53 quotient de . . . . . . . . . . . . . . . . . . . 26, 75
de conduction . . . . . . . . . . . . . . . . . . . . 53 Runge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
de masse . . . . . . . . . . . . . . . . . . . . . . 54, 55 méthode de Runge-Kutta . 51–53, 55,
de raideur . . . . . . . . . . . . . . . . . . . . . . . . 16 62, 84
de rigidité . . . . . . . . . . . . . . 13, 14, 38, 55 Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . 83
tangente, 64
de rotation . . . . . . . . . . . . . . . . . . . . . . . 27 S
jacobienne . . . . . . . . . . . . . . . . 38–40, 63 Schur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80, 82
pleine . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 complément de . . . . . . . . . . . . . . . . . . . 11
semi-définie . . . . . . . . . . . . . . . . . . . . . . 17 condensation de . . . . . . . . . . . . . . . . . . 13
sous-matrice . . . . . . . . . . . . . . . . . . 11, 12 Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
symétrique . . . . . . . . . . . . . . . . . . . . . . . 13 algorithme de Gauss-Seidel . . . . . . . 19
symétrique définie positive . . . . . . . 15 méthode de Gauss-Seidel . . . . . 18, 19
triangulaire supérieure . . . . . . . . . . . 11 Southwell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
valeurs propres . . . . . . . . . . . . . . . 20, 25 système
aux valeurs propres . . . . 25, 26, 75, 76
N d’exploitation . . . . . . . . . . . . . . . . . . . . 10
Newmark différentiel
méthode de . . . . . . . . . . . . . . . 55, 57, 61 du premier ordre, 67
INDEX 87

non régulier, 67
dynamique . . . . . . . . . . . . . . . . . . . . . . . 67
linéaire . . . . . . . . . . . . . . . 5, 8, 26, 32, 39
creux, 8
triangulaire, 10, 12
non-linéaire . . . . . . . . . . . . . . . . . . . 18, 37

T
trapèze(s)
méthode des . . . . . . . . . . . 47, 48, 61, 67
schéma des . . . . . . . . . . . . . . . . . . . . . . . 48

W
Wielandt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
méthode de . . . . . . . . . . . . . . . . . . . . . . 30
Wilkinson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
test de . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36