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. Mont-
pellier, France. 2008, pp.98. �cel-00351713v2�

HAL Id: cel-00351713

[Link]
Submitted on 3 Jan 2014

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.

Distributed under a Creative Commons Attribution - NonCommercial - ShareAlike| 4.0

International License
 Méthodes numériques
appliquées à la conception
par éléments finis

David Dureisseix
Université de Montpellier 2 — 2013
Ce document est sous licence Creative Commons 3.0 France:
— paternité;
— pas d’utilisation commerciale;
— partage des conditions initiales à l’identique;
[Link]
 Table des matières

1 Résolution des systèmes linéaires

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

2 Systèmes aux valeurs propres

2.1 Préliminaires 33
2.1.1 Rappels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.1.2 Domaine d’emploi des différentes méthodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.2 Méthode de Jacobi 34
2.3 Méthode des puissances 36
2.4 Méthode des itérations inverses 38
 Table des matières 4

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

2.5.1 Principe et convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.5.2 Exemples d’utilisation dans un code éléments finis . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.6 Transformation de Householder 40
2.6.1 Factorisation QR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.6.2 Détermination des valeurs propres par la méthode QR . . . . . . . . . . . . . . . . . . . . . . . 41
2.6.3 Forme de Hessenberg . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.7 Exemple d’utilisation dans un code éléments finis 44
2.7.1 Réduction dans un sous-espace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.7.2 Technique itérative . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.8 Méthode de Lanczos 44

3 Résolution de systèmes non-linéaires particuliers

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

4 Résolution de systèmes différentiels

4.1 Équation scalaire du premier ordre 53
4.1.1 Problème de Cauchy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.1.2 Ré-écriture du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.1.3 Schémas explicites à un pas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.1.4 Quelques problèmes pratiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.1.5 A-stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.1.6 Méthode des trapèzes généralisée . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.1.7 Méthode à pas multiples : Adams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.1.8 Méthode de prédiction-correction basée sur Adams . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.1.9 Méthode de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.1.10 Autres méthodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.2 Systèmes différentiels du premier ordre 64
4.2.1 Emploi des méthodes précédentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.2.2 Influence de l’ordre d’intégration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.3 Systèmes différentiels du deuxième ordre 66
4.3.1 Méthode de Newmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.3.2 Méthode de Gear implicite à deux pas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.3.3 -méthode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.3.4 ˛-HHT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.3.5 Autres méthodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
5 Table des matières

4.4 Intégration de modèles de comportement 72

4.4.1 Exemple de modèle de plasticité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.4.2 Méthodes avec retour radial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.4.3 Modèle de viscoplasticité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.5 Cas de la dynamique non régulière 77

Bibliographie

A Méthodes particulières
A.1 Factorisation de cholesky et orthogonalisation d’un sous-espace 87
A.2 Technique d’accélération d’Aitken 87

B Bornage par les valeurs propres élémentaires

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

C Repères biographiques

Index
 Introduction

Ce document de travail (avec renvoi à des ressources de références s’il y a besoin d’appro-
fondir) est utile à ceux qui sont intéressés ou amenés à utiliser des méthodes de simulation
(essentiellement 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
différents mais ce document ne se veut pas exhaustif. Il ne traite en particulier pas les cas
suivants :
— la recherche de racine de f .U/ D 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éplacements 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 à [49].
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.

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 xi (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 A1Wn;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 ex-
posant (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.
 Introduction 8

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 à [6, 28] 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) [2] :
— [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]
Sa description est celle de la FAQ (Frequently Asked Questions) :
LAPACK provides routines for solving systems of simultaneous linear equa-
tions, least-squares solutions of linear systems of equations, eigenvalue pro-
blems, and singular value problems. The associated matrix factorizations (LU,
Cholesky, QR, SVD, Schur, generalized Schur) are also provided, as are rela-
ted 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 com-
plex matrices, in both single and double precision.
Citons aussi sa version parallélisée, ScaLAPACK (Scalable LAPACK, [8]) :
The ScaLAPACK project is a collaborative effort involving several institutions
(Oak Ridge National Laboratory, Rice University, University of California, Ber-
keley, University of California, Los Angeles, University of Illinois, University
9 Introduction

of Tennessee, Knoxville), and comprises four components : dense and band ma-
trix software (ScaLAPACK), large sparse eigenvalue software (PARPACK and
ARPACK), sparse direct systems software (CAPSS and MFACT), preconditio-
ners for large sparse iterative solvers (ParPre). The ScaLAPACK (or Scalable
LAPACK) library includes a subset of LAPACK routines redesigned for dis-
tributed 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 minimize the frequency of data move-
ment 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 pro-
cessor.) The fundamental building blocks of the ScaLAPACK library are dis-
tributed memory versions (PBLAS) of the Level 1, 2 and 3 BLAS, and a set
of Basic Linear Algebra Communication Subprograms (BLACS) for commu-
nication tasks that arise frequently in parallel linear algebra computations. In
the ScaLAPACK routines, all interprocessor communication occurs within 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 pos-
sible.
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

eigenvalue problems. The package is designed to compute a few eigenvalues
and corresponding eigenvectors of a general n  n matrix A. It is most appro-
priate 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 n2
floating point operations. This software is based upon an algorithmic variant
of the Arnoldi 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 technique that is suitable for large scale problems. For many stan-
dard 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,
 Introduction 10

portable, and widely available, they are commonly used in the development of
high quality linear 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 D b pour la simulation sur ordinateur ; par exemple en éléments finis [54], une analyse
statique conduit à résoudre un système que l’on notera Ku D 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 ou la recherche d’une base modale sont souvent
basé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 D 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 ¤ 0
— A symétrique , A D A> avec 8x; y; x
y D x> y; x> Ay D y> A> x
— A positive , 8x; x> Ax > 0
— A définie , x> Ax D 0 ) x D 0
— A mal conditionnée, typiquement cond K D 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 à [4, 59].

1.1 Conditionnement
1.1.1 Définition
Pour A symétrique et inversible,
1
cond A D kAk
kA k (1.1)
Avec une norme matricielle cohérente avec celle définie pour les vecteurs :
kAxk
kAk D sup (1.2)
kxk¤0 kxk
P
où kxk2 D x> x D i xi2 est la norme euclidienne, on a :
jjmax
cond A D (1.3)
jjmin
 Section 1.1 | Conditionnement 12

où les  sont les valeurs propres de A. En général, il est très coûteux d’atteindre le condition-
nement 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éci-
sion (cumul d’arrondis numériques). Pour les illustrer, on est amené à considérer un système
perturbé (par des perturbations généralement supposées petites) AQ Q x D b.
Q

Exemple 1.1 Considérons le système Ax D b suivant :

      
1 1 x1 2 2
D ) xD (1.4)
1 1 C " x2 2 0

Le polynôme caractéristique 2 .2 C "/ C " D 0 permet de déterminer les valeurs propres, à

savoir :
p
2 C " ˙ 4 C "2
1;2 D (1.5)
2
Avec " petit, le conditionnement est donc très grand puisque cond A  4=". Considérons mainte-
nant le système AQQ x D bQ perturbé tel que :
      
1 1 x1 2 2 ˛"
D ) xD ˛ (1.6)
1 1 C " x2 2C˛ "

Ce résultat entraîne la discussion suivante :

— si ˛  "  1 alors x  Œ2 0> ;
— si ˛  "  1, ce qui est tout à fait possible avec les arrondis numériques, la solution est
très différente puisque x  Œ1 1> .

Exemple 1.2 Examinons le système suivant :

     
1;2969 0;8648 0;8642 0;8642
AD bD xD (1.7)
0;2161 0;1441 0;1440 0;1440

Le conditionnement de A est 3;3
108 ; le résidu est petit puisque b Ax D Œ 10 8 10 8 >

alors que la solution exacte est très différente de x puisque c’est Œ2 2> . 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.1= h2 /. La figure 1.1 illustre cela dans un cas bidimensionnel très
simple. Le conditionnement 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.

1.1.2 Influence des perturbations

Considérons le système AQx D b.Q Si bQ D b C ıb avec kıbk  kbk, a-t-on kıxk  kxk
pour xQ D x C ıx ? Calculons simplement, en utilisant l"inégalité de Milkowski :

Ax D b ) kAx D bk ) kAk
kxk > kbk (1.8)

13 Chapitre 1 | Résolution des systèmes linéaires

105
L
h
104

cond K
103

102

10
10 2 10 1 1
h=L
(a) maillage (b) conditionnement

Figure 1.1 — Conditionnement en éléments finis

autrement dit :
1 kAk
6 (1.9)
kxk kbk
et maintenant :
1 1
[Link]ıx/ D bCıb ) Aıx D ıb ) ıx D A ıb ) kıxk 6 kA k
kıbk (1.10)

soit :
1
kıxk 6 kA k
kıbk (1.11)

par conséquent :

kıxk kıbk
6 cond A (1.12)
kxk kbk

Q x D b. Si A
Considérons maintenant le système perturbé AQ Q D A C ıA et xQ D x C ıx, alors :

.A C ıA/.x C ıx/ D b (1.13)

implique :

Aıx C ıA.x C ıx/ D 0 (1.14)

soit :
1
ıx D A ıA.x C ıx/ (1.15)
 Section 1.2 | Techniques de stockage et coût 14

ce qui entraîne finalement :

1
kıxk 6 kA k
kıAk
kx C ıxk (1.16)

autrement dit :
kıxk kıAk
6 cond A (1.17)
kx C ı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
système est n D 103 , le stockage de A pleine nécessite 8 Mo de mémoire ; avec n D 106 ,
8 000 Go sont nécessaires.
Pire encore, le nombre d’opérations pour résoudre le système Ax D 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 flot-
tante de nop  .n C 1/Š, ce qui devient rapidement inextricable, même avec les ordinateurs
les plus puissants [15, 48].
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, lb ) ;
— 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 :
— Cuthill-MacKee inverse [12] ;
— nested dissection [24] ;
— degré minimum [25, 64]...

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
1. en anglais, un système creux est dit sparse
15 Chapitre 1 | Résolution des systèmes linéaires

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 problè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 deux opérations en fait : addition et multiplica-
tion), dscal (série de multiplications), dgemv (calcul d’un produit matrice A par vecteur x,
pondéré par un coefficient ˇ et ajouté à un précédent vecteur y, lui même multiplié par un
coefficient ˛) sont appelées avec la librairie BLAS, déjà mentionnée. L’opération dénommée
˛y C ˇAx est donc l’équivalent de dgemv, programmé directement en Fortran. À partir de
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’utilisa-
tion de routines standards augmente la lisibilité du programme, permet de programmer à un
niveau 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
2. les performances du système, du compilateur, les accès aux disques
 Section 1.3 | Méthodes directes 16

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 C ˇAx 4,95 (5,2) 3,4 (3,9) 0,84
dgemv - 2 (2,3) 1,1

Tableau 1.1 — Estimation du coût des opérations élémentaires

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 [17].

1.3.1 Méthode d’élimination de Gauss

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

Élimination de Gauss par bloc

Écrivons le système à résoudre sous une forme de blocs :

    
A11 A12 x1 b1
D (1.18)
A21 A22 x2 b2

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

x1 D A111 .b1 A12 x2 / (1.19)

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

comme seule inconnue :

A22 A21 A111 A12 x2 D b2 A21 A111 b1 soit Sx2 D c (1.20)

17 Chapitre 1 | Résolution des systèmes linéaires

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.19). Le système de départ (1.18) est équivalent à :
    
A11 A12 x1 b1
D (1.21)
0 S x2 c

et on peut le réécrire :
      
I 0 A11 A12 x1 I 0 b1
D (1.22)
A21 A111 I 0 S x2 A21 A111 I c

Utilisation pour la factorisation LU

Si on prend comme bloc 11 la première équation scalaire du système, il faut avoir A11 ¤ 0
(c’est le premier pivot). En continuant sous la forme de l’équation (1.22), récursivement
sur les sous-matrices qui restent, on obtient la factorisation A D 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/ .

.k 1/
Conditions d’emploi Si à chaque étape de la factorisation, le pivot Akk est non nul,
on peut progresser sans problè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 assemblé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 :

jAl;k j D max jAm;k j (1.23)

mDk;:::;n

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 :

jAl;c j D max jAm;p j (1.24)

m;pDk;:::;n
 Section 1.3 | Méthodes directes 18

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 b0 ,
on résoud :

LUx0 D b0 (1.25)

en considérant, dans un premier temps, Ly0 D b0 étape appelée descente, puis puis Ux0 D y0 ,
appelée montée.
À chaque fois, il ne s’agit plus que de résoudre un système triangulaire (supérieur ou
infé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 D [Link]/ D det L det U D det U D Uk;k (1.26)
kD1;n

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

1.3.2 Factorisation de Cholesky

Cette technique a été décrite pour la première fois dans [7] en 1924. Elle considère la facto-
risation LU quand A est symétrique et inversible, c’est-à-dire A D B> 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 B> B requière environ n3 =6 opérations. Pour A

bande, la factorisation B> B requière approximativement nlb2 =2 opérations, une montée (ou
une descente) coûte de l’ordre de nlb opérations. Le stockage nécessite aux alentours de
nlb réels (généralement 8 octets). De façon courante, on a lb  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.1.

1.3.3 Factorisation de Crout

C’est la méthode la plus employée. La factorisation se met sous la forme A D LDL> avec L
triangulaire inférieure à diagonale unitaire et D diagonale. Cette factorisation s’applique
pour A symétrique, définie positive, mais traite aussi les cas où les valeurs propres sont
négatives ou nulles.
lb 1 1
3. n D 10 à 100
 19 Chapitre 1 | Résolution des systèmes linéaires

boucle sur j D 1; 2; : : : ; n
boucle sur i D 1; 2; : : : ; j 1
vi Aj;i Ai;i
fin
vj Aj;j Aj;1Wj 1 v1Wj 1
Aj;j vj
1
Aj C1Wn;j vj
.Aj C1Wn;j Aj C1Wn;1Wj 1 v1Wj 1 /
fin

Algorithme 1.1 — Factorisation LDL stockée pleine

Complexité 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.1 présente cette factorisation dans le cas d’une matrice symétrique A sto-
cké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 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

zone d’intérêt particulier

i
S

Figure 1.2 — Utilisation d’un super-élément

notion de super-élément. Considérons une partie d’une structure qui représente une zone
d’intérêt particulier (par exemple, une zone qui reste en élasticité alors que le reste de la
structure plastifie, où une zone dans laquelle des fonctions de bases particulières ont été
ajoutées comme dans la méthode de partition de l’unité...) et qui a été maillée par éléments
finis, figure 1.2. 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 Kib
KD (1.27)
Kbi Kbb
 Section 1.4 | Quelques particularités en éléments finis 20

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

S D Kbb Kbi Ki i1 Kib (1.28)

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.

Conservation de la bande

La figure 1.3 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 LDM> . D’autres variantes des techniques précédentes permettent de tirer
mieux parti des conditions de stockage, en particulier pour optimiser les échanges de don-
nées avec la partie des informations stockées sur disque, comme la méthode frontale [38].
D’autres versions encore permettent de réaliser des factorisations en parallèle de façon assez
performante, comme les méthodes multi-frontales [16].

k
coefficients
modifiés à
l’étape k

Figure 1.3 — É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.
21 Chapitre 1 | Résolution des systèmes linéaires

Séparation des degrés de liberté

Le problème correspondant de minimisation du potentiel sous contrainte est le suivant :
1
min J1 .u/ D u> Ku f> u (1.29)
CuDud 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.18) en séparant les de-
grés de liberté imposés u2 D Cu D ud . Le problème (1.29) est ainsi équivalent à :
1
min u> K11 u1 f> >
1 u1 C u1 K12 ud C constante (1.30)
u1 2 1

dont l’équation d’Euler correspondante est K11 u1 D 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 terme K12 , et qui entraîne un surcoût important lorsque le nombre de degrés de
liberté bloqués est grand.

Multiplicateurs de Lagrange
Une autre technique consiste à reformuler le problème (1.29) 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/ l> .Cu ud / (1.31)

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 C> u f
D (1.32)
C 0 l ud
„ ƒ‚ …
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 D (1.33)

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.
4. ce sont les réactions à l’appui qui permettent d’avoir un déplacement bien égal à ud
 Section 1.4 | Quelques particularités en éléments finis 22

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 :
0 1> 2 30 1
l 1 1 1 l1  >  
1 @ 1A 4 ud l1
extr J1 .u/ C l2 1 1 15 @ l2 A (1.34)
u;l1 ;l2 2 ud l2
u 1 1 0 u2
„ ƒ‚ …
1
.l1 l2 /2 C .l1 C l2 /> .Cu ud /
2

dont les équations d’Euler sont :

Ku D f C> .l1 C l2 /
l1 D l2 (1.35)
Cu D ud

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

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’ensemble structure + ressort extérieur est alors :
1
J1 .u/ C .Cu ud /> [Link] ud / (1.36)
2
Les nouvelles équations d’Euler de cette formulation sont alors :

.K C C> kC/u D f C C> kCud (1.37)

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 D 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.
 23 Chapitre 1 | Résolution des systèmes linéaires

1.4.3 Système singulier : semi défini positif

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 >
u Ku D 0 (1.38)
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 D 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 R> b D 0. Exemple : les efforts ne sollicitent pas les modes
de solide rigide.
— 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 C
R˛ où ˛ est quelconque.

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

AC telle que AAC A D A. Si A est inversible, il n’y a qu’une pseudo-inverse AC D A 1 .
Pour A symétrique, toutes les pseudo-inverses sont de la forme AC C B> R> C RB, avec B
quelconque [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.18), où
A11 est la plus grande sous-matrice régulière. On a alors S D A22 A21 A111 A12 D 0. Une
pseudo-inverse de A est :
 1 
A11 0
(1.39)
0 0

et la solution est :
   1   
x1 A11 b1 A12
D C x2 soit x D x0 C R˛ (1.40)
x2 0 Id

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
 Section 1.5 | Méthodes itératives stationnaires 24

une combinaison linéaire des précédentes. L’inconnue correspondante est repérée et fixée
à 0 [voir x0 dans (1.40)], 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.40)].

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). Éventuelle-
ment, comme en éléments finis, on peut le faire sous forme élémentaire :
X
Aw D Ae we (1.41)
e

— on construit une suite de solutions approchées :

x.k/ ! A 1
b (1.42)
k!1

à partir d’une récurrence linéaire x.k/ D Bx.k 1/ C 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 D Bx C c ) .x.k/ x/ D B.x.k 1/
x/ (1.43)

Pour avoir kx.k/ xk !k!1 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
asymptotique est défini par :

D ln .B/ (1.44)

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
.x.k/ x/ D Bk .x.0/ x/ ) 6 kBk k ) ln 6 k (1.45)
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 [29].
5. Cette caractéristique n’existe pas pour les méthodes non stationnaires.
 25 Chapitre 1 | Résolution des systèmes linéaires

1.5.1 Méthode de Jacobi

Pour les systèmes linéaires, cette méthode, d’algorithme 1.2, est décrite pour mémoire seule-
ment, car elle possède un faible taux de convergence et dans nos applications, est remplacée
par d’autre techniques plus efficaces. Elle illustre cependant bien les techniques de point

Vecteur de départ x.0/

tant que R > "
itération sur k
boucle sur i D 1; n
.k 1/ .k 1/
xN i Ai;1Wi 1 x1Wi 1 C Ai;i C1Wn xi C1Wn
xN i .bi xN i /=Ai;i
fin
x.k/ xN
fin

Algorithme 1.2 — Méthode de Jacobi

fixe : si on suppose connues toutes les inconnues sauf l’inconnue i, cette dernière peut être
déterminée à l’aide de l’équation i. Procédant ainsi sur toutes les inconnues, on obtient un
nouvel itéré xk à partir de l’ancien 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 D .I D 1
A/xk 1
CD 1
b (1.46)

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

ici B D I D 1 A. On peut noter l’utilisation d’un vecteur de stockage xN supplémentaire.

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

h2
J D 1 C O.h4 / (1.47)
2
d’où un taux de convergence faible lorsque le maillage est raffiné (donc lorsque h est faible) :

h2
 (1.48)
2

1.5.2 Méthode de Gauss-Seidel

L’idée est similaire à celle de la méthode de Jacobi, à ceci près que l’on utilise les nouvelles
valeurs estimées dès qu’elles sont disponibles, et non pas à la fin de l’itération. Sous forme
matricielle, cela se traduit par :

xk D .D E/ 1
.F xk 1
C b/ (1.49)
 Section 1.5 | Méthodes itératives stationnaires 26

où E est la partie triangulaire inférieure de A, et F sa partie triangulaire supérieure. La

matrice d’itération est ici :
1 1 1 1
B D .D E/ F D .I D E/ D F (1.50)
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.

Vecteur de départ x.0/

tant que R > "
itération sur k
boucle sur i D 1; 2; : : : ; n
.k 1/
xi 0
.k 1/
ˇ Ai;1Wn x1Wn
.k/
xi .bi ˇ/=Ai;i
fin
fin

Algorithme 1.3 — Méthode de Gauss-Seidel

1.5.3 Méthode de relaxation

Cette technique, SOR pour successive over relaxation, est attribuée à Southwell [63] et
date de 1946. Elle utilise comme base la méthode de Gauss-Seidel ainsi qu’une technique
maintenant souvent utilisée avec d’autre approches : la technique de relaxation.
Le principe, expliqué par l’algorithme 1.4, consiste à modifier (par relaxation) l’itéré
x fourni par l’algorithme de Gauss-Seidel à chaque itération. Pour cela, il faut conserver
k

l’itéré précédent xk 1 et modifier l’itéré courant comme suit :

xk !xk C .1 !/xk 1
(1.51)
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 :
1 1 1
B D .I !D E/ Œ.1 !/I C !D F (1.52)
Le théorème de Kahan donne une condition nécessaire de convergence : 0 < ! < 2. Si
! D 0, l’algorithme ne progresse plus. Le théorème d’Ostrowski-Reich permet de conclure
à la convergence 8! 2 0 I 2Œ quand A est symétrique, définie positive.
Le choix de ! est fait à partir d’heuristiques, comme ! D 2 O.h/ pour une discrétisa-
tion éléments finis dont la taille de maille est h. Il existe aussi des techniques d’optimisation
de ! en cours de route, comme la technique d’accélération d’Aitken, voir annexe A.2.
On peut aussi s’inspirer de résultats dans des cas particuliers. Par exemple dans un cas
tridiagonal bloc, on montre que la valeur optimale de ! est :
2
!opt D q (1.53)
1C 1 J2
 27 Chapitre 1 | Résolution des systèmes linéaires

Vecteur de départ x.0/

tant que R > "
itération sur k
boucle sur i D 1; 2; : : : ; n
.k 1/
˛ xi
.k 1/
xi 0
.k 1/
ˇ Ai;1Wn x1Wn
x.k/
i .1 !/˛ C !ˇ
fin
fin

Algorithme 1.4 — Méthode SOR

La valeur de J est ainsi à estimer sur le problème en question. On a alors SOR D !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 de l’algorithme 1.5.

Vecteur de départ x.0/

tant que R > "
itération sur k
boucle sur i D 1; 2; : : : ; n
.k 1/ .k 1/
xN i Ai;1Wi 1 x1Wi 1 C Ai;i C1Wn xi C1Wn
xN i .bi xN i /=Ai;i
fin
x.k/ x.k 1/ C !.Nx x.k 1/ /
fin

Algorithme 1.5 — Méthode de Jacobi relaxée

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
D (1.54)
C > 0 l g

Pour que le problème continue d’être inversible, il suffit que C soit injective.
La méthode d’Uzawa développée en 1958 consiste à résoudre itérativement de façon
découplée, comme le montre l’algorithme 1.6 où ˛ est un paramètre de la méthode. Pour
analyser cet algorithme, on peut condenser la variable x.k/ sur la variable l.k/ pour trouver :

l.k/ D .I ˛C> K 1
C/l.k 1/
C ˛.C> K 1
f g/ (1.55)
 Section 1.6 | Méthodes itératives non-stationnaires 28

Vecteurs de départ x.0/ , l.0/

tant que R > ", itérer sur k
résoudre Kx.k/ D b Cl.k 1/
réactualiser l.k/ l.k 1/ C ˛.C> x.k/ c/
fin

Algorithme 1.6 — Méthode d’Uzawa

avec B D I ˛C> K 1 C. Si on appelle 1 ; : : : ; p , les valeurs propres de la matrice

C> K 1 C, on en déduit l’équivalence suivante :
%.B/ < 1 , 8i; j1 ˛i j < 1 (1.56)
Avec C> K 1C symétrique, définie positive, i > 0 donc il faut prendre :
2
0<˛< (1.57)
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 :
rDb Ax (1.58)
Il s’agit maintenant d’avoir r D 0. L’idée est ici de se trouver avec un problème de plus

b Ax

x
U

Figure 1.4 — Réduction dans un sous-espace

petite taille, de structure plus simple (par exemple, tridiagonal) par plusieurs étapes de pro-
jection.
Considérons un sous-espace U de l’espace dans lequel la solution est cherchée. Le pro-
blème approché est défini par x D QOx 2 U avec b Ax ? U , comme décrit sur la figure 1.4.
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 à :
Q> .b AQOx/ D 0 ) .Q> AQ/Ox D Q> b; soit O x D bO
AO (1.59)
Une méthode de projection peut donc se schématiser de la façon suivante :
 29 Chapitre 1 | Résolution des systèmes linéaires

1. r b Ax
O
2. xO D A 1 Q> r

3. x x C QOx
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.kC1/ x.k/ C .k/ d.k/ (1.60)

où d.k/ est la direction de recherche dans la direction du gradient :

d.k/ D Ax.k/ b (1.61)

.k/ est le paramètre d’optimalité de la nouvelle solution dans la direction d.k/ . Il se calcule
comme :
d.k/> .b Ax.k/ / d.k/> d.k/
.k/ D D (1.62)
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
choisit pour Q le sous-espace généré par le seul vecteur résidu Q D 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 [32] en 1952. Elle est utilisée
dans le cas où A est symétrique, définie positive.
L’initialisation consiste à choisir un vecteur de départ x.0/ , à calculer le résidu associé
r.0/ D b Ax.0/ et de prendre comme première direction de recherche d.0/ D r.0/ . Ensuite,
les itérations sont de la forme :
1. itéré x.kC1/ D x.k/ C .k/ d.k/ avec le pas optimal :

r.k/> d.k/ r.k/> r.k/

.k/ D D (1.63)
d.k/> Ad.k/ d.k/> Ad.k/

2. résidu r.kC1/ D b Ax.kC1/ ;

3. direction de recherche d.kC1/ D r.kC1/ C .k/ d.k/ qu’on orthogonalise par rapport
à la précédente :

r.kC1/> Ad.k/ r.kC1/> r.kC1/

d.kC1/> d.k/ D 0 ) .k/ D D (1.64)
d.k/> Ad.k/ r.k/> r.k/
4. itération suivante.
 Section 1.6 | Méthodes itératives non-stationnaires 30

Propriétés Si i < j , alors r.i/> r.j / D 0 D d.i/> r.j / D d.i/> Ad.j / D r.i/> Ad.j / . On
construit donc des vecteurs orthogonaux au sens de A, d.0/ ; d.1/ ; : : : qui forment 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 et on a la propriété de convergence suivante :
h pcond A 1 i.k/
.k/
kx xkA 6 2 p kx.0/ xkA (1.65)
cond A C 1
d’où, en éléments finis :
 1 
 DO p D O.h/ (1.66)
cond A

1.6.4 Préconditionnement
Si on peut avoir facilement une approximation spectrale M de A, l’idée est de chercher à
résoudre :
M 1 Ax D M 1 b (1.67)
où M est symétrique, définie positive. Ce problème devient non symétrique, mais on peut
lui substituer :
1=2 1=2
.M AM /M1=2 x D M 1=2
b (1.68)
Le choix du préconditionneur dépend du problème à traiter. Dans le meilleur des cas, il est
facile à inverser et représente bien la matrice A. En fait, si M D I il est facile à inverser,
mais est très différent de A, et si M D 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 
 DO p (1.69)
cond M 1 A
Comme exemples de préconditionneurs classiques, on a la factorisation de Cholesky incom-
plète, SSOR, etc. L’algorithme 1.7 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 [61] ;
— pour A symétrique mais non définie : MinRES (résidu minimum, minimal residual),
SymmLQ (symmetric LQ) [55] ;
— pour A non symétrique connaissant A> : QMR (quasi minimal residual) [22] ;
— pour A non symétrique : CGS (conjugate gradient square) [62], BiCGStab (bi-
conjugate gradient stabilized) [65], GMRES (generalized minimal residual) [60],
GMRES restarted...
 31 Chapitre 1 | Résolution des systèmes linéaires

Vecteur de départ x.0/

Résidu d’équilibre r.0/ b Ax.0/
tant que R > ", itérer sur k
résoudre Mz.k 1/ D r.k 1/
>
k 1 r.k 1/ z.k 1/
si k D 1 alors
d.1/ z.0/
sinon
k 1
ˇk 1 k 2
d.k/ z.k 1/ C ˇk 1 d.k 1/
fin
q.k/ Ad.k/
k 1
˛k .k/ > .k/
d q
x.k/ x.k 1/ C ˛k d.k/
r.k/ b.k/ Ax.k/ D r.k 1/
˛k q.k/
fin

Algorithme 1.7 — Gradient conjugué préconditionné

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/ D x.k/ x soit
suffisamment petit. Malheureusement, x est la solution à convergence et on ne peut atteindre
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/ D b Ax.k/ D Ae.k/ (1.70)
Un critère souvent employé est le suivant :
kr.k/ k
<" (1.71)
kAx.k/ k C 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 6 kA 1
k (1.72)
donc :
r.k/ e.k/
<" ) < " cond A (1.73)
b x
et, en tout état de cause, le paramètre " devrait être choisi aussi fonction du conditionnement
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,
.k/
kr 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 carac-
téristiques propres, c’est-à-dire les couples .; x/ où x ¤ 0 tels que Ax D 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 D ! 2 Mx, K symétrique,
définie positive et M symétrique positive, avec x ¤ 0. Dans ce cas, on peut utiliser la fac-
torisation de Cholesky de M, soit M D LL> et avec y D L> x, on est ramené au problème
précédent L 1 KL > y D ! 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 [26,
27, 58], flambage...). Dans le cas où il y a de l’amortissement, le système à résoudre est :

MRx C CPx C Kx D f.t/ (2.1)

La matrice d’amortissement, C, est souvent mal maîtrisée et si l’amortissement est peu

influent à basse fréquence, ce qui n’est pas évident pour tous les cas, on est ramené à :

MRx C Kx D 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 [3].

2.1 Préliminaires
2.1.1 Rappels
Matrice quelconque

 est valeur propre de A ssi det.A I/ D 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.

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 ;
— 9B symétrique, définie positive telle que A D 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.
 Section 2.2 | Méthode de Jacobi 34

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

Le quotient de Rayleigh est défini comme étant, pour x ¤ 0 :
x> Ax x> Kx
R.x/ D ou R.x/ D (2.3)
x> x x> Mx
Il vient alors l’équivalence suivante :

R est extremum en x ¤ 0
.; x/ est caractéristique propre de A , (2.4)
R.x/ D 
La i e valeur propre i est le minimum de R.x/; x ¤ 0 sur le sous espace orthogonal aux
précédents vecteurs propres. De plus, en dimension finie, la valeur propre maxi est n D
maxx¤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 D 250 à n D 25 000.
La première méthode est celle de l’équation caractéristique. On peut en effet rechercher
les zéros du polynôme caractéristique de la matrice, mais si les coefficients de ce poly-
nô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
employé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 permettent 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.
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 [39], 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.
35 Chapitre 2 | Systèmes aux valeurs propres

tant que R > ", itérer sur k

choix du terme à annuler Ap;q
.k 1/ .k 1/
Aq;q Ap;p
 .k 1/
2Ap;q
calcul de la tangente t p
calcul du cosinus c 1= 1 C t 2
calcul du sinus s ct
A.k/ A.k 1/
A.k/
p;1Wn
.k 1/
cAp;1Wn .k 1/
sAq;1Wn
A.k/
q;1Wn
.k 1/
sAp;1Wn .k 1/
C cAq;1Wn
A.k/
p;p
.k 1/
Ap;p .k 1/
tAp;q
A.k/
q;q
.k 1/ .k 1/
Aq;q C tAp;q
A.k/
p;q 0
A.k/
q;p 0
fin

Algorithme 2.1 — Méthode de Jacobi

Dans cette méthode, on cherche à ramener la matrice A à une forme diagonale pour trou-
ver tous ses valeurs propres et vecteurs propres, et ce par une séquence de transformations
orthogonales :
>
A.kC1/ D R.k/ A.k/ R.k/ (2.5)
>
avec les rotations R.k/ vérifiant R.k/ R.k/ D I, transformations qui sont des changements
de base de type isométries ne modifiant pas les valeurs propres.
On commence par chercher l’élément hors diagonale de module maximal :
.k/
jA.k/
p;q j D max jAi;j j (2.6)
i¤j

On construit alors la matrice de rotation d’un angle  ad hoc dans le plan .ep ; eq / de façon
à annuler dans l’itéré suivant les termes A.kC1/
q;p et A.kC1/
p;q :

R.k/ D 1 C .cos  1/.ep e> >

p C eq eq / C .sin  1/ep e>
p .sin  C 1/eq e>
q (2.7)
avec :
.kC1/ 1   1
Ap;q D sin 2 A.k/
p;p A.k/
q;q C cos 2 A.k/
p;q D 0 (2.8)
2 2
ce qui entraîne :
.k/ .k/
Ap;p Aq;q
cotan 2 D .k/
(2.9)
Ap;q
Seuls changent les termes sur les lignes ou colonnes p ou q. Dans l’algorithme 2.1, la
tangente t est calculée comme suit :

1 p si  D 0
tD (2.10)
 C sgn./ 1 C  2 sinon
 Section 2.3 | Méthode des puissances 36

Au cours des itérations, un terme nul peut redevenir non nul mais :

.k/ 2
X
a.k/ D Ai;j ! 0 (2.11)
k!1
i¤j

En effet :
X .kC1/
X  .kC1/ .kC1/
 X  .kC1/ .kC1/

a.kC1/ D ŒAi;j 2 C ŒAi;p 2 C ŒAi;q 2 C ŒAp;j 2 C ŒAq;j 2
i ¤j;p
„ ƒ‚ … i ¤p;q „ ƒ‚ … j ¤p;q „ ƒ‚ …
.k/
q;j ¤p;q ŒAi;j 2 ŒA.k/ 2 .k/ 2
ŒA.k/ 2 .k/ 2
i;p  C ŒAi;q  p;j  C ŒAq;j 

D a.k/ 2ŒA.k/
p;q 
2
(2.12)

donc a.kC1/ < 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.kC1/ 6 1 a.k/ et  (2.13)
n.n 1/ n2

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

x .0/ , kx .0/ k D 1
tant que R > ", itérer sur k
v Kx .k 1/
résoudre Mx.k/ D v
˛ kx.k/ k
.k/
 R.x.k/ / D .x.k/ > Kx.k/ /=.x.k/ > Mx.k/ /
normalisation x.k/ x.k/ =˛
fin

Algorithme 2.2 — Méthode des puissances

départ x.0/ puis les itérés successifs par récurrence de la manière suivante :

x.kC1/ D Ax.k/ (2.14)

Si .vi /i est la base propre, on peut écrire les itérés dans cette base. En particulier :
X
x.0/ D ˛i vi (2.15)
i

alors :
X  n
X   k 
.k/ i
x D ki ˛i vi D k1 ˛1 v1 C ˛i vi (2.16)
1
i iD2
37 Chapitre 2 | Systèmes aux valeurs propres

et si ˛1 ¤ 0 puisque ji j=j1 j < 1 alors :

x.k/ ! k1 ˛1 v1 (2.17)

k!1

En pratique, il est utile de normer x.k/ à chaque itération pour éviter l’explosion et dans
ce cas obtenir
P x.k/ ! v1 . Pour comprendre le comportement de la convergence, prenons
x.0/ D i vi , alors :

n 
x.k/ X i k
D v1 C vi (2.18)
k1 iD2
1

donc :
2 n 31=2
x.kC1/ X i 2.kC1/
kC1 v1 6
7

6 iD2 1
1 6 7 2
(2.19)
7
x.k/ D 6 n 7 !
6 X i 2k 7 k!1 1
 
v1
4 5
k
1 1
iD2

Le taux de convergence est alors estimé comme suit :

2 1
 ln D 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
combinaison linéaire des vecteurs propres associés.
Dans le cas de valeurs propres proches, 2 D 1 Cı, 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.kC1/ k=kx.k/ k
se stabilise, et pour la convergence sur les vecteurs propres, kx.kC1/ x.k/ k < " ou :

kAx.kC1/ .kC1/ x.kC1/ k

< "0 (2.21)
kAx.kC1/ C .kC1/ x.kC1/ k

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

 Section 2.4 | Méthode des itérations inverses 38

v1
x

Figure 2.1 — Projection orthogonale à v1

Obtention des modes suivants

La technique utilisée est la déflation orthogonale : on prend x.0/ orthogonal à v1 . Pour éviter
les pertes d’orthogonalité par cumul d’arrondis, on réorthogonalise à chaque itération. Pour
cela, on utilise un projecteur comme décrit sur la figure 2.1 :

v1 v>
1
P1 D I (2.22)
1 v1
v>
d’où :
v>
1x
P1 x D x v1 (2.23)
1 v1
v>
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.

2.4 Méthode des itérations inverses

Comme décrit dans l’algorithme 2.3, on choisit encore un x.0/ arbitraire : les solutions clas-
siques sont de prendre un vecteur de composantes aléatoires, normé, quitte à essayer avec
plusieurs vecteurs de départ pour vérifier la convergence vers la même valeur. On résoud

x.0/ , kx.0/ k D 1
tant que R > ", itérer sur k
v Mx.k 1/
résoudre Kx.k/ D v
˛ kx.k/ k
.k/
 R.x.k/ / D .x.k/ > Kx.k/ /=.x.k/ > Mx.k/ /
normalisation x.k/ x.k/ =˛
fin

Algorithme 2.3 — Itérations inverses

ensuite Ax.kC1/ D x.k/ et on norme x.kC1/ . Il s’agit bien de la méthode des puissances,
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 :
jn 1 j
  ln (2.24)
jn j
 39 Chapitre 2 | Systèmes aux valeurs propres

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? D A I, on s’intéresse donc à A? x D ? x et ? D  .
Comme on converge vers la valeur propre j? j la plus petite, on va donc trouver i la plus
proche de . Le taux de convergence est :
 h j j ji j i
i
D ln inf ; (2.25)
ji 1 j jiC1 j

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.kC1/ D 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 vecteur 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 conver-
gence cubique :

kx.kC1/ vk D [Link].k/ vk3 / (2.26)

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

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

On considère le système Kx D ! 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 !02 M
2. factorisation LDL> 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 D ? Mx pour
obtenir .? ; x/
4. ! 2 !02 C ?
5. le numéro du mode trouvé est N? si ! 2 < !02 ou N? C 1 si ! 2 > !02 .
 Section 2.6 | Transformation de Householder 40

Problème 2 On veut trouver tous les modes dans l’intervalle Œ!1 ; !2 . On procède alors
de la façon suivante :
1. pour trouver le nombre de modes dans Œ!1 ; !2 , on compte le nombre de termes
strictement négatifs de D dans les factorisations de K !12 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

2.6.1 Factorisation QR
On cherche à mettre A sous la forme A D QR où R est triangulaire supérieure et Q or-
thogonale, par applications successives de transformations de Householder à A ! HA [34,
67]. H est de la forme H D I 2uu> avec kuk D 1. Il s’agit cette fois d’une symétrie par
rapport à un plan de normale u comme illustré sur la figure 2.2. En conséquence, H est sy-
métrique, orthogonale, et conserve les valeurs propres. Ces transformations de Householder

u
a1

a1
e1

Figure 2.2 — Transformation de Householder de a1

sont appliquées sur les colonnes de A D Œa1 a2


an  successivement. Par exemple, sur la
première, on veut Ha1 D ˛e1 d’où le choix :

˛ D ka1 k D kHa1 k > 0 (2.27)

alors a1 2u.u> a1 / D ˛e1 qui implique a1 ˛e1 D 2u.u> a1 /. En prenant le produit

scalaire par a1 , on obtient :

˛2 ˛a> >
1 e1 D 2.u a1 /
2
(2.28)

d’où :
r
> 1
u a1 D ˛.˛ 1 e1 /
a> (2.29)
2

On peut noter que le terme .˛ a> 1 e1 / est non nul si a1 n’est pas déjà parallèle à e1 .
Finalement, on en déduit le vecteur cherché :
a1 ˛e1
uD (2.30)
2u> a1
 41 Chapitre 2 | Systèmes aux valeurs propres

Le procédé se poursuit en considérant ensuite la deuxième colonne de A et H telle que :

 
I 0
HD (2.31)
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 remplissage de A lors des transformations mais c’est
assez délicat. En général, on ne stocke pas la matrice Q ; l’application de Q à un vecteur
est coûteuse (4n opérations), cette méthode n’est donc pas appliquée à la résolution de
systèmes linéaires. Cet algorithme se trouve aussi sous le nom de Householder-Businger-
Golub. La méthode est présentée par l’algorithme 2.4. En ce qui concerne le stockage, vk

boucle sur k D 1; 2; : : : ; n 1
ı A.k/> A.k/
pkWn;k kWn;k
˛ ı
1
ˇ .k/
ı ˛Ak;k
.k/ .k/
vkWn AkWn;k
.k/ .k/
vk vk ˛
.k/> .k/
zkC1Wn
>
ˇvkWn AkWn;kC1Wn
.kC1/
AkWn;kC1Wn A.k/
kWn;kC1Wn
v.k/ z>
kWn kC1Wn
Rk;k ˛
Rk;kC1Wn A.kC1/
k;kC1Wn
fin

Algorithme 2.4 — Transformation de Householder pour la factorisation QR

écrase A1Wk;k en mémoire, RkC1;kC1Wn écrase AkC1;kC1Wn . Rk;k doit être stocké ailleurs
pour chaque itération, d’où un 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 2n3 =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 D x avec des valeurs propres simples. On procède alors de la façon suivante :
1. A.1/ D A ;
2. factorisation QR de A.k/ D Q.k/ R.k/ ;
3. calcul de A.kC1/ D R.k/ Q.k/ ;
4. itération suivante.
 Section 2.6 | Transformation de Householder 42

Propriétés
— si les valeurs propres sont toutes simples, A.k/ tend vers la matrice diagonale des
valeurs propres. Q.k/ tend vers la matrice des vecteurs propres ;
— 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 de l’ordre de 2n3 =3 opérations et un produit
QR, de l’ordre de n3 =2 opérations.
On peut aussi appliquer des translations (décalages) :

A.k/ ˛ .k/ I D Q.k/ R.k/

(2.32)
A.kC1/ D R.k/ Q.k/ C ˛ .k/ I

Les matrices restent toujours semblables, donc conservent les valeurs propres. Concernant
.k/
la convergence de la méthode, avec des décalages, et en prenant ˛ .k/ D An;n :

jn A.k/
n;n j ! 0 (2.33)
k!1

alors :

k.A kC1 I/q.k/ k D O.k.A k I/q.k 1/ 3

k / (2.34)

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 [31] et est dédiée aux matrices A réelles symétriques. On cherche

0
ak

Figure 2.3 — Matrice en cours de tri-diagonalisation

cette fois à rendre A tri-diagonale. Il s’agit de la même technique colonne par colonne que
précédemment. On prend ici une transformation qui ne s’applique qu’à la deuxième colonne
pour la première étape :
 
I 0
HD Q (2.35)
0 H
43 Chapitre 2 | Systèmes aux valeurs propres

mais comme A est symétrique, il faut appliquer les transformations à gauche et à droite :
1
H AH D H> AH D HAH (2.36)
A.k/ prend donc la forme de la figure 2.3 et le principe consiste à travailler sur le vecteur
.k/ .k/
ŒAk;kC1


Ak;n > pour le rendre colinéaire à ekC1 . Cette méthode a la propriété intéres-
sante 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é :
 
A1;1 aQ >
A.1/ D 1
Q .1/ (2.37)
aQ 1 A
on prend :
Q .1/ D 1
H 2uQ uQ > (2.38)
en souhaitant avoir :
Q .1/ aQ 1 D ˛Qe1
H (2.39)
En fait, ici eQ 1 D e2 et on procède comme précédemment pour avoir u.
Q On cherche donc

boucle sur k D 1; 2; : : : ; n 1
.k/> .k/
ı AkC1Wn;k AkC1Wn;k
p
˛ ı
ˇ 1=.ı ˛A.k/kC1;k /
v.k/
kC1Wn
A.k/
kC1Wn;k
.k/ .k/
vkC1 vkC1 ˛
zkC1Wn ˇA.k/ v.k/
kC1Wn;kC1Wn kC1Wn
ˇ=2v.k/> z
kC1Wn kC1Wn
zkC1Wn zkC1Wn v.k/
kC1Wn
A.kC1/
kC1Wn;kC1Wn
A.k/
kC1Wn;kC1Wn
v.k/ z>
kC1Wn kC1Wn
zkC1Wn v.k/>
kC1Wn
fin

Algorithme 2.5 — Transformation de Householder

Q .1/ D I ˇ yQ yQ > avec ˛ D aQ > aQ 1 , ˇ D 1=.˛.˛

H A2;1 // et yQ D aQ 1 ˛Qe1 . On calcule
1
ensuite la transformée par :
H Q .1/ D .1 ˇ yQ yQ > /A
Q .1/ H
Q .1/ A Q .1/ .1 ˇ yQ yQ > /
(2.40)
DA Q .1/ ˇ yQ .Qy> AQ .1/ / ˇ.A Q .1/ yQ /Qy> C ˇ 2 yQ .Qy> A
Q .1/ yQ /Qy>

Q .1/ yQ et Qt D zQ ˇ >
et en utilisant les intermédiaires zQ D ˇ A y zQ /Qy
2 .Q :
Q .1/ A
H Q .1/ D A
Q .1/ H Q .1/ yQ zQ > zQ yQ > C ˇ.Qz> yQ /QyyQ >
(2.41)
DAQ .1/ yQ Qt> QtyQ >
C’est la technique mise en œuvre dans l’algorithme 2.5, 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.
 Section 2.7 | Exemple d’utilisation dans un code éléments finis 44

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

On considère Kx D 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 > m,
défini par p vecteurs rangés dans la matrice de base V D Œv1


vp . Éventuellement, les
vecteurs vi sont orthonormés au sens de K, soit V> KV D K O D I.
On cherche alors les vecteurs propres x sous la forme x D VOx, xO est de dimension p. Le
O
problème est alors KVOx D MVO x, d’où :
O > MVOx
V> KVOx D V , O x D O MO
KO Ox (2.42)

Le problème (2.42) a été construit par une réduction de Ritz. Sa résolution complète permet
O xO /. On retourne alors dans l’espace de départ où
de trouver p caractéristiques propres .;
x D VOx est l’approximation d’un vecteur propre, O 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.7.2 Technique itérative

Rapidement, il s’agit de :
1. choisir p vecteurs de départ ;
O x D O MO
2. réduire à KO O 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 D min.2m; m C 8/ (2.43)

La convergence pour une valeur propre i est proportionnelle à ji j=jpC1 j.

2.8 Méthode de Lanczos

Cette méthode est due à Lanczos [44] en 1950. On travaille sur un sous-espace généré par
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’orthogonalité sur les vecteurs successifs (elle se
produit à chaque fois qu’une valeur propre converge), d’où des instabilités numériques. En
45 Chapitre 2 | Systèmes aux valeurs propres

v.0/ 0, ˇ .0/ 0
Vecteur de départ v.1/ , kv.1/ k D 1
tant que R > ", itérer sur k
résoudre Az D v.k/
˛ .k/ v.k/ > z
d z ˛v p
.k/
ˇ .k 1/ v.k 1/
ˇ .k/ 1= d> d
v .kC1/
ˇ .k/ d
recherche de .Ox; / O tel que TOx D O
Ox
fin

Algorithme 2.6 — Méthode de Lanczos

pratique, pour pouvoir traiter de grands problèmes, il est nécessaire d’utiliser des stratégies
de ré-orthogonalisation. Ces dernières ne seront pourtant pas décrites ici.
En termes simples, l’algorithme se résume à choisir un vecteur de départ v.0/ puis de
résoudre le système Az.k/ D v.k/ , d’orthogonaliser z.k/ par rapport aux v.j / , j 6 k puis
de normer v.kC1/ z.k/ =kz.k/ k.
Pour démontrer le taux de convergence, il suffit de chercher v.kC1/ sous la forme :

ˇ .k/ v.kC1/ D A 1 .k/

v ˛ .k/ v.k/ ˇ .k 1/ .k 1/
v (2.44)

ce que l’on peut montrer par récurrence, avec v.k/> v.j / D 0 pour j D k 1 et j D k 2,
et v.k/> v.k/ D 1, pour avoir orthogonalité par rapport à tous les v.j / , j < k. Ces conditions
donnent les valeurs de ˛ .k/ D v.k/ A 1 v.k/ et ˇ .k 1/ D v.k/ A 1 v.k 1/ . Matriciellement,
en notant V de taille .n; k/ la matrice qui contient tous les v.j / disponibles à l’itération k,
la relation (2.44) est équivalente à :
2 .0/ 3
˛ ˇ .0/ 0
6 .0/ : : :: 7
1
6ˇ : : 7  
A V V6 6
:: ::
7 D 0


0 ˇ .k/ v .kC1/
7 (2.45)
4 : : ˇ .k 1/ 5
0 ˇ .k 1/ ˛ .k/
„ ƒ‚ …
T

ce qui implique :

V> A 1
V TD0 (2.46)

Par conséquent, T a les mêmes valeurs propres que V> A 1 V (voir le principe des sous-
espaces). On cherche ensuite les caractéristiques propres de T, TOx D O O x, O sont des ap-
proximations des valeurs propres de A et VOx, celles des vecteurs propres associés.
1

Afin d’illustrer le comportement de cette méthode, la figure 2.4 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
valeurs propres très proches (les cercles sont les valeurs propres exactes). On voit que la
multiplicité et/ou proximité des valeurs propres échappe à l’algorithme, et que les valeurs
propres se mettent à converger quasiment l’une après l’autre.
 Section 2.8 | Méthode de Lanczos 46

12

10

valeurs propres 6

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

Figure 2.4 — 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 D A 1 VOx O x avec :
VO
 
A 1 VOx D VTOx C 0


0 ˇ .k/ v .kC1/ xO (2.47)

O x, d’où r D ˇ .k/ v.kC1/ xO k .

avec TOx D O
L’algorithme 2.6 décrit cette méthode, pour lequel les ˛ .k/ et les ˇ .k/ sont les compo-
santes 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 parti-
culiers, 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ésenta-
tion dédiée aux éléments finis, le lecteur est invité à se reporter, par exemple, aux nombreux
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 [41], décrites ici.
Pour plus d’informations, on pourra consulter [11].

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é)
" D ".U/ (3.1)
— les équations d’admissibilité statique : le champ de contraintes  doit être en équi-
libre avec les charges extérieures ;
— les relations de comportement [45]. 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 :
" D "e C "p (3.2)
Les lois d’état reliant les variables cinématiques et statiques sont associées à l’éner-
gie libre (qui dépend de l’état) ."e ; X/ :
@
D
@"e (3.3)
@
YD
@X
 Section 3.1 | Non-linéarités matérielles de type plasticité 48

L’évolution de l’état, pour un modèle standard associé, se déduit du potentiel de

dissipation ˆ.; Y/ par la loi de normalité :

@ˆ
"Pp D P
@ (3.4)
@ˆ
P P
XD
@Y

avec, dans le cas non visqueux, ˆ 6 0, P > 0 et ˆP D 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/ D K.u/u fD0 (3.5)

r.u/ est le gradient à annuler. La matrice de rigidité dépend de la solution u, par l’in-
termédiaire de la relation de comportement du matériau. Les méthodes de type Newton
sont habituellement utilisées : si u0 est une approximation de la solution, cette dernière est
u D u0 C ıu et un développement de Taylor de r donne au voisinage de u :

@r
.u0 /ıu D r.u0 / (3.6)
@u

ce qui suggère un algorithme itératif de la forme :

uiC1 ui D si Hi ri (3.7)

où ri D [Link] /, Hi est une approximation de l’inverse de la matrice jacobienne [Link] / D

@r
.u /. si est une amplitude destinée à réduire riC1 1 .
@u i
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.

3.1.3 Méthode de la plus grande pente

Dans ce cas :

HDI (3.8)

1. Dans un premier temps, on peut choisir si D 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
 49 Chapitre 3 | Résolution de systèmes non-linéaires particuliers

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 D 1
di D (3.9)
ri C i di 1 pour i > 1

et :

uiC1 D ui si di (3.10)

i est une scalaire utilisé pour conjuguer les directions de recherche, i.e. d>
i di 1 D 0, pour
lequel deux expressions classiques sont :
1. Fletcher-Reeves [21] :

r>
i ri
i D (3.11)
r>
i 1 ri 1

2. Polak-Ribière [57] :

r>
i .ri ri 1 /
i D (3.12)
ri 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 [14].

3.1.5 Méthode de Newton-Raphson

Ici, si D 1 et :
" # 1 " # 1
1
@r @K
Hi D [Link] / D .ui / D [Link] / C .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 D J.u? / 1
(3.14)
 Section 3.1 | Non-linéarités matérielles de type plasticité 50

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.
HiC1 D Hi C ıHi (3.15)
en s’assurant que la condition sécante :
HiC1 .uiC1 ui / D riC1 ri (3.16)
est satisfaite. Elle correspond à une approximation au premier ordre de la matrice jacobienne
au voisinage de ui .
Si Hi est mis à jour, on parle de mise à jour inverse, alors qu’une mise à jour directe
modifie Hi 1 .
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 D ab> où a et b sont deux vecteurs, la formule de Sherman-
Morrison-Woodbury [18, 30, 35] permet d’exprimer directement l’inverse, tant que H n’est
pas singulier et que 1 C b> H 1 a ¤ 0 :
H 1 ab> H 1
.H C ab> / 1
DH 1
(3.17)
1 C b> H 1 a
L’algorithme générique consiste à itérer sur :
— direction de recherche di D Hi ri ;
— minsi >0 kriC1 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éra-
tions, en utilisant l’information des gradients successifs. Par exemple, avec la précédente
correction de rang 1, en vérifiant la symétrie et la condition sécante, on obtient :
.ı i Hi i /.ı i H i i />
HiC1 D Hi C (3.18)
>
i .ı i Hi i /
où ı i D uiC1 ui et i D riC1 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 d>
i 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 [20] préserve le caractère défini positif. En démar-
rant avec H0 symétrique, défini positif, la mise à jour est :
ı i ı> .Hi i /.Hi i />
HDFP
iC1 D Hi C
i
(3.19)
ı>
i ıi i 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 D J 1 . Si de plus H0 D I, il s’agit de la méthode du gradient
conjugué.
51 Chapitre 3 | Résolution de systèmes non-linéaires particuliers

Méthode de Broyden-Fletcher-Goldfarb-Shanno
Cette méthode de rang 2, proposée dans [9], préserve aussi le caractère défini positif. En
démarrant avec H0 symétrique, défini positif, la mise à jour est :
! !>
ıi Hi i ıi Hi i
HBFGS DFP >
iC1 D HiC1 C . i Hi i /
ı>
i i i Hi i
> ı>
i i i Hi i
>
! (3.20)
> >
i Hi i ı i ı i
> i Hi C Hi i ı i
ıi >
D Hi C 1 C
ı>
i i ı>
i i ı>
i 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 :

HiC1 D ˛HBFGS
iC1 C .1 ˛/HDFP
iC1 (3.21)

˛ étant tel que le caractère défini positif est assuré (0 6 ˛ 6 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 à :

HiC1 D Vi H0 (3.22)

avec :

j ı>
j
Vj D I (3.23)
ı>
j j

et la méthode BFGS, à :
i
!> i
! i i
!> i
!
Y Y X Y ık ı >
k
Y
HiC1 D Vj H0 Vj C Vj Vj (3.24)
j D0 j D0 j D0 kDj C1
ı>
k k kDj C1

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
limited memory de BFGS (L-BFGS) [42, 47].
4 Résolution de systèmes
différentiels

Il est possible de consulter sur le sujet les références [36, 40]. 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 2 Œ0; T, telle que :
dx
D f .x; t/
dt (4.1)
x.0/ D u0

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 :

9k > 0; 8t 2 Œ0; T; 8.y; z/ 2 R2 ; jf .z; t/ f .y; t/j 6 kjz yj (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 :
Z t
x D u0 C f .x; t/ dt (4.3)
0

Une discrétisation de l’intervalle de temps Œ0; T est la séquence d’instants 0 D t0 < t1 <
: : : < tm D T. Le pas de temps est h D tiC1 ti . Le principe est de chercher à approcher la
valeur de la fonction aux piquets de temps, [Link] /, par la valeur xi , ce qui permet de définir
l’erreur ei D [Link] / xi j. Dans la suite, on notera fi D f .xi ; t/. Le problème incrémental
est, ayant la solution approchée jusqu’à ti , de la déterminer à tiC1 avec :
Z t
x D xi C f .x; t/ dt (4.4)
ti

Exemple 4.1 — Schéma d’Euler explicite. La figure 4.3(a) illustre l’estimation de l’inté-

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

 Section 4.1 | Équation scalaire du premier ordre 54

grale cherchée. Le schéma s’écrit :

xi C1 D xi C hfi
(4.5)
x0 D u0

4.1.3 Schémas explicites à un pas

Ce sont les schémas de la forme :
xiC1 D xi C hˆ.xi ; ti ; h/
(4.6)
x0 D u0

Convergence Si on a la propriété supi jei j !h!0 0 (en l’absence d’erreur d’arrondi !)

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

xQ iC1 D xQ i C hˆ.xQ i ; ti ; h/ C "iC1

(4.7)
xQ 0 D u0 C "0

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

m
X
9S > 0; sup jxQ i xi j 6 S j"i j (4.8)
i iD0

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


ciC1 D x.tiC1 / [Link] / C hˆ.[Link] /; ti ; h/ (4.9)

Le schéma est consistant ssi :

X
jci j ! 0 (4.10)
h!0
i

Ordre On dit d’un schéma qu’il est d’ordre p quand 9M > 0, supi jci j 6 MhpC1 .
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é

Figure 4.1 — Propriétés de base d’un schéma d’intégration

 55 Chapitre 4 | Résolution de systèmes différentiels

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/ D f .x; t/ , d’ordre p 6 1
— d’ordre p ) jei j 6 Nhp

Exemple 4.2 — Schéma d’Euler explicite. Dans ce cas, ˆ.x; t; h/ D f .x; t/ et par consé-
quent, le schéma est consistant et stable, donc convergent :


ci C1 D [Link] C1 / [Link] / C hf .[Link] /; ti / D O.h2 / (4.11)

P i / C O.h2 / et x.t
or [Link] C1 / D [Link] / C hx.t P i / D f .[Link] /; ti / 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 :

xQ iC1 D xQ i C hˆ.xQ i ; ti ; h/ C "iC1

(4.12)
xQ 0 D u0 C "0

Si le schéma est convergent, ei D [Link] / xi j ! 0. Qu’en est-il de eQi D [Link] / xQ i j ?

h!0

eQi D [Link] / xi C xi xQ i j 6 [Link] / xi j C jxi xQ i j (4.13)

„ ƒ‚ … „ ƒ‚ …
① ②

avec, dans l’équation (4.13) :

— le schéma est d’ordre p donc le terme ① 6 Nhp ;
P Pm 
m
— le schéma est stable donc le terme ② 6 S iD0 j"j C iD1 hj"j .
Au bilan :

eQi 6 SŒ.m C 1/j"j C hmj"j C Nhp (4.14)

avec hm D T, d’où :
 T
eQi 6 Sj"j 1 C T C C Nhp (4.15)
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.
 Section 4.1 | Équation scalaire du premier ordre 56

Problème numériquement mal posé

Considérons le problème suivant :
dx
D 4x 1
dt (4.16)
x.0/ D 1=4
dont la solution est x.t/ D 1=4. Perturbons maintenant les conditions initiales, soit x.0/ D
1=4 C ". La nouvelle solution x.t/ D 1=4 C "e4t est alors très différente de la précédente.
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
D ! 2x
dt (4.17)
x.0/ D 1
dont la solution exacte est :
!2t
x.t/ D e (4.18)
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 (ou A-stable pour abso-
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 D 1=10 (b) h D 1=7

Figure 4.2 — Stabilité du schéma d’Euler explicite

lute stability) [13] se définit pour les équations linéaires. On teste alors le schéma sur un
problème du type :
dx
D cx (4.19)
dt
2. On peut d’ailleurs se poser des questions quant à la notion même de solution pour un tel problème.
 57 Chapitre 4 | Résolution de systèmes différentiels

où c est une constante complexe et pour lequel où le schéma est tel que xiC1 D [Link]/xi ,
où r est un polynôme ou une fraction rationnelle. On définit alors le domaine de A-stabilité
comme suit :

Dstab D fz 2 C; jr.z/j 6 1g (4.20)

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

fz 2 C; <.z/ 6 0g  Dstab (4.21)

Exemple 4.3 — Schéma d’Euler explicite. La relation est la suivante :

xi C1 D xi C hcxi D .1 C hc/xi ) r.z/ D 1 C z (4.22)

La condition de stabilité est donc :

j1 C hcj 6 1 (4.23)

Pour c D ! 2 2 R , elle s’écrit h 6 2=! 2 . Il s’agit donc d’un schéma conditionnellement stable.
Pour l’exemple précédent, avec c D 15 (la solution tend vers 0), le critère de stabilité donne
0 6 h 6 1=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 approximation bornée.

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

Cette méthode est décrite avec un paramètre ˛ 2 Œ0; 1 par :

xiC1 D xi C h Œ.1 ˛/fi C ˛fiC1 

(4.24)
x0 D u0

Le calcul de xiC1 fait intervenir ˛fiC1 D ˛f .xiC1 ; tiC1 /. Il s’agit donc d’une méthode
implicite si ˛ ¤ 0.
Le cas ˛ D 0 conduit à la méthode d’Euler explicite ; dans le cas où ˛ D 1, il s’agit
de la méthode d’Euler implicite ; enfin si ˛ D 1=2, 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.

Stabilité
On considère le problème :

dx
D cx (4.25)
dt
le schéma donne :

xiC1 D xi C hŒ.1 ˛/cxi C ˛cxiC1  (4.26)

 Section 4.1 | Équation scalaire du premier ordre 58

f f f

t t t
ti ti C1 ti ti C1 ti ti C1

(a) ˛ D 0 (b) ˛ D 1 (c) ˛ D 0;5

Figure 4.3 — Différents types d’intégration

donc :

1
xiC1 D .1 ˛hc/ .1 C .1 ˛/hc/ xi (4.27)
„ ƒ‚ …
B

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 :

j1 C .1 ˛/hcj
61 (4.28)
j1 ˛hcj

— si ˛ > 1=2 le schéma est inconditionnellement stable ;

— si ˛ < 1=2 avec c D ! 2 , la condition de stabilité s’écrit h 6 2=..1 2˛/! 2 /.
Les différents domaines de stabilité correspondants sont tracés sur la figure 4.4.

=.hc/
˛ D 1=2

˛ > 1=2
˛ < 1=2 ˛D
1
˛D0
<.hc/
1 1

Figure 4.4 — Domaines de stabilité pour la méthode des trapèzes généralisés

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 ˛ ¤ 1=2, et
qu’il est d’ordre 2 pour ˛ D 1=2, c’est-à-dire pour Crank-Nicolson.
 59 Chapitre 4 | Résolution de systèmes différentiels

4.1.7 Méthode à pas multiples : Adams

L’idée des méthodes à pas multiple est d’utiliser non seulement l’évaluation fi mais aussi
un certain nombre de valeurs passées fi 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.

Adams explicite
On utilise pour cela un développement de Taylor (avant) de x.tiC1 / :

h2
x.tiC1 / D [Link] / C hx.t
P i/ C R i/ C : : :
x.t (4.29)
2Š

avec x.t R i / D df
P i / D f .[Link] /; ti / et x.t dt .[Link] /; ti /. On garde ensuite les termes d’un ordre
suffisant 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 xiC1 D xi Chfi CO.h2 /, c’est Euler explicite !
Pour Adams explicite d’ordre 2, on conserve les termes en O.h2 / qui impliquent la
dérivée première de f ; pour l’estimer on utilise encore le développement d’une fonction g :

g.tiC1 / D [Link] / P i / C O.h2 /

hg.t (4.30)

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

1
P i/ D
g.t Œ[Link] / g.tiC1 / C O.h/ (4.31)
h
df
Avec ce procédé appliqué à dt , on obtient la méthode cherchée :

h
xiC1 D xi C Œ3fi fi 1 C O.h3 / (4.32)
2
La démarche se prolonge pour les ordres de schéma supérieurs. On obtient ainsi les formules
d’Adams-Bashforth :

xi D xi 1 C hfi 1
h
xi D xi 1 C .3fi 1 fi 2 /
2
h
xi D xi 1C .23fi 1 16fi 2 C 5fi 3 /
12
h
xi D xi 1C .55fi 1 59fi 2 C 37fi 3 9fi 4 / (4.33)
24
h
xi D xi 1C .1901fi 1 2774fi 2 C 2616fi 3 1274fi 4 C 251fi 5 /
720
h
xi D xi 1C .4277fi 1 7923fi 2 C 9982fi 3 7298fi 4 C 2877fi 5 475fi 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.
lac
Section 4.1 | Équation scalaire du premier ordre 60

3
1;5
2
1 kD4

0;5 1
kD3
kD2 kD3 kD4
=.z/

kD1 kD2

=.z/
0 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

Figure 4.5 — Domaine de stabilité dans le plan complexe

Adams implicite
L’idée est la même que précédemment, en utilisant un développement de Taylor arrière de :
h2
[Link] / D x.tiC1 h/ D x.tiC1 / P iC1 / C
hx.t R iC1 / C : : :
x.t (4.34)
2Š
La méthode d’Adams implicite d’ordre 1 est ainsi :
xi D xiC1 hfiC1 C O.h2 / (4.35)
d’où :
xiC1 D xi C hfiC1 C O.h2 / (4.36)
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 :
xi D xi C hfi
1
h
xi D xi 1 C .fi C fi 1 /
2
h
xi D xi 1 C .fi C 8fi 1 fi 2 /
12
h
xi D xi 1 C .9fi C 19fi 1 5fi 2 C fi 3 / (4.37)
24
h
xi D xi 1 C .251fi C 646fi 1 264fi 2 C 106fi 3 19fi 4 /
720
h
xi D xi 1 C .475fi C 1427fi 1 798fi 2 C 482fi 3 173fi 4 C 27fi 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.
 61 Chapitre 4 | Résolution de systèmes différentiels

Amorçage du schéma
Aux premiers pas de temps, on ne dispose pas forcément de tous les fi 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 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 xiC1 ,
.0/
on évalue f .xiC1 ; tiC1 / qui permet, en utilisant la formule d’intégration d’Adams impli-
.1/
cite, d’atteindre une valeur xiC1 . Cette dernière étape est appelée correction. On itère ainsi
.0/
jusqu’à une convergence souhaitée. Pour avoir la première valeur xiC1 , 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
xiC1 xi C .3fi fi 1/ (4.38)
2
2. modificateur :
.0/ .0/ 1 .0/

xO iC1 xiC1 xi xi (4.39)
12
3. correcteur (AI2) :
h O

xiC1 xi C fiC1 C fi 1 (4.40)
2
.0/
où fOiC1 D f .xO iC1 ; tiC1 /. Le coefficient 1
12 est celui issu de la formule de (AI3) ? ? ?.

4.1.9 Méthode de Runge-Kutta

Historiquement présentée dans [43] en 1901, cette méthode est de la forme :
xiC1 D xi C hˆ.x; t; h/
(4.41)
x0 D u0
avec :
q
X
ˆ.x; t; h/ D al kl .x; t; h/
lD1
k1 .x; t; h/ D f .x; t/ (4.42)
l 1
X .l/

kl .x; t; h/ D f x C ˇj kj .x; t; h/; t C ˛l h
j D1
 Section 4.1 | Équation scalaire du premier ordre 62

L’idée est d’intégrer en utilisant des valeurs en des points particuliers, avec des poids d’in-
.l/
tégration particuliers. Les constantes al , ˛l et ˇj sont ajustées pour obtenir l’ordre sou-
haité.
Exemple 4.4 — Schéma à l’ordre 2. On prend ˆ D af .x; t/ C bf .t C ˛h; x C ˇhf .x; t//
et on détermine a, b, ˛, ˇ pour avoir :
1


x.t C h/ x.t/ ˆ x.t/; t; h D O.h2 / (4.43)
h
„ ƒ‚ … „ ƒ‚ …
① ②
avec, dans l’équation (4.43) :

h  @f @f 
①Df C C f C O.h2 /
2 @t @x
 @f @f 
② D af C b f C ˛h C ˇhf C O.h2 /
@t @x
Il faut donc avoir :
h h
0D1 a b 0D b˛h 0D bˇh (4.44)
2 2
1
d’où les coefficients cherchés a D 1 b, ˛ D ˇ D 2b et le schéma :
 h h
ˆ.x; t; h/ D .1 b/f .x; t/ C bf x C f .x; t/; t C (4.45)
2b 2b
où b est une constante arbitraire non nulle. Des cas particuliers existent :
— b D 1 : il s’agit de la méthode d’Euler modifiée, d’ordre 2 ;
— b D 21 : il s’agit de la méthode de Heun, d’ordre 2.

Exemple 4.5 — Schéma à l’ordre 4. Bien entendu, cette technique se poursuit pour les ordres
supérieurs. Un ordre souvent utilisé avec cette méthode est l’ordre 4. Le schéma se présente alors
comme suit :
1
ˆ.x; t; h/ D .k1 C 2k2 C 2k3 C k4 /
6
k1 D f .x; t/
 h h
k2 D f x C k1 ; t C (4.46)
2 2
 h h
k3 D f x C k2 ; t C
2 2
k4 D f .x C hk3 ; t C h/

Bien entendu, des méthodes de prédiction-correction basées sur la méthode de Runge-Kutta peuvent
aussi être mises en place.
 63 Chapitre 4 | Résolution de systèmes différentiels

4 3 2 1 0 1

Figure 4.6 — Domaine de stabilité pour Runge-Kutta (ordre 1 à 8)

Cette méthode requière beaucoup d’évaluations 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 [53].
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 D 0 (4.47)
0 dt
L’idée est de pouvoir utiliser des discrétisations en temps qui tolèrent a priori des disconti-
nuités aux pas de temps : à t D ti , on peut définir une valeur à gauche xi et une valeur à
droite xiC différentes. Le saut est noté JxKi D xiC xi . L’intégrale précédente est donc à
prendre en un sens généralisé, qui fait « travailler » les sauts 3 :
Z T X Z ti C1  
dx 
? dx 
? ?
f .x; t/ x dt D f .x; t/ x dt C JxKi xi D 0 (4.48)
0 dt ti dt
i

La fonction test est choisie dans le même espace que x, et on définit sa valeur en ti comme
xi? D ˛xi?C C .1 ˛/xi? . Dans le cas standard, on prend ˛ D 1 4 . Le problème peut alors
3. On suppose cependant f .x; t/ suffisamment régulière.
4. one-sided upwinding
 Section 4.2 | Systèmes différentiels du premier ordre 64

s’écrire sur chaque intervalle de temps :

Z ti C1 
dx 
f .x; t/ x ? dt C JxKi xi? D 0 (4.49)
ti dt

Cette méthode est très ressemblante à une méthode de Runge-Kutta à un pas implicite (A-
stable) d’ordre 2l C 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 :

MPx C Kx D g.t/
(4.50)
x.0/ D 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 4.6 — Schéma d’Euler explicite. Le problème précédent peut s’écrire :

1 1
xP D M Kx C M g (4.51)

On applique alors le schéma voulu :

1 1
xi C1 D xi C h. M Kxi C M gi / (4.52)

Pour éviter de faire intervenir directement M 1

, on préfère le mettre sous la forme :

Mxi C1 D Mxi C [Link] C gi /

(4.53)
x0 D 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 pulsation ! la plus grande. Dans le cas d’Euler explicite, on a donc le critère de
stabilité h 6 2=!M 2 et on peut borner la pulsation propre maxi du système par la pulsation

propre maxi élémentaire !M 6 !max e , voir annexe B.

 65 Chapitre 4 | Résolution de systèmes différentiels

Pour avoir un vrai schéma explicite, il faut ne pas avoir à résoudre de système linéaire
couplé dans (4.53). 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 [10], on peut se poser la question de la conservation des
propriétés du schéma.

4.2.2 Influence de l’ordre d’intégration

Si on considère que le coût est directement lié au nombre d’évaluations de la fonction f ,
on peut tracer le coût en fonction de l’erreur obtenue pour différents schémas. Choisissons
le cas test suivant :
       
xP y x.0/ u0
D avec D (4.54)
yP .1 x 2 /y x y.0/ 0

Il s’agit d’un système différentiel du premier ordre, non linéaire, pour t 2 Œ0; T. Si ˛ 
2;008619861986087484313650940188 et avec T  6;6632868593231301896996820305,
la solution est périodique, de période T. Une mesure de l’erreur peut donc être de prendre à

106
Euler explicite

Heun
104

Runge-Kutta 4

102
10 8 10 6 10 4 10 2

Figure 4.7 — Influence de l’ordre d’un schéma sur le coût

la fin de l’intervalle d’étude :

   
x x0
eDk T k (4.55)
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 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.
 Section 4.3 | Systèmes différentiels du deuxième ordre 66

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

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

MRx C CPx C Kx D g.t/

xP .0/ D v0 (4.56)
x.0/ D 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 D 0. L’équilibre
est donc :

MRx C Kx D g.t/ (4.57)

4.3.1 Méthode de Newmark

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

xiC1 D xi C hPxi C h2 .1 2ˇ/Rxi C 2ˇ xR iC1
2
(4.58)
xP iC1 D xP i C h .1 /Rxi C xR iC1

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

membre. On peut aussi remarquer que si ˇ D D 0, on retrouve des développements de
Taylor puisque le schéma de Newmark ressemble à un développement au troisième ordre :

h2 h3 xR iC1 xR i
xiC1 D xi C hPxi C xR i C 6ˇ
2 6 h (4.59)
h2 xR iC1 xR i
xP iC1 D xP i C hRxi C 2
2 h
Il reste à écrire l’équilibre (4.57) à l’instant tiC1 avec le schéma précédent. On peut alors
prendre l’accélération comme inconnue :

 1  
ˇh2 K C M xR iC1 D giC1 K xi C hPxi C h2 ˇ xR i (4.60)
2
et on obtient le système à résoudre pour atteindre xR iC1 ; mais on peut aussi l’écrire avec le
déplacement comme inconnue :

 1  
ˇh2 K C M xiC1 D ˇh2 giC1 C M xi C hPxi C h2 ˇ xR i (4.61)
2
qui est le système à résoudre pour atteindre xiC1 . Dans ce dernier cas, cependant, un pro-
blème subsiste si ˇ D 0. En effet, soit on utilise l’équilibre pour revenir à l’accélération :

1
xR iC1 D M .giC1 KxiC1 / (4.62)
67 Chapitre 4 | Résolution de systèmes différentiels

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

1 1 1 2ˇ
xR iC1 D .xiC1 xi / xP i xR i (4.63)
ˇh2 ˇh 2ˇ
mais on s’aperçoit bien du problème quand ˇ D 0.

Cas particuliers
1
— D 2, ˇ D 14 : méthode des accélérations moyennes ;
1
— D 2, ˇ D 16 : méthode des accélérations linéaires ;
— D 1
2, ˇ D 0 : méthode des différences centrées, qui peut être explicite 5 ;
1 1
— D 2, ˇ D 12 : méthode de Fox-Goodwin.

Ordres
1
On peut montrer qu’il s’agit d’une méthode d’ordre 1 pour ¤ 2 et d’ordre 2 pour D 21 .

Stabilité
Dans le cas d’une seule équation scalaire, xR D ! 2 x, le schéma donne :
     
1 C .h!/2 ˇ 0 xiC1 1 .h!/2 . 21 ˇ/ h xi
D (4.64)
h! 2 ˛ 1 xP iC1 h! 2 .1 ˛/ 1 xP i
„ ƒ‚ … „ ƒ‚ … „ ƒ‚ … „ƒ‚…
A1 yi C1 A2 yi

d’où la relation yiC1 D A1 1 A2 yi et la matrice d’itération B D A1 1 A2 . La condition de

stabilité est donc %.B/ < 1 ou bien %.B/ 6 1 avec des valeurs propres simples. L’équation
caractéristique permettant de déterminer les valeurs propres de B est 2 2b C c D 0
avec :
 1  .h!/2  1  .h!/2
2b D 2 ˛C c D 1 ˛ (4.65)
2 1 C ˇ.h!/2 2 1 C ˇ.h!/2
les conditions de stabilité sont donc :
c61 1 2b C c > 0 1 C 2b C c > 0 (4.66)
sauf en ce qui concerne les points isolés .b; c/ D . 1; 1/ et .1; 1/. D’où les conditions en
terme de coefficients de la méthode :
1 h!.2ˇ /
> 2C >0 (4.67)
2 1 C ˇ.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
h6 p (4.68)
! =2 ˇ
La figure 4.8 illustre ces résultats.
5. Cette méthode est aussi appelée Velocity Verlet algorithm ou Leap Frog en dynamique moléculaire !
 Section 4.3 | Systèmes différentiels du deuxième ordre 68

c ˇ

1 IS
racines complexes
conjuguées
1

1 1 0 1=2
racines
réelles

Figure 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 6 :

xR C 2! xP C ! 2 x D 0 (4.69)

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

   1=2
1 1 1 2
 2 C 2 ˇ C  2 . 2/ 1
hN 6 1
(4.70)
2 ˇ !

Dans le cas où > 21 , l’amortissement a un effet stabilisant puisque sur les conditions de
stabilité, on a hN < 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
ˇ6 C (4.71)
4 2
Le point de dissipation Haute Fréquence maximal est déterminé par le couple .ˇ D 1,
D 23 /, 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
> ˇD C (4.72)
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.
69 Chapitre 4 | 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

Figure 4.9 — Parasites hautes fréquences (d’après [46])

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

les instants ti et tiC1 [37]. L’écriture du schéma de Newmark (4.58) en différence entre ces
deux instants, en utilisant la notation suivante :

1
hxi i D .xiC1 C xi /
2 (4.73)
Œxi  D xiC1 xi

conduit à :

h2
Œxi  D h hPxi i C .2ˇ /ŒRxi 
2 (4.74)
 1
ŒPxi  D h hRxi i C h ŒRxi 
2
De même, la différence des deux équilibres (4.57) devient :

MŒRxi  C KŒPxi  Œgi  D 0 (4.75)

6. Cette écriture dans la base propre d’un système entier correspond à un modèle d’amortissement de Ray-
leigh.
 Section 4.3 | Systèmes différentiels du deuxième ordre 70

contraintes (MPa) contraintes (MPa)

200
200

100 100

0
0
0 200 400 600 800 0 200 400 600 800
temps temps
(a) g D 0;55 (b) g D 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 D 0;75 (d) g D 1;5

Figure 4.10 — Influence de l’amortissement numérique (d’après [46])

pertes énergétiques (%)

0
g D 0;55

g D 0;6

g D 0;75
0;05

g D 1;5

0;1
0 200 400 600 800
temps

Figure 4.11 — Influence de l’amortissement numérique (suite) (d’après [46])

En multipliant cette dernière relation par ŒPxi T et en utilisant la propriété hxi iT Œxi  D Œ 12 x2i ,
 71 Chapitre 4 | Résolution de systèmes différentiels

15
3

2 10

1 5
kD2 kD3 kD5
=.z/

kD4 kD6

=.z/
kD1
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/

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

avec des matrices M et K positives, on obtient le bilan :

   
1 T ? 1 T 1  1
xR i M xR i C xP i KPxi D ŒPxi T Œgi  ŒRxi T M? ŒRxi  (4.76)
2 2 h 2
où M? D M C 12 h2 .2ˇ /K. On trouve donc, dans l’ordre, un terme lié à la variation de
l’énergie cinétique, un terme lié à la variation de l’énergie interne, un terme lié au travail
des 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 D 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 [23] 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
suivant :
1
xP iC1 D .3xiC1 4xi C xi 1 /
2h (4.77)
1
xR iC1 D .3PxiC1 4Pxi C xP i 1 /
2h
De la même façon, l’équilibre (4.57) à l’instant tiC1 donne le système à résoudre pour
atteindre xR iC1 :
 4  4 1 8 2 
M C h2 K xR iC1 D giC1 K xi xi 1 C hPxi hPxi 1 (4.78)
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 [5], cette méthode estime une accélération li-
néaire entre les instants ti et tiC D ti C h,  > 0 étant un paramètre de la méthode. On a
 Section 4.4 | Intégration de modèles de comportement 72

donc :
xR iC xR i
xR .ti C / D  C xR i (4.79)
h
et par intégration, on a les expressions :
1
xP iC D xP i C [Link] C xR iC /
2 (4.80)
1 1 
xiC D xi C hPxi C  2 h2 xR i C xR iC
3 6
Cette fois-ci, l’équilibre (4.57) à l’instant tiC conduit au système en xiC :
 6   6 6 
K C 2 2 M xiC D giC C M 2 2 xi C xP i C 2Rxi (4.81)
 h  h h
—  > 1;37 : la méthode est inconditionnellement stable ;
—  D 1 : méthode des accélérations linéaires.

4.3.4 ˛-HHT
Cette méthode [33] 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 :




MRxiC1 C .1C˛/CPxiC1 ˛CPxi C .1C˛/KxiC1 ˛Kxi D .1C˛/fiC1 ˛fi (4.82)

Le démarrage est réalisé avec x0 et xP 0 données ainsi que xR 0 D M 1 .f0 CPx0 Kx0 /.
Pour le choix de paramètres de la forme ˇ D .1 ˛/2 =4 et D 1=2 ˛, le schéma est
du second ordre pour les basses fréquences, il est inconditionnellement stable pour 1=3 6
˛ 6 0, la dissipation en hautes fréquences est maximum pour ˛ D 1=3, enfin, si on
impose de plus ˛ D 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 [56].

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/ D 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 D Tr ".u/ d (4.83)


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

 73 Chapitre 4 | Résolution de systèmes différentiels

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 D." "p /
(4.84)
Y D 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
plasticité, 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 P à cause
du caractère non différentiable du potentiel de dissipation :
@f
"Pp D P (4.85)
@
P D P @f
X (4.86)
@Y
On a alors les conditions de complémentarité de Khun-Tucker :

f 6 0; P > 0; f P D 0 (4.87)

L’expression du multiplicateur est donnée dans le cas de l’écoulement plastique par la condi-
tion de cohérence :
@f @f P
fP D 0 D P C Y (4.88)
@ @Y

(dans le cas contraire P D 0 et le comportement est localement élastique). En remplaçant

P par leurs expressions issues de la dérivation des lois d’état, et en utilisant les lois
P et Y
d’évolution, on obtient l’expression du multiplicateur :
@f
D"P
P D @f
@
(4.89)
D C @f
@f
@ @
@G @f
@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) :
  "  
"Pp p
P D g ; "; P
" (4.90)
X 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
 Section 4.4 | Intégration de modèles de comportement 74

fait que l’évolution de la déformation totale est imposée. Dans le système précédent, " et "P
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
structure, 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 avec 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 cela, on prend en compte la propriété de séparation en partie élastique et plastique
des modè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 C "p D ", avec "p D 0 et X D 0, les
lois d’état (4.84) donnent sur cet incrément :
    
 D 0 "
D (4.91)
Y 0 0 0

Si, au bout de cet incrément, le nouvel état . C ; "p ; X; Y C Y/ vérifie f 6 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 C "p D 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 D 0 à la fin de la correction. Cette phase de correction peut être réalisée par uti-
lisation des méthodes d’intégration déjà citées. On peut aussi utiliser des méthodes
comme celle du retour radial d’Ortiz et Simo [52, 53].

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 . C ; Y C Y/ D 0 (4.92)

75 Chapitre 4 | Résolution de systèmes différentiels

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 . C ı; Y C ıY/  f .; Y/ C ı C ıY (4.93)
@ @Y

pour lesquels on cherche à rester sur le seuil :

@f @f
f C ı C ıY D 0 (4.94)
@ @Y

Comme ı" D 0 lors de la correction, les équations (4.85) et (4.86) donnent alors la valeur
d’un « multiplicateur plastique » pour le sous-pas :

f
ı D @f
(4.95)
D @f
@ @
C @f @G @f
@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/.

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 D .; Y; "p ; X/ (4.96)

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 C s D . C ; Y C Y; "p C "p ; X C X/ (4.97)

Si la rigidité tangente continue à l’état s est DT telle que P D DT "P, la matrice tan-
gente cohérente DC est définie par ı D DC ı". Elle relie la variation de correction de
contrainte obtenue, à la variation de pilotage en incrément de déformation totale : piloter
avec " C ı" à partir du même état s conduira à un nouvel état suivant s C ıs C ıs.
On peut obtenir DC par une méthode de perturbation (on procède à plusieurs intégrations
avec des perturbations différentes ı" sur le pilotage) mais cela est assez coûteux numéri-
quement, 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.
 Section 4.4 | Intégration de modèles de comportement 76

4.4.3 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 P : les conditions de complémentarité de Kuhn-
Tucker (4.87) sont alors remplacées par :
D f .; Y/ E n
P D (4.98)
m C

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

x C jxj
hxiC D max.0; x/ D (4.99)
2
et :
x jxj
hxi D min.0; x/ D D h xiC (4.100)
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 fN D f hf iC , alors fN D 0 correspond à f > 0 c’est-à-dire à un écoulement
viscoplastique, et fN < 0 correspond à f 6 0 c’est-à-dire à une évolution élastique, et dans
ce cas fN D f . On peut aussi se ramener à une relation de complémentarité : si on note

Ayant les incréments donnés " et t,

prédiction élastique avec (4.101)
prédiction du seuil f .; Y/
si f > 0 alors
correction viscoplastique
 D 0, fQ D f
tant que fQ > 0
P avec (4.105)
estimer "Pp et X
calculer la correction ı avec (4.106)
mise à jour de , Y, "p , X avec (4.104)
mise à jour du multiplicateur   C ı
mise à jour du seuil corrigé fQ
fin
fin

Algorithme 4.1 — Intégration du comportement viscoplastique

 D , P on doit en effet avoir  > 0, fQ D f m1=n 6 0 et fQ D 0. La principale

différence par rapport au cas de la plasticité est la dépendance avec le temps au travers
des vitesses des quantités cherchées. L’intégration de la loi de comportement ne sera plus
pilotée par le seul incrément de déformation ", mais aussi par l’incrément de temps t.
Le principe consiste encore à prédire une évolution élastique :

 D D"; Y D 0
(4.101)
X D 0; "p D 0
77 Chapitre 4 | Résolution de systèmes différentiels

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

@fQ @fQ @fQ

fQ. C ı; Y C ıY;  C ı/  fQ.; Y; / C ı C ıY C ı (4.102)
@ @Y @

pour lesquels on cherche à rester sur le seuil corrigé :

@f @f m 1=n
f .; Y/ m1=n C ı C ıY  1
ı D 0 (4.103)
@ @Y n
@G
avec ı D Dı"p , ıY D @X
ıX et les incréments :

ı"p @f
D . C ı/ "Pp
t @ (4.104)
ıX @f P
D . C ı/ X
t @Y
où les dérivées sont estimées par rapport au précédent pas de temps (et doivent être réévaluée
à chaque itéré de la correction) :

"p "pi 1 
"Pp D "Ppi
t 
(4.105)
P D X Xi 1  P
X Xi
t 
L’équation (4.102) permet alors d’obtenir l’incrément de multiplicateur viscoplastique :

f m1=n
t
C @f
@
D"Pp @f @G P
@Y @X
X
ı D @f @f @f @G @f
(4.106)
m 1=n 1 1
D C @Y @X @Y C n 
@ @ t

La valeur courante du multiplicateur doit alors aussi être mise à jour par  C ı, à chaque
itéré de la correction, qui est poursuivie jusqu’à retourner sur le seuil corrigé fQ.
Connaissant la solution à ti , 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.1. La solution à
tiC1 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,
 Section 4.5 | Cas de la dynamique non régulière 78

lorsque les structures sont soumises à des chocs ou impacts. Dans ce dernier cas, en particu-
lier 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 7 . 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 vnC ) est opposée
en direction à celle avant impact (notée vn ) :

vnC D evn (4.107)

le saut de vitesse lors de l’impact étant vnC vn . Le coefficient matériau e, appelé coefficient
de restitution d’énergie, traduit le comportement de l’impact ; 0 6 e 6 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 granu-
laires ou soumis à la fragmentation... Dans ce cas, une technique event driven est inappli-
cable à 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’emploie à 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 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 :

mxR D f (4.108)

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 8 . Pour cela, commençons par intégrer
l’équation du mouvement sur un pas de temps Œti ; tiC1  ; on obtient :
Z ti C1
[Link] iC1 xP i / D f .t/ dt (4.109)
ti

L’intégrale de la force extérieure (régulière) peut être approximée par un schéma d’inté-
gration, par exemple la méthode des trapèzes généralisés, voir section 4.1.6. Le problème
7. 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.
8. On travaillera alors à la fois sur les positions xi et les vitesses xP i aux piquets de temps ti .
79 Chapitre 4 | Résolution de systèmes différentiels

s’écrit alors :

mxP iC1 D mxP i C .1 ˛/hfi C ˛hfiC1 (4.110)

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

Z ti C1
xiC1 D xi C xP dt  xi C .1 ˛/hxP i C ˛hxP iC1 (4.111)
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 :

mvP D f
(4.112)
xP D 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 ti et tiC1 va s’écrire :
Z ti C1
[Link] iC1 xP i / D f .t/ dt C R (4.113)
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 D 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 6 x ? F > 0, qui est une
version concise de x > 0, F > 0, x
F D 0. Le pendant de ce comportement en dynamique
est la condition de Signorini-vitesse issue du Lemme de Moreau :
— si x > 0, R D 0
— si x D 0, 0 6 xP ? R > 0
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, RiC1 D 0
— si x p 6 0, 0 6 xP iC1 ? RiC1 > 0
Un prédicteur classique (saute-mouton) est x p D xi C h2 xP 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 xP iC1 par
une vitesse dite formelle :
xP iC1 xP i
vQ D (4.114)
1Ce
 Bibliographie

[1] A. A ITKEN. « On the iterative solution of a system of linear equations ». Proceedings of the
Royal Society of Edinburgh 63 1950, p. 52–60. 87

[2] E. A NDERSON, Z. BAI, C. B ISCHOF, S. B LACKFORD, J. D EMMEL, J. D ONGARRA, J. D U

C ROZ, A. G REENBAUM, S. H AMMARLING, A. M C K ENNEY et D. S ORENSEN. LAPACK
Users’ Guide. 3e édition. Philadelphia : Society for Industrial et Applied Mathematics, 1999.
ISBN : 0-89871-447-8. 8

[3] Z. BAI, J. D EMMEL, J. D ONGARRA, A. RUHE et H. van der VORST, éditeurs. Templates for
the Solution of Algebraic Eigenvalue Problems : A Practical Guide. Draft. 1999. 33

[4] R. BARRETT, M. B ERRY, T. C HAN, J. D EMMEL, J. D ONATO, J. D ONGARRA, V. E IJKHOUT,

R. P OZO, C. ROMINE et H. Van der VORST. Templates for the Solution of Linear Systems :
Building Blocks for Iterative Methods. 2e édition. Philadelphia : SIAM, 1994.
URL : [Link] 11

[5] K. BATHE et E. W ILSON. « Stability and accuracy analysis of direct integration methods ».
Earthquake Engineering & Structural Dynamics 1(3) 1972, p. 283–291. ISSN : 1096-9845. 71
DOI : 10.1002/eqe.4290010308

[6] K. BATHE et E. W ILSON. Numerical Analysis in Finite Element Analysis. 1re édition. Prentice-
Hall, 1976. ISBN : 9780136271901. 8

[7] Commandant B ENOÎT. « Note sur une méthode de résolution des équations normales prove-
nant 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
1924, p. 5–77. 18

[8] L. B LACKFORD, J. C HOI, A. C LEARY, E. D ’A ZEVEDO, J. D EMMEL, I. D HILLON, J. D ON -

GARRA , S. H AMMARLING , G. H ENRY , A. P ETITET , K. S TANLEY , D. WALKER et R. W HA -
LEY . ScaLAPACK Users’ Guide. Philadelphia : Society for Industrial et Applied Mathematics,
1997. ISBN : 0-89871-397-8. 8

[9] C. B ROYDEN. « The Convergence of a Class of Double-rank Minimization Algorithms. Ge-

neral Considerations ». IMA Journal of Applied Mathematics 6(3) 1970, p. 76–90. 51
DOI : 10.1093/imamat/6.1.76

[10] J.-P. C OMBE. « Sur le contrôle des calculs en dynamique rapide. Application aux problèmes
d’impact ». Thèse de doctorat. Laboratoire de Mécanique et Technologie : École Nationale
Supérieure de Cachan, 2000. 65

[11] M. C RISFIELD. Non-Linear Finite Element Analysis of Solids and Structures. Wiley, 1996.
ISBN : 9780471970590. 47

[12] E. C UTHILL et J. M C K EE. « Reducing the Bandwidth of Sparse Symmetric Matrices ». Pro-
ceedings of the 24th National Conference. New York, USA : ACM, 1969, p. 157–172
DOI : 10.1145/800195.805928. 14

[13] G. DAHLQUIST. « A special stability problem for linear multistep methods ». BIT. Numerical
Mathematics 3(1) 1963, p. 27–43. 56
DOI : 10.1007/BF01963532

[14] Y. DAI et Y. Y UAN. « Global convergence of the method of shortest residuals ». Numerische
Mathematik 83(4) 1999, p. 581–598. ISSN : 0029-599X. 49
DOI : 10.1007/s002119900080
 Bibliographie 82

[15] J. D ONGARRA, H. M EUER et E. S TROHMEIER. « TOP 500 supercomputer sites ». Procee-

dings of Supercomputer ’99 Conference. 17. Mannheim, 1999.
URL : [Link] 14

[16] I. D UFF. « Parallel implementation of multifrontal schemes ». Parallel Computing 3(3) 1986,
p. 193–204. ISSN : 0167-8191. 20
DOI : 10.1016/0167-8191(86)90019-0

[17] I. D UFF, A. E RISMAN et J. R EID. Direct Methods for Sparse Matrices. Oxford : Oxford
University Press, 1990. 16

[18] W. D UNCAN. « Some devices for the solution of large sets of simultaneous linear equations ».
Philosophical Magazine. 7e série 35(249) 1944, p. 660–670. 50
DOI : 10.1080/14786444408520897

[19] C. FARHAT et D. R IXEN. « Encyclopedia of Vibration ». Academic Press, 2002. Chapitre Li-
near Algebra, p. 710–720. ISBN : 9780122270857. 23

[20] R. F LETCHER et M. P OWELL. « A Rapidly Convergent Descent Method for Minimization ».

The Computer Journal 6(2) 1963, p. 163–168. 50
DOI : 10.1093/comjnl/6.2.163

[21] R. F LETCHER et C. R EEVES. « Function minimization by conjugate gradients ». The Compu-

ter Journal 7(2) 1964, p. 149–154. 49
DOI : 10.1093/comjnl/7.2.149

[22] R. F REUND et N. NACHTIGAL. « QMR : a quasi-minimal residual method for non-Hermitian

linear systems ». Numerische Mathematik 60(1) 1991, p. 315–339. ISSN : 0029-599X. 30
DOI : 10.1007/BF01385726

[23] C. G EAR. Numerical Initial Value Problems in Ordinary Differential Equation. Prentice-Hall,
1971. ISBN : 9780136266068. 71

[24] A. G EORGE. « Nested dissection of a regular finite-element mesh ». SIAM Journal on Nume-
rical Analysis 10 1973, p. 345–363. 14
DOI : 10.1137/0710032

[25] A. G EORGE et W. L IU. « The Evolution of the Minimum Degree Ordering Algorithm ». SIAM
Review 31(1) 1989, p. 1–19. ISSN : 0036-1445. 14
DOI : 10.1137/1031001

[26] M. G ÉRADIN et D. R IXEN. Théorie des Vibrations. Application à la dynamique des structures.
Physique Fondamentale et Appliquée. Masson, 1996. ISBN : 9782225851735. 33

[27] M. G ÉRADIN et D. R IXEN. Mechanical Vibrations. Theory and Application to Structural

Dynamics. 2e édition. Wiley, 1997. ISBN : 9780471975465. 33
e
[28] G. G OLUB et C. Van L OAN. Matrix Computations. 2 édition. Baltimore : The Johns Hopkins
University Press, 1989. ISBN : 9780801830105. 8

[29] L. H AGEMANA et D. YOUNG. Applied Iterative Methods. New York : Academic Press, 1981.
ISBN : 9780486434773. 24

[30] W. H AGER. « Updating the Inverse of a Matrix ». SIAM Review 31(2) 1989, p. 221–239. ISSN :
0036-1445. 50
DOI : 10.1137/1031049

[31] G. H ESSENBERG. Transzendenz von e und  : ein Beitrag zur hoheren Mathematik vom
elementaren Standpunkte aus. Berlin : B. G. Teubner, 1912. 42

[32] M. H ESTENES et E. S TIEFEL. « Methods of conjugate gradients for solving linear systems ».
Journal of research of the National Bureau of Standards 49 1952, p. 409–436. 29
83 Bibliographie

[33] H. H ILBER, T. H UGHES et R. TAYLOR. « Improved numerical dissipation for time integration
algorithms in structural dynamics ». Earthquake Engineering & Structural Dynamics 5(3)
1977, p. 283–292. ISSN : 1096-9845. 72
DOI : 10.1002/eqe.4290050306

[34] A. H OUSEHOLDER. « Unitary Triangularization of a Nonsymmetric Matrix ». Journal of the

Association of Computing Machinery 5(4) 1958, p. 339–342. ISSN : 0004-5411. 40
DOI : 10.1145/320941.320947

[35] C. H UANG et G. V ERCHERY. « An exact structural static reanalysis method ». Communica-

tions in Numerical Methods in Engineering 13(2) 1997, p. 103–112. ISSN : 1099-0887. 50
DOI : 10.1002/(SICI)1099-0887(199702)13:2<103::AID-CNM36>[Link];2-D

[36] T. H UGHES. The Finite Element Method, Linear Static and Dynamic Finite Element Analysis.
Englewood Cliffs : Prentice-Hall, 1987. ISBN : 9780133170252. 53

[37] T. H UGHES et T. B ELYTSCHKO. Nonlinear Finite Element Analysis. Zace sevices Ltd - ICE
Division, 1995. 69

[38] B. I RONS. « A frontal solution program for finite element analysis ». International Journal
for Numerical Methods in Engineering 2(1) 1970, p. 5–32. 20
DOI : 10.1002/nme.1620020104

[39] C. JACOBI. « Über ein leichtes Verfahren, die in der Theorie der Säcularstörungen vorkom-
menden Gleichungen numerish aufsulösen ». Journal für die reine und angewandle Mathe-
matik 30 1846, p. 51–94. 34

[40] D. J ORDAN et P. S MITH. Nonlinear Ordinary Differential Equations. Oxford applied mathe-
matics et computing science series, 1987. 53

[41] C. K ELLEY. Solving Nonlinear Equations with Newton’s Method. Tome 1. Fundamentals
of Algorithms. Philadelphia : Society for Industrial et Applied Mathematics, 2003. ISBN :
9780898715460. 47

[42] T. KOLDA, D. O’L EARY et L. NAZARETH. « BFGS with Update Skipping and Varying Me-
mory ». SIAM Journal on Optimization 8(4) 1998, p. 1060–1083. ISSN : 1052-6234. 51
DOI : 10.1137/S1052623496306450

[43] W. K UTTA. « Beitrag zur naherungsweisen Integration totaler Differentialgleichunge ». Zeit-

schrift für Mathematik und Physik 46 1901, p. 435–453. 61

[44] C. L ANCZOS. « 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 1950, p. 225–280. 44

[45] J. L EMAITRE et J.-L. C HABOCHE. Mechanics of Solid Materials. Cambridge University

Press, 1994. 47

[46] H. L EMOUSSU. « Approche non incrémentale pour le calcul de choc avec contact unilaté-
ral avec frottement ». Thèse de doctorat. Laboratoire de Mécanique et Technologie : École
Nationale Supérieure de Cachan, 1999. 69, 70

[47] D. L IU et J. N OCEDAL. « On the Limited Memory BFGS Method for Large Scale Optimiza-
tion ». Numerische Mathematik 45(3) 1989, p. 503–528. ISSN : 0025-5610. 51
DOI : 10.1007/BF01589116

[48] H. M EUER, E. S TROHMAIER, J. D ONGARRA et H. S IMON. « Top500 Supercomputer Sites ».

16th International Supercomputer Conference : Applications — Architectures — Trends. 2001.
URL : [Link] 14

[49] B. M OHAMMADI et J.-H. S AIAC. Pratique de la simulation numérique. Dunod / Industries et

Technologies, 2003. ISBN : 210006407X. 7
 Bibliographie 84

[50] L. NAZARETH. « A relationship between BFGS and conjugate gradient algorithms and its
implications for new algorithms ». SIAM Journal on Numerical Analysis 16 1979, p. 794–
800. 51
DOI : 10.1137/0716059

[51] N. N EWMARK. « A method for computation of structural dynamics ». Proceedings of the

American Society of Civil Engineers. Sous la direction d’EMS. Tome 85. 1959, p. 67–94. 66

[52] M. O RTIZ et E. P OPOV. « Accuracy and Stability of Integration Algorithms for Elastoplastic
Constitutive Relations ». International Journal for Numerical Methods in Engineering 21(9)
1985, p. 1561–1576. 74
DOI : 10.1002/nme.1620210902

[53] M. O RTIZ et J. S IMO. « An analysis of a new class of integration algorithms for elasto-plastic
constitutive relations ». International Journal for Numerical Methods in Engineering 23(3)
1986, p. 353–366. 63, 74
DOI : 10.1002/nme.1620230303

[54] H. O UDIN. Méthode des éléments finis. École Centrale de Nantes, France, 2008.
OAI : [Link]:cel-00341772. 11

[55] C. PAIGE et M. S AUNDERS. « Solution of sparse indefinite systems of linear equations ».

SIAM Journal on Numerical Analysis 12 1975, p. 617–629. 30
DOI : 10.1137/0712047

[56] K. PARK. « An Improved Stiffly Stable Method for Direct Integration of Nonlinear Structural
Dynamic Equations ». Journal of Applied Mechanics 42(2) 1975, p. 464–470. ISSN : 0021-
8936. 72
DOI : 10.1115/1.3423600

[57] E. P OLAK et G. R IBIÈRE. « Note sur la convergence de méthodes de directions conjuguées ».

Revue française d’informatique et de recherche opérationnelle 3(1) 1969, p. 35–43. ISSN :
0764-583X. . 49
URL : [Link]

[58] J. R AO. Advanced Theory of Vibration. John Wiley & Sons, 1989. 33

[59] Y. S AAD. Iterative Methods for Sparse Linear Systems. 2000. 11

[60] Y. S AAD et M. S CHULTZ. « GMRES : A Generalized Minimal Residual Algorithm for Sol-
ving Nonsymmetric Linear Systems ». SIAM Journal on Scientific and Statistical Computing
7(3) 1986, p. 856–869. ISSN : 0196-5204. 30
DOI : 10.1137/0907058

[61] G. S HORTLEY. « Use of Tschebyscheff-polynomial operators in the numerical solution of

boundary-value problems ». Journal of Applied Physics 24 1953, p. 392–396. 30
DOI : 10.1063/1.1721292

[62] P. S ONNEVELD. « CGS, A Fast Lanczos-Type Solver for Nonsymmetric Linear systems ».
SIAM Journal on Scientific and Statistical Computing 10(1) 1989, p. 36–52. 30
DOI : 10.1137/0910004

[63] R. S OUTHWELL. Relaxation Methods in Theoretical Physics : A Continuation of the Treatise,

Relaxation Methods in Engineering Science. Tome 1. The Oxford engineering science series.
The Clarendon Press, 1946. 26

[64] W. T INNEY et J. WALKER. « Direct solutions of sparse network equations by optimally orde-
red triangular factorization ». 55(11) 1967, p. 1801–1809. ISSN : 0018-9219. 14
DOI : 10.1109/PROC.1967.6011
85 Bibliographie

[65] H. van der VORST. « BI-CGSTAB : A Fast and Smoothly Converging Variant of BI-CG for
the Solution of Nonsymmetric Linear Systems ». SIAM Journal on Scientific and Statistical
Computing 13(2) 1992, p. 631–644. ISSN : 0196-5204. 30
DOI : 10.1137/0913035

[66] J. W ILKINSON. The Algebraic Eigenvalue Problem. New York : Oxford University Press,
1965. ISBN : 9780198534037.
[67] R. W ILLOUGHBY. « Collection of articles honouring Alston S. Householder ». ACM 18 1975,
p. 3–58. 40
A Méthodes particulières

A.1 Factorisation de cholesky et 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 U> U D I.
Pour ce faire, on peut utiliser une technique basée sur la factorisation précédente : si
on construit la matrice M D V> V, pleine et de taille réduite .n; n/ et sa factorisation
M D LL> , alors :

V> V D LL> ) L 1
V> VL >
D I ) U> D L 1
V> (A.1)

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

LU> D V> (A.2)

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

A.2 Technique d’accélération d’Aitken

Cette technique [1] 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 ,N construite par un
algorithme 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  .n 1 N
/
(A.3)
n 1 N  .n 2 N
/

Ces deux équations permettent de calculer une approximation de , N qui va servir pour
construire le nouvel itéré 0n , 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
0n D n (A.4)
n 2n 1 C n
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 D Mx, et on défini les matrices diagonales par bloc
contenant les quantités élémentaires de la façon suivante :
2 3 2 3
Ke Me
N D6
K 4 Ke 0 5 et M
7 N D6 4 Me 0 7
5 (B.1)
:: ::
: :
L’assemblage consiste à utiliser formellement une matrice d’assemblage P booléenne telle
que :
N >
K D PKP (B.2)
Comme M est issue de la même discrétisation, la même matrice d’assemblage est utilisée :
N >
M D PMP (B.3)
Considérons le système, de taille supérieure à n :
N x D N MN
KN Nx (B.4)
Il englobe le système initial car il suffit de l’écrire dans le sous-espace où xN est issu du
désassemblage d’un vecteur x de raille n : xN D P> x 1 . Avec le quotient de Rayleigh :
Nx
xN > KN
N x/ D
R.N (B.5)
Nx
xN > MN
on sait que :
N max D max R.N
N x/ (B.6)
xN ¤0

donc :
Nx
xN > KN
N max > (B.7)
x¤0
N Nx
xN > MN

Au vu du problème désassemblé, N 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 xN D P> x donc :
N > xN
xN > PKP x> Kx
N max > D > D R.x/ (B.8)
x¤0
N N > xN
xN > PMP x Mx
1. Les forces généralisées, elles, s’assembleraient de la façon suivante : f D PNf.
 Section B.2 | Système aux valeurs propres 90

donc :

N max > max R.x/ D max (B.9)

x¤0

B.2 Système aux valeurs propres

Cette fois-ci, le problème est de la forme Ax D 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 D 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 D 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 D max R.x/ (B.10)

x¤0

où :

x> Ax
R.x/ D (B.11)
x> x
donc :

x> Ax
max > (B.12)
x¤0 x> 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 x>
e xe D x x et les contributions élémentaires x Ax D
> >
N > x D x>
x> PAP e Ae xe C a où, si toutes les matrices élémentaires Ae sont positives, a > 0,
alors :
x>
e Ae xe C a x>
e Ae xe
max > > (B.13)
xe ¤0 x>
e xe x>
e xe

donc :

max > N 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 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, Ait-
ken 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.
— 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, combinatorics, and
matrix theory.
 92

— 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 mathematician,
but the Italian Encyclopedia refers to him as an Italian mathematician. They
certainly have some justification in this claim since Lagrange was born in
Turin and baptised in the name of Giuseppe Lodovico Lagrangia. The pa-
pers by Lagrange cover a variety of topics : beautiful results on the calculus
of variations, work on the calculus of probabilities... In a work on the foundations of dyna-
mics, 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 summarized 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 trans-
formed mechanics into a branch of mathematical analysis. « Lagrange, in one of the later
years of his life, imagined that he had overcome the difficulty (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 impor-
tant 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 determinant widely known.
— 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 diop-
tics and mathematical analysis. His work on lens identified mathematically
five coefficients describing the aberration of a lens, now called ’Seidel sums’. These Seidel
sums correspond to spherical aberration, coma, astigmatism, Petzval curvature and distor-
tion. He also introduced the concept of nonuniform convergence and applied probability to
astronomy.
— Born : 24 Oct 1821 in Zweibrœcken, Germany.
— Died : 13 Aug 1896 in Munich, Germany.
93 Annexe C | Repères biographiques

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 num-
bers 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.

Helmut Wielandt
Wielandt entered the University of Berlin in 1929 and there he studied
mathematics, 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 appli-
cations of mathematics, particularly applications of mathematics to biology.
In 1946, after the end of the war, Householder joined the Mathematics Di-
vision 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 impor-
tant books : The theory of matrices in numerical analysis.
 94

— 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 me-
thods 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.

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 equa-
tions in which the roots were expressed as infinite series of rational functions of the coeffi-
cients.
— Born : 30 Aug 1856 in Bremen, Germany.
— Died : 3 Jan 1927 in Göttingen, Germany.
95 Annexe C | Repères biographiques

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 calcu-
lus when we think of him. He certainly did develop the calculus by building
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 envisioned 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 centuries.
— 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 theorem 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
became 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 con-
dition in aerodynamics.
— Born : 3 Nov 1867 in Pitschen, Upper Silesia (now Byczyna), Poland.
— Died : 25 Dec 1944 in Frstenfeldbruck, Germany.
 Index

A méthode de . . . . . . . . . . . . . . . . . . . . 14
accélération Crout
d’Aitken . . . . . . . . . . . . . . . . . . . 26, 87 factorisation . . . . . . . . . . . . . . . . . . . 18
méthode de l’accélération linéaire 67, factorisation de . . . . . . . . . . . . . 19–21
72 Cuthill-MacKee . . . . . . . . . . . . . . . . . . . . 14
méthode de l’accélération moyenne
67 E
schéma de Newmark . . . . . . . . . . . 67 éléments finis . . 7, 11, 12, 19, 47, 89, 90
Aitken . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 code . . . . . . . . . . . . . . . . . . . . . . . . . . 39
technique d’accélération . . . . . 26, 87 conditionnement . . . . . . . . . . . . . . . 13
conditions aux limites . . . . . . . . . . 20
B convergence . . . . . . . . . . . . . . . . . . . 25
Broyden discrétisation . . . . . . . . . . . . . . . 14, 33
Méthode BFGS . . . . . . . . . . . . . . . . 51 problème de mémoire . . . . . . . . . . 51
super-élément . . . . . . . . . . . . . . . . . 19
C système différentiel du premier ordre
Cauchy 64
problème de . . . . . . . . . . . . . . . .53, 63 équation
Cholesky . . . . . . . . . . . . . . . . . . . . . . . . . . 91 caractéristique . . . . . . . . . . . . . . 34, 68
factorisation de . . . 18, 19, 30, 33, 87 d’équilibre . . . . . . . . . . . . . . . . . . . . 78
coût d’état . . . . . . . . . . . . . . . . . . . . . . . . . 75
CPU . . . . . . . . . . . . . . . . . . . . . . . . . . 15 d’admissibilité cinématique . . . . . 47
d’un algorithme . . 14–16, 18, 21, 42, d’admissibilité statique . . . . . . . . . 47
63, 65 d’Euler . . . . . . . . . . . . . . . . . . . . 21, 22
d’une méthode . . . . . . . . . . . . . . . . . 24 de comportement . . . . . . . . . . . . . . 78
fonction . . . . . . . . . . . . . . . . . . . . . . . 50 de la dynamique . . . . . . . . . . . . . . . 79
fonction (optimisation) . . . . . . . . . 48 différentielle . . . . . . . . . . . . . . . . . . . 79
complexité . . . . . . . . . . . . . . 15, 18, 19, 41 différentielle non-linéaire . . . . . . . 74
conditionnement . . . . . . . . . . . . 12, 31, 34 du mouvement . . . . . . . . . . . . . . . . . 79
convergence . . . . . . . . . . . . . . . . . . . . 24, 26 Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
asymptotique . . . . . . . . . . . . . . . . . . 24 équation d’ . . . . . . . . . . . . . . . . . 21, 22
cubique . . . . . . . . . . . . . . . . . . . . . . . 39 schéma d’
gradient conjugué préconditionné 30 explicite, 54, 56, 58
linéaire . . . . . . . . . . . . . . . . . . . . . . . 42 implicite, 58, 61, 74
méthode de Jacobi . . . . . . . . . . . . . 36
méthode de Lanczos . . . . . . . . . . . . 44 F
méthode des itérations inverses . . 39 factorisation
méthode des puissances . . . . . . . . . 37 de Cholesky . . . . . 18, 19, 30, 33, 87
méthode du gradient . . . . . . . . . . . .29 de Crout . . . . . . . . . . . . . . . . . . . 18–21
quadratique . . . . . . . . . . . . . . . . . . . 75 de Gauss . . . . . . . . . . . . . . . . . . . 19, 21
schéma d’Euler explicite . . . . . . . . 54 LDL . . . . . . . . . . . . . . . . . . . . . . 19, 39
taux de . . . . . . . . . . . . . . . . . . . . 25, 87 LDM> . . . . . . . . . . . . . . . . . . . . . . . . 20
valeurs propres multiples . . . . . . . 37 LU . . . . . . . . . . . . . . . . . . . . . . . . 17–19
Cramer QR . . . . . . . . . . . . . . . . . . . . . . . . 40, 41
97 Index

G récurrence(s) . . . . . . . . . . . . . . . . . . 24
Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 relation(s) . . . . . . . . . . . . . . . . . . . . . 22
élimination de . . . . . . . . . . . . . . . . . 16 système(s) 11, 25, 34, 41, 49, 65, 66
algorithme de Gauss-Seidel . . . . . 26
factorisation de . . . . . . . . . . . . . 19, 21 M
méthode de Gauss-Seidel . . . . 25, 26 matrice . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
bande . . . . . . . . . . . . . . . . . . . . . . . . . 21
H booléenne . . . . . . . . . . . . . . . . . . . . . 21
Hessenberg condensée . . . . . . . . . . . . . . . . . . . . . 17
forme de . . . . . . . . . . . . . . . . . . . . . . 42 d’amortissement . . . . . . . . . . . . 33, 66
Householder . . . . . . . . . . . . . . . . . . . . . . . 93 d’assemblage . . . . . . . . . . . . . . . . . . 89
algorithme de Householder-Businger- d’itération . . . . . . . . . . . . . . . . . 26, 68
Golub . . . . . . . . . . . . . . . . . . . . . 41 décalée . . . . . . . . . . . . . . . . . . . . . . . 39
factorisation QR . . . . . . . . . . . . . . . 41 de capacité . . . . . . . . . . . . . . . . . . . . 64
méthode de . . . . . . . . . . . . . . . . . . . . 34 de conduction . . . . . . . . . . . . . . . . . 64
transformation de . . . . . . . . . . . 40, 43 de masse . . . . . . . . . . . . . . . . . . . 65, 66
de raideur . . . . . . . . . . . . . . . . . . . . . 22
J de rigidité . . . . . . . . . . . . . . 20, 48, 66
Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
tangente, 75
algorithme de . . . . . . . . . . . . . . . . . . 35
de rotation . . . . . . . . . . . . . . . . . . . . 35
méthode de . . . . . . . . . . 25, 27, 34, 37
jacobienne . . . . . . . . . . . . . . 48–50, 75
convergence, 36
pleine . . . . . . . . . . . . . . . . . . . . . . . . . 18
K semi-définie . . . . . . . . . . . . . . . . . . . 23
Kahan sous-matrice . . . . . . . . . . . . . . . . . . . 17
théorème de . . . . . . . . . . . . . . . . . . . 26 symétrique . . . . . . . . . . . . . . . . . . . . 19
Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 symétrique définie positive . . . . . . 21
Kutta-Joukowski . . . . . . . . . . . . . . . 95 triangulaire supérieure . . . . . . . . . . 17
méthode de Runge-Kutta 62–64, 66, valeurs propres . . . . . . . . . . . . . 28, 33
74
Runge-Kutta . . . . . . . . . . . . . . . . . . .94 N
Zhukovsky-Kutta . . . . . . . . . . . . . . 95 Newmark
méthode de . . . . . . . . . . . . . 66, 69, 72
L schéma de . . . . . . . . . . . . . . . . . 67, 69
Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . 92 stabilité . . . . . . . . . . . . . . . . . . . . . . . 68
multiplicateurs de . . . . . . . . . . . 21, 27 Newton . . . . . . . . . . . . . . . . . . . . . . . . 92, 95
Lanczos . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 fluxions . . . . . . . . . . . . . . . . . . . . . . . 95
itérations de . . . . . . . . . . . . . . . . . . . 46 loi de . . . . . . . . . . . . . . . . . . . . . . . . . 78
méthode de . . . . . . . . . . . . . 34, 44, 45 méthode de . . . . . . . . . . . . . 47, 48, 75
Leibniz . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 modèle de . . . . . . . . . . . . . . . . . . . . . 80
linéaire Newton-Raphson . . . . . . . . . . . 49, 75
convergence . . . . . . . . . . . . . . . . . . . 42 quasi-newton . . . . . . . . . . . . . . . . . . 50
linéaire(s) norme . . . . . 11–13, 24, 30, 31, 37, 38, 46
équation(s) . . . . . . . . . . . . . . . . . 56, 65 infinie . . . . . . . . . . . . . . . . . . . . . . . . .37
combinaison(s) . . . . . . . . . 16, 24, 37 matricielle . . . . . . . . . . . . . . . . . . . . . 11
méthode des accélérations . . . 67, 72 normer . . . . . . . . . . . . . . . . . . . . . . . . 37, 45
 Index 98

O linéaire . . . . . . . . . . 11, 14, 34, 41, 49

orthogonalisation creux, 14
d’un sous-espace . . . . . . . . . . . 18, 87 triangulaire, 16, 18
de Gram-Schmidt . . . . . . . . . . . . . . 87 non-linéaire . . . . . . . . . . . . . . . . 25, 47
ré-orthogonalisation de grand problème
45 T
trapèze(s)
P méthode des . . . . . . . . . 57–59, 73, 79
pivot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 schéma des . . . . . . . . . . . . . . . . . . . . 58
Crout . . . . . . . . . . . . . . . . . . . . . . . . . 21
facrotisation LU . . . . . . . . . . . . . . . 17 W
partiel . . . . . . . . . . . . . . . . . . . . . . . . . 17 Wielandt . . . . . . . . . . . . . . . . . . . . . . . . . . 93
stratégie . . . . . . . . . . . . . . . . . . . . . . . 17 méthode de . . . . . . . . . . . . . . . . . . . . 39
total . . . . . . . . . . . . . . . . . . . . . . . . . . 17 Wilkinson . . . . . . . . . . . . . . . . . . . . . . . . . 93
préconditionnement . . . . . . . . . . . . . . . . 30 test de . . . . . . . . . . . . . . . . . . . . . . . . 45
gradient conjugué préconditionné 31,
51
préconditionneur . . . . . . . . . . . . . . . 30

R
Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . 94
Newton-Raphson . . . . . . . . . . . 49, 75
Rayleigh . . . . . . . . . . . . . . . . . . . . . . . . . . 93
amotissement de . . . . . . . . . . . . . . . 68
itérations inverses de . . . . . . . . . . . 39
quotient de . . . . . . . . . . . . . . . . . 34, 89
Runge . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
méthode de Runge-Kutta 62–64, 66,
74, 95
Runge-Kutta . . . . . . . . . . . . . . . . . . .94

S
Schur . . . . . . . . . . . . . . . . . . . . . . . . . . 91, 93
complément de . . . . . . . . . . . . . . . . 17
condensation de . . . . . . . . . . . . . . . .19
Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
algorithme de Gauss-Seidel . . . . . 26
méthode de Gauss-Seidel . . . . 25, 26
Southwell . . . . . . . . . . . . . . . . . . . . . . . . . 26
système
aux valeurs propres . . 33, 34, 89, 90
d’exploitation . . . . . . . . . . . . . . . . . 15
différentiel
du premier ordre, 79
non régulier, 79
dynamique . . . . . . . . . . . . . . . . . . . . 79