# Animation dynamique de corps déformables continus

Olivier Nocent

► To cite this version:

Olivier Nocent. Animation dynamique de corps déformables continus: Application à la simulation de textiles tricotés. Informatique [cs]. Université de Reims Champagne Ardenne (URCA), 2001. Français. NNT: 2001REIMS026. tel-04426043

HAL Id: tel-04426043

https://hal.univ-reims.fr/tel-04426043

Submitted on 30 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
THESE DE DOCTORAT

pour obtenir le titre de
**Docteur de l'Université de Reims Champagne-Ardenne**

présentée par
**Olivier Nocent**

Le 20 décembre 2001

**Animation dynamique de corps déformables continus**

Application à la simulation de
textiles tricotés

**Composition du jury**

Claude

SECROUN

Président

Professeur, URCA

Bruno

ARNALDI

Rapporteur

Professeur INSA, Université de Rennes

Christophe

CHAILLOU

Rapporteur

Professeur EUDIL en détachement INRIA, Lille

Jean-Marie

BACHMANN

Examinateur

Ingénieur, IFTH

Yannick

REMION

Directeur

Professeur, URCA
<figure><img src="image_1.png" /></figure>
REMERCIEMENTS

En premier lieu, je tiens à remercier les membres de mon jury de thèse :

BRUNO ARNALDI pour sa fidélité indéfectible aux travaux de recherche de l'équipe d'animation dynamique du LERI, illustrée par sa présence aux soutenances de thèse de DIDIER GILLARD, de JEAN-MICHEL NOURRIT et la mienne, ainsi qu'à la soutenance d'habilitation à diriger des recherches de YANNICK REMION.

CHRISTOPHE CHAILLOU pour son intérêt enthousiaste à l'égard des mêmes travaux, ayant donné lieu à de nombreuses rencontres avec les membres de l'équipe GRAPHIX du LIFL. Le fruit de ces discussions a permis entre autres d'intégrer notre solution au simulateur d'actes chirurgicaux SPORE développé au LIFL.

JEAN-MARIE BACHMANN pour son expertise éclairée en matière de confection de textiles tricotés, grâce à laquelle il aura su sensibiliser les profanes que nous sommes à la complexité structurelle et la diversité des tricots techniques.

CLAUDE SECROUN pour sa bienveillance et sa confiance de tous les instants. Je lui suis infiniment reconnaissant de m'avoir ouvert les portes de ce domaine merveilleux qu'est l'informatique graphique en facilitant mon inscription en DUIIC (Diplôme d'Université en Image Infographie Communication) alors que je ne disposais d'aucun bagage informatique.

YANNICK REMION pour son exigence et sa rigueur scientifique. Je tiens particulièrement à manifester mon admiration et mon respect pour YANNICK REMION qui, de par sa patience et son sens didactique, m'a aidé progresser et à tirer le meilleur de moi même durant ses trois années de thèse.

Je profite de cette occasion pour saluer mes collègues du LERI qui entretiennent au quotidien cette ambiance conviviale, propice aux échanges d'idées et à l'épanouissement personnel : « double respects » à 2R. Je remercie surtout JEROME CUTRONA, un de mes anciens étudiants devenu un collègue mais surtout un grand ami, pour sa disponibilité sans faille et son épaule solide lorsque j'étais au creux de la vague (et Dieu sait que j'en ai traversé des tempêtes !).

Une pensée à l'égard de mes proches qui ont suivi de près ou de loin cette longue aventure : ma mère FRANÇOISE, mes sœurs NATHALIE et ISABELLE, NAT&MANU, HALIMA&JEROME, GRONIKO, mat@stones.com, RICO DEL LULU, MURLOCK, BILBO, les petits noirs serrés de SEGAFREDO® (compagnons privilégiés de mes nuits blanches), les rythmes atmosphériques de THIEVERY CORPORATION.

Enfin toute ma reconnaissance s'adresse à ma compagne DOMINIQUE, dont l'amour et la patience infinie auront été mis à rude épreuve pendant ces derniers mois : « *gosto de ti, meu coração* »
<figure><img src="image_2.png" /></figure>
# TABLE DES MATIERES

NOTATIONS 5

## INTRODUCTION GENERALE

1 Préambule 9

2 Contexte 9

2.1 Contexte scientifique 9

2.2 Contexte applicatif 10

3 Objectifs 11

## CAISSE A OUTILS THEORIQUES POUR L'ANIMATION DYNAMIQUE

### CHAPITRE ①

### EQUATIONS DE LAGRANGE : DE LEUR JUSTIFICATION A LEUR MISE

EN ŒUVRE INFORMATIQUE 17

1 Introduction 17

2 Préambule historique 18

3 Définitions 20

3.1 Configuration géométrique 20

3.2 Contraintes supplémentaires 22

3.3 Cinématique 24

3.4 Dynamique 25

4 Approche intégrale 27

4.1 Lagrangien 27

4.2 Principe de moindre action 28

4.3 Equations de Lagrange 30

4.4 Expression du lagrangien 32

5 Approche différentielle 33

5.1 Mouvements virtuels 33

5.2 Vitesse et puissance virtuelles 34

5.3 Principe des puissances virtuelles 36

5.4 Principe de d'Alembert 36

5.5 Equations de Lagrange 37

5.6 Prise en compte des contraintes 38

5.6.1 Contraintes holonomes 38

5.6.2 Contraintes cinématiques 40

5.6.3 Contraintes unilatérales 41

5.6.4 Multiplicateurs de Lagrange 43

6 Petit guide à l'attention de l'informaticien 44

6.1 Quelques bonnes raisons pour privilégier le formalisme lagrangien 44

6.1.1 Généricité 44

6.1.2 Prise en compte naturelle des liaisons 45

6.2 Mise en œuvre informatique 45

6.3 Formulation équivalente des équations de Lagrange 46
6.3.1 Equations de Lagrange 46

6.3.2 Equations des contraintes 47

6.4 Architecture logicielle 49

## CHAPITRE ②

### METHODES D'INTEGRATION NUMERIQUE 51

1 Introduction 51

2 Définitions 52

2.1 Discrétisation du problème 52

2.2 Méthodes d'intégration 52

2.3 Consistance, stabilité et convergence des méthodes d'intégration 53

2.3.1 Consistance 53

2.3.2 Stabilité 54

2.3.3 Convergence 56

2.4 Classification des problèmes d'intégration numérique 56

3 Méthodes d'intégration à un pas 57

3.1 Méthodes explicites 57

3.1.1 Critères de convergence 57

3.1.2 La méthode d'Euler (1768) 58

3.1.3 La méthode de Runge (1895) 59

3.1.4 Construction des méthodes explicites à un pas (1901) 60

3.1.5 Classification des méthodes explicites à un pas 62

3.1.5.1 Méthode d'ordre 1 62

3.1.5.2 Méthodes d'ordre 2 62

3.1.5.3 Méthodes d'ordre 3 63

3.1.5.4 Méthodes d'ordre 4 63

3.1.5.5 Méthodes d'ordre 5 63

3.1.5.6 Méthodes d'ordre 6 64

3.1.6 Récapitulatif sur les méthodes explicites à un pas 64

3.1.7 Contrôle de l'erreur, sélection du pas 65

3.2 Méthodes implicites 66

3.2.1 Résolution par itérations 67

3.2.2 Linéarisation 67

3.2.3 Méthode Prédiction-Evaluation-Correction 69

4 Méthodes d'intégration à pas multiples 69

4.1 Méthodes explicites 70

4.1.1 Méthode d'Adams-Bashforth (1855) 70

4.2 Méthodes implicites 71

4.2.1 Méthode d'Adams-Moulton (1855) 71

4.2.2 La méthode BDF (Backward Differentiation Formulas) 73

5 Synthèse 73

## CHAPITRE ③

### CORPS DEFORMABLES EN ANIMATION DYNAMIQUE. ETAT DE L'ART

1 Introduction 75

2 Modèles classiques 76

2.1 Mécanique du point matériel 76

2.1.1 Modèle masses-ressorts 76

2.1.2 Modèle à base de particules 77

2.1.3 Modèle SPH 78

2.2 Mécanique des milieux continus 79

2.2.1 Modèle de Terzopoulos 80

2.2.2 Méthode des éléments finis 82

2.2.3 Modèle « dynamique modale » 84

2.2.4 Modèle « dynamique globale » 84

2.3 Synthèse 85
## VARIETES SPLINES DE DIMENSION D

### CHAPITRE ①

GEOMETRIE DES VARIETES SPLINES 91

1 Introduction 91

2 Courbe spline 91

2.1 Courbe d'Hermite 94

2.2 Courbe de Catmull-Rom 95

2.3 Courbe de Kochanek (i.e. courbe TCB) 95

2.4 Courbe de Bézier 96

3 Extension aux dimensions supérieures 96

3.1 Surface et volume spline 96

3.2 Multi-indice de dimension d 98

3.3 Variété spline de dimension d 98

3.4 Coordonnées généralisées 100

3.5 Fonction de positionnement 100

4 Codage informatique 101

### CHAPITRE ②

CINEMATIQUE ET DYNAMIQUE DES VARIETES SPLINES 105

1 Cinématique 105

1.1 Fonction des configurations 105

1.2 Fonction de vitesse 106

2 Dynamique 107

2.1 Répartition massique 107

2.2 Energie cinétique, équations de Lagrange 112

2.3 Efforts 113

2.3.1 Force ponctuelle générique 114

2.3.2 Densité d'efforts : amortissement 114

2.3.3 Energie potentielle : gravitation 115

3 Codage informatique 116

4 Exemples de variétés splines 117

4.1 Plongeoir 117

4.2 Cube élastique 118

4.3 Balançoire 119

### CHAPITRE ③

ENERGIE CONTINUE DE DEFORMATION DES VARIETES SPLINES 121

1 Introduction 121

2 Théorie de l'élasticité 122

2.1 Déformations 122

2.1.1 Déformations en description lagrangienne 122

2.1.2 Elasticité linéaire 125

3 Contraintes 126

4 Loi de comportement 127

4.1 Matériaux élastiques 127

4.2 Matériaux élastiques isotropes 128

4.3 Matériaux de St Venant-Kirchhoff 129

4.4 Matériaux hyper élastiques 129

4.5 Loi de Hooke 130
5 Application aux courbes splines matérielles 131

5.1 Aspects informatiques de l'énergie continue de déformation 137

5.2 Loi de ressort équivalente 138

5.3 Energie de déformation vs. ressorts 140

6 Conclusion 142

## CHAPITRE ④

### REDUCTION DE PARAMETRAGE : APPLICATION AUX TEXTILES

TRICOTES 143

1 Introduction 143

2 Réduction de paramétrage 144

3 Application à l'animation des textiles tricotés 147

3.1 Mise en œuvre 148

3.1.1 Volume paramétrique englobant 148

3.1.2 Plongement du liage 150

3.2 Validation qualitative 151

3.3 Construction du système différentiel « réduit » 153

3.3.1 Matrice des masses généralisées 153

3.3.2 Puissance virtuelle des efforts 154

3.4 Exemple numérique 154

4 Conclusion 156

## ■ CONCLUSION GENERALE

BIBLIOGRAPHIE 165
NOTATIONS

Les quantités scalaires sont désignées à l'aide de lettres italiques telles que le temps t.

Les quantités vectorielles ou matricielles sont désignées par des lettres en caractère gras comme un vecteur **v** ou une matrice **M**. Afin de distinguer la numérotation de quantités vectorielles ou matricielles avec l'indexation de leurs coordonnées, on choisit d'adopter la convention suivante :

□ Les numéros en exposant désignent le numéro du vecteur (resp. matrice) dans un ensemble ordonné. Exemple : $\mathbf{p}^i$ $1 \le i \le n$.

□ Les numéros en indice désignent le numéro de la coordonnée du vecteur (resp. matrice). Exemple : $\mathbf{p} \in \mathbb{R}^3 : \mathbf{p} = (p_1, p_2, p_3)$.

Pour désigner des sous-éléments d'un vecteur ou d'une matrice, on utilise la notation symbolique suivante :

$\mathbf{M}_{i0}$ correspond au vecteur ligne n° i de la matrice $\mathbf{M}$.

$\mathbf{M}_{0i}$ correspond au vecteur colonne n° i de $\mathbf{M}$.

$\mathbf{v}^*$ correspond au vecteur $\mathbf{v}$ privé de son premier élément.

$\mathbf{M}^*$, $\mathbf{M}^{**}$, $\mathbf{M}^{**}$ correspondent à la matrice $\mathbf{M}$ privée respectivement de sa première ligne, de sa première colonne, ou des 2.

Les quantités tensorielles sont quant à elles désignées par des lettres cursives comme $\mathcal{I}$, $\mathcal{F}(\mathbf{x}^0, t)$ ou $\mathcal{E}(\mathbf{x}^0, t)$.

Concernant les dérivées partielles par rapport à une quantité scalaire ou même vectorielle, on utilise les notations équivalentes suivantes :

$$ \frac{\partial f}{\partial x_i}(\mathbf{x}, t) = \partial_{x_i} f(\mathbf{x}, t) \quad \text{et} \quad = \frac{\partial f}{\partial \mathbf{x}}(\mathbf{x}, t) = \partial_{\mathbf{x}} f(\mathbf{x}, t) = \operatorname{grad}_{\mathbf{x}} f(\mathbf{x}, t) $$

Enfin par souci de concision, on désigne plusieurs espaces matriciels usuels à l'aide des notations suivantes :

$\mathbf{M}_n(\mathbb{R})$ désigne l'ensemble des matrices carrées de taille $n$ à coefficients réels.

$\mathbf{M}_n^+(\mathbb{R})$ désigne le sous-ensemble de $\mathbf{M}_n(\mathbb{R})$ des matrices positives.

$\mathbf{M}_n^>(\mathbb{R})$ désigne le sous-ensemble de $\mathbf{M}_n(\mathbb{R})$ des matrices définies positives.

$\mathbf{S}_n(\mathbb{R})$ désigne le sous-ensemble de $\mathbf{M}_n(\mathbb{R})$ des matrices symétriques.

$\mathbf{O}_n(\mathbb{R})$ désigne le sous-ensemble de $\mathbf{M}_n(\mathbb{R})$ des matrices orthogonales.
<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>

# INTRODUCTION GENERALE
<figure><img src="image_5.png" /></figure>
# 1 Préambule

Depuis plusieurs décennies, les sociétés dites modernes sont sujettes à de profondes mutations sous l'impulsion d'une informatisation massive. En effet, l'informatique a d'ores et déjà investi de nombreux secteurs d'activité demandeurs de toujours plus de performance, d'automatisation et de capacité de communication. Les conséquences de cette conquête triomphante de l'ordinateur sont palpables dans les comportements sociaux qui reflètent globalement une aspiration à une certaine forme d'instantanéité dans la communication et la transmission du savoir. Dans ce nouveau contexte, l'image devient plus que jamais le vecteur principal de l'information. A la fois riche et concise, l'image numérique (qu'elle soit issue de la réalité ou d'une réalité virtuelle) devient un objet d'intérêt central pour un large éventail de professions comme les métiers liés à l'audiovisuel, la médecine, l'industrie, la défense, etc.

En effet, l'aura de l'informatique graphique éclaire de nombreuses disciplines « classiques », leur offrant de nouvelles voies de recherche grâce à l'apport de nouvelles possibilités de manipulation, d'organisation et de visualisation de l'information. De façon quelque peu naïve, l'émergence de l'informatique graphique et l'apparition d'une cohorte de nouveaux outils semblent laisser supposer que l'on puisse à terme, représenter et simuler les moindres aspects de notre environnement, des plus extraordinaires aux plus quotidiens.

Un des défis que l'animation dynamique cherche à relever s'inscrit dans cette logique de simulation des nombreux aspects de l'environnement en choisissant de s'appuyer sur les lois de la physique pour reproduire le mouvement. La tâche est de grande ampleur si l'on considère l'infinité des variétés de mouvements ou d'évolutions des formes qui nous entourent. Qu'il s'agisse d'imiter le rebond d'une balle, l'enroulement de volutes de fumée, la trace d'empreintes de pas dans le sable pour les besoins de films d'animations ou bien de proposer des environnements virtuels pour la simulation d'actes chirurgicaux, la confection de vêtements, la modélisation informatique de matériaux déformables devient un enjeu de taille. Le travail présenté dans ce mémoire se propose de développer la construction d'un modèle physique pour représenter l'évolution de corps déformables continus. L'implantation informatique de ce modèle de corps déformables vient apporter une réponse alternative aux modèles déjà existants en animation dynamique.

Par ailleurs, l'habillement est depuis quelques années déjà un de ces secteurs industriels fortement marqués par l'intégration des techniques informatiques de pointe. En effet, la modélisation informatique de matières textiles offre de nombreux avantages tels que la faculté de concevoir, d'étudier et de visualiser des vêtements virtuels sans avoir recours à leur confection réelle sur un métier à tisser ou à tricoter. Je m'emploie alors à appliquer mon modèle de corps déformables aux textiles tricotés afin d'apporter quelques éléments de réponse aux nombreuses attentes formulées par l'industrie textile.

# 2 Contexte

## 2.1 Contexte scientifique

Le contenu de mon travail a tout d'abord été influencé par un contexte scientifique dans la mesure où il s'inscrit dans la politique scientifique à long terme de l'équipe d'animation
dynamique du LERI. Le début de ma thèse coïncida en fait avec une importante restructuration des travaux de l'équipe. Dans le cadre de leurs de thèses respectives, Didier Gillard et Jean-Michel Nourrit avaient été amenés à développer des moteurs d'animation dynamique dédiés respectivement à l'animation de corps rigides articulés et de textiles tricotés. Sous l'impulsion de Yannick Remion, l'équipe entreprit l'élaboration d'une plate-forme d'animation dynamique unifiée et générique afin de rassembler l'ensemble des savoir faire acquis au fil des travaux précédents. Mais cet effort de restructuration fut aussi l'occasion d'approfondir nos connaissances pour asseoir notre maîtrise de compétences annexes mais pourtant nécessaires à l'animation dynamique.

La première et la plus essentielle de ces disciplines étrangères à l'informatique mais indispensable est la mécanique. Pour tout informaticien que nous sommes, la compréhension des équations mathématiques du mouvement est un passage obligatoire pour la conception d'une application d'animation dynamique. Cette compréhension nous donne entre autres la faculté de choisir le formalisme physique le plus approprié pour simuler un phénomène donné. Plus encore, cette lucidité vis-à-vis des équations du mouvement nous permet d'envisager une implantation informatique performante exploitant les nombreuses propriétés structurelles des équations du mouvement retenues.

L'intégration numérique est une discipline issue de l'analyse numérique qui permet de déduire des solutions approchées des équation différentielles que sont les équations du mouvement. Ainsi l'animateur se doit de connaître les résultats principaux dans ce domaine : au moins doit-il avoir en tête une classification des méthodes les plus usuelles en maîtrisant les avantages et les inconvénients de chacune. De cette manière, il est à même, face à un problème imposé, de choisir parmi les méthodes à sa disposition celle qui offrira la solution approchée la plus acceptable.

## 2.2 Contexte applicatif

Ce travail de recherche<fnref n="1" /> est en fait partie intégrante d'une collaboration industrielle avec l'antenne de l'IFTH (Institut Français du Textile et de l'Habillement<fnref n="2" />) basée à Troyes, spécialiste dans l'étude et la confection de matériaux tricotés. Ce partenariat liant l'IFTH et le LERI (Laboratoire de Recherches et d'Etudes Informatiques) a vu le jour en 1995 sous la forme d'une première thèse réalisée par Jean-Michel Nourrit.

Les attentes de l'industrie textile sont aussi diverses que d'envergure dans la mesure où elle commence à peine à faire appel aux possibilités offertes par l'informatique graphique. Parmi les nombreuses applications de l'informatique à la modélisation de textiles, on réalise aisément les avantages qu'offrirait une plate-forme de « prototypage virtuel ». Cette application permettrait en l'occurrence d'étudier les caractéristiques visuelles et le comportement mécanique sous différentes conditions d'expérimentation d'un échantillon textile sans avoir à le confectionner. Cette démarche permettrait d'une part d'effectuer de nombreux essais en un temps moindre puisqu'il suffirait en fait de modifier quelques paramètres caractéristiques du textile et de

<footnote n="1">¹ Les travaux figurant dans ce mémoire de thèse de doctorat ont été financés par le conseil régional de Champagne-Ardenne.</footnote>

<footnote n="2">² http://www.ifth.org.</footnote>
réitérer la simulation. D'autre part, l'impact sur le coût financier serait loin d'être négligeable dans la mesure où la confection d'échantillons nécessite de monopoliser plusieurs métiers à tisser où à tricoter qui consistent en des ressources d'importance amputées à l'activité de production.

Mais l'apparition de l'informatique graphique dans la chaîne de production textile a aussi des conséquences à d'autres niveaux. Doté d'un modèle informatique fidèle à la fois à l'aspect et au comportement d'un textile, l'industrie serait à même de produire des catalogues de collections de vêtements virtuels sans avoir à les produire physiquement. Ainsi équipée de ce nouveau moyen de présentation, elle pourrait adapter de façon plus souple et plus performante l'offre à la demande, en ne confectionnant que les articles ayant les faveurs de la clientèle.

A l'extrémité de la chaîne de production, le comportement du client serait lui aussi modifié de façon drastique par l'intégration d'outils issus de l'informatique graphique. Disposant d'un modèle informatique paramétrable de vêtement, le client modifierait son comportement de consommation en achetant des vêtements sur mesure depuis son domicile. La disparition des nombreux intermédiaires dans la chaîne de distribution engendrerait une impression de proximité dans la relation qu'entretiennent producteurs et consommateurs. L'informatique devenant l'unique interface entre ces deux mondes, instaurerait une confiance née d'une personnalisation apparente du service au consommateur qui serait un moteur pour le dynamisme de l'industrie textile.

Enfin, un dernier aspect de l'apport de l'informatique graphique consiste en la proposition de modèles scientifiques de description géométrique et physique des matériaux textiles. Cette contribution de l'informatique graphique vise à combler les lacunes dans la caractérisation des matériaux tricotés surtout. Contrairement au textiles tissés dont la structure relativement simple est constituée de croisements perpendiculaires de fils, les tricots constituent un véritable défi pour la modélisation géométrique au vu de la diversité des liages que les confectionneurs ont déjà envisagés. Cette différence fondamentale entre textiles tissés et textiles tricotés est une des raisons pour laquelle il n'existe pas à ce jour de critères objectifs pour caractériser les textiles tricotés. La réalisation d'un modèle dynamique de tricots entièrement paramétrable pourrait permettre à partir de simulations de dégager des critères caractéristiques des différents modes du tricot, à l'image des tests de Kawabata qui permettent de caractériser et de classifier l'ensemble des textiles tissés.

## 3 Objectifs

Face à la diversité des contextes qui constituaient mon environnement de travail, j'orientais mes travaux dans plusieurs directions ayant pour trait commun l'enrichissement de nos réalisations en animation dynamique. Ce document reflète d'ailleurs cette diversité de centres d'intérêt par sa structure composée de deux grandes parties distinctes.

La première partie intitulée « Caisse à outils théoriques pour l'animation dynamique » regroupe les différents travaux de nature bibliographique motivés par la volonté de l'équipe de constituer un ensemble de documents de référence traitant de problèmes liés à l'animation dynamique. La première étude bibliographique de cette partie traite d'une justification des équations de Lagrange. Ces équations, pierre angulaire de la mécanique analytique, constituent le formalisme
physique pour lequel nous avons opté au sein de l'équipe d'animation dynamique. Cette étude a
pour vocation d'illustrer les nombreux concepts intervenant dans les équations de Lagrange
pour en dégager le sens physique et de les organiser pour préparer un terrain favorable à une
implantation informatique ultérieure. La deuxième étude bibliographique, rédigée dans le même
état d'esprit didactique que la première, traite des problèmes liés à la résolution numérique des
équations différentielles telles que les équations d'évolution régissant le mouvement des corps
soumis à un ensemble d'efforts. J'expose dans cette étude plusieurs méthodes usuelles
répertoriées selon des critères objectifs de comparaison afin d'éclairer le lecteur sur les
avantages et les inconvénients des méthodes proposées. Enfin, je termine cette partie consacrée
aux études bibliographiques par un état de l'art sur la modélisation des corps déformables en
animation dynamique. Ce dernier chapitre plus conventionnel dresse un inventaire des solutions
déjà proposées me permettant de situer le modèle innovant que je propose en le comparant aux
modèles existants.

Ce travail bibliographique, constituant une partie importante de ce mémoire, découle d'une
volonté forte de rassembler, d'approfondir et de figer l'expérience commune du groupe dans la
maîtrise de compétences étrangères à l'informatique afin de constituer un ensemble de
documents de référence à l'intention des membres de l'animation dynamique et j'ose l'espérer à
l'intention de tout informaticien désireux d'entreprendre la mise en œuvre d'une plate-forme
d'animation dynamique.

La conception d'une plate-forme d'animation dynamique unifiée m'a amené à proposer une
refonte du modèle informatique de textiles tricotés présenté par Jean-Michel Nourrit afin de
répondre aux impératifs du nouveau moteur générique. Cette circonstance opportune fut
l'occasion de proposer à la fois une unification et une abstraction de la démarche de conception
du modèle mécanique précédent.

Un échantillon de textile tricoté est représenté par les fils qui le constituent. Ainsi, le modèle
que je qualifie d'initial pour le distinguer de l'extension que j'en propose utilise des courbes
splines matérielles (c'est-à-dire munies d'une répartition massique et d'une élasticité interne)
pour représenter les fils. La structure du tricot déterminée par l'agencement des mailles est
traduite par des contraintes géométriques de contact entre courbes splines. En résumé, le modèle
mécanique sous-jacent est complètement déterminé par la définition d'une courbe spline
matérielle.

La seconde partie du document intitulée « Variétés splines de dimension *d* » retrace étapes par
étapes la construction d'un modèle mécanique pour les corps déformables qui consiste en une
généralisation du modèle de courbe spline matérielle. En premier lieu, je propose une
expression d'une courbe spline indépendante de la nature de ses fonctions de mélange. Cette
première unification vient enrichir le nombre de type de courbes prises en compte par le modèle.
En effet, pour des raisons pratiques, Jean-Michel Nourrit avait délibérément limité son travail au
cas de courbes de Catmull-Rom. Au vu de la forte similarité d'expression mathématique des
courbes, surfaces et volumes splines, je propose une deuxième abstraction visant à unifier
l'expression des objets splines vis-à-vis de leur dimension intrinsèque. Le modèle de variétés
splines est le fruit de ces deux étapes d'abstraction qui rassemblent dans un même contexte
théorique une grande variété de formes.
Ce modèle de nature géométrique est alors envisagé dans un contexte cinématique dans lequel on est amené à définir les degrés de liberté d'une variété qui conditionnent son évolution au fil du temps. L'ajout d'une répartition massique compatible avec la définition géométrique permet d'envisager une formulation particulière des équations du mouvement. En effet, la nature très spécifique d'une variété spline, dont la position des points est issue d'une combinaison linéaire de points de contrôle pondérée par des fonctions d'influence, permet d'opérer des simplifications formelles drastiques pour déboucher sur une expression équivalente optimale des équations de Lagrange.

L'implantation informatique des variétés splines tire parti des avantages de la conception objet tels que le polymorphisme et la généricité pour proposer un modèle informatique fidèle à l'expression mathématique unifiée mais aussi performant en exploitant les propriétés particulières de fonctions de mélange spécifiques.

Ce travail s'inscrit dans une démarche globale et générique visant à proposer un modèle physique et informatique polyvalent pour la représentation de corps déformables à priori quelconques.

Les textiles tricotés constituent alors une des applications envisageables de ce modèle. Le dernier chapitre de la seconde partie développe la façon dont j'utilise le modèle de variété spline pour envisager une modification du modèle initial de tricots afin de concilier rigueur de représentation et performance.
<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /></figure>

CAISSE A OUTILS
THEORIQUES
POUR L'ANIMATION
DYNAMIQUE

L'animation dynamique fait appel à des disciplines scientifiques autres que l'informatique telles que les mathématiques, la mécanique, l'analyse numérique, etc. Dans ces conditions, l'informaticien ayant enfilé la combinaison de « l'animateur » doit au préalable se doter d'une caisse à outils théoriques s'il veut aboutir avec succès à l'élaboration d'une application de production d'animations régies par les lois de la physique.

Cette nécessité pour l'informaticien de maîtriser des savoirs hors de sa compétence initiale s'imposa très tôt aux membres de l'équipe d'animation dynamique du LERI. Ce n'est pas faire injure à leurs auteurs que de reconnaître que les thèses successives de Didier Gillard et Jean-Michel Nourrit auront permis de débroussailler et d'identifier les problèmes-clefs suscités par la réalisation d'une plate-forme d'animation dynamique. Mais pour des raisons

pragmatiques privilégiant des résultats intermédiaires plutôt qu'une longue et fastidieuse étude théorique, nombre de problèmes restaient en suspens.

Mon arrivée en thèse coïncida avec une profonde restructuration des travaux de l'équipe d'animation qui aboutit à la mise en œuvre d'une plate-forme générique et robuste pour l'animation dynamique. Cette restructuration s'accompagnait d'une forte volonté de rassembler et d'organiser l'ensemble des savoirs faire acquis au fil des années et d'approfondir les thèmes jusque là non encore explorés. Le contenu de cette première partie est le fruit de ma contribution à cette volonté de cristalliser la maîtrise de compétences nécessaires à l'animation dynamique à l'aide d'un document de « référence » pour notre groupe d'animation et plus vastement, je l'espère, pour notre communauté.
<figure><img src="image_8.png" /></figure>
# EQUATIONS DE LAGRANGE : DE LEUR JUSTIFICATION A LEUR MISE EN ŒUVRE INFORMATIQUE

## 1 Introduction

Une étude bibliographique sur les équations de Lagrange peut, au premier abord, paraître hors sujet au sein d'un mémoire de thèse de doctorat en informatique graphique. Mais en observant l'évolution des travaux de recherche en animation par ordinateur, on réalise que l'usage des équations de la dynamique est de plus en plus fréquent et ne se cantonne plus désormais à la simple utilisation de la loi fondamentale de la dynamique régissant le mouvement d'un point matériel. En effet, les informaticiens n'hésitent plus à recourir à des formalismes physiques de plus en plus complets afin de simuler des phénomènes de plus en plus complexes [O'Brien 1999] [Yngve 2000]. Ainsi l'informaticien, désireux de proposer des modèles informatiques de plus en plus riches pour l'animation, doit maîtriser ces nombreux concepts : au moins doit-il appréhender le sens physique des formalismes qu'il a à sa disposition afin de faire les choix les plus pertinents pour la modélisation informatique d'un phénomène donné. Sous l'impulsion de Yannick Remion, l'équipe d'animation dynamique du LERI a, dès ses premiers travaux [Gillard 1997] [Nourrit 1999], opté pour le formalisme lagrangien : un formalisme puissant permettant de déduire les équations qui conditionnent le mouvement d'une grande variété de systèmes matériels tels que les systèmes de particules, les solides rigides et les solides déformables. Au début de ma thèse j'ai dû, comme les autres doctorants de l'équipe d'animation dynamique avant moi, faire face aux difficultés que représentaient alors les équations de Lagrange. Pour satisfaire ma curiosité intellectuelle et l'envie de devenir plus qu'un « utilisateur » du formalisme lagrangien, j'entrepris à cette époque une étude bibliographique sur les équations de Lagrange qui n'avait initialement pas vocation à être publiée. Assez rapidement, je me trouvais dans un état de grande confusion devant la diversité des approches adoptées par les nombreux ouvrages consacrés à ce sujet [Fox 1967] [Germain 1986] [Arnold 1989] [Calkin 1996] [Bérest 1997]. En plus d'une terminologie rarement unique, les auteurs proposaient des justifications des équations de Lagrange de natures différentes. Il aura fallu attendre le moment où je m'hasardais à la lecture d'un ouvrage d'épistémologie [Barroso Filho 1994] pour commencer à cerner les grandes lignes directrices de la démarche de Joseph Louis Lagrange, me permettant
d'organiser les nombreuses pistes déjà explorées en un raisonnement cohérent. Sur le conseil
avisé de Yannick Remion, je commençais la rédaction d'une synthèse des idées principales
glanées à droite et à gauche afin de constituer un document de référence pour l'équipe
d'animation dynamique du LERI.

Ce document retrace donc la chronologie des étapes majeures du développement des équations de Lagrange, étapes trop souvent éparpillées dans la littérature scientifique contemporaine, en mentionnant les contributions respectives d'illustres savants. J'insiste sur l'existence de deux justifications des équations de Lagrange issues de deux approches intellectuelles distinctes proposées à des époques différentes. Mais il ne s'agit pas de présenter ici des démonstrations mathématiques rigoureuses de ces équations, en effet mon ambition est tout autre : avec ce document, j'espère guider le lecteur au travers des concepts touffus de ce formalisme en lui évitant les travers que j'ai pu rencontrer. Ainsi, les éléments constitutifs des équations de Lagrange sont introduits progressivement, justifiés systématiquement et détaillés dans une optique d'implantation informatique. Je profite aussi de cette étude bibliographique pour fixer la terminologie ainsi que les notations qui seront intensivement utilisées par la suite.

Dans un souci de cohérence, je commence par présenter une chronologie des évènements majeurs qui ont permis d'aboutir à la construction de la mécanique analytique. Puis, après avoir introduit les définitions des éléments fondateurs, je présente successivement les deux justifications des équations de Lagrange en insistant sur le principe dont elles découlent et leur champ d'application. Je tiens à mentionner que, pour les besoins de son habilitation à diriger des recherches [Remion 2000], Yannick Remion a prolongé ce travail en exposant un développement plus rigoureux de la seconde preuve des équations de Lagrange.

## 2 Préambule historique

Dans la chronologie des étapes majeures du développement des équations de Lagrange (Figure 1), je choisis comme point de départ un mémoire d'optique géométrique rédigé par Fermat en 1657 traitant de la réflexion et de la réfraction de la lumière dans des milieux non homogènes. Dans ce document, Fermat faisait la remarque suivante : « ... la nature agit toujours de manière simple et par moindre dépense ... ». Cette idée de moindre dépense s'imposa à l'esprit des scientifiques du siècle des lumières comme un principe fondateur. Mauperthuis fut le premier, en 1744, à utiliser ce principe dans le cadre de la mécanique sans réussir à justifier son argumentation. Euler, convaincu du bien-fondé du raisonnement de Mauperthuis, rédigea un mémoire en 1753 dans lequel il présentait les premiers essais d'une formalisation mathématique du principe de moindre action.

Parallèlement aux efforts accomplis en mécanique, les mathématiques étaient aussi en plein
essor grâce notamment à l'introduction par Leibniz du calcul différentiel en 1684. Ce calcul
offrait une méthode rigoureuse et systématique pour la résolution de problèmes de recherche
d'extrema. Le XVIIIème siècle fut marqué par une vive démarche de formalisation
mathématique dans le but de faire converger des disciplines à l'époque distinctes telles que la
géométrie et l'algèbre. Lagrange fut un acteur majeur de cette démarche de formalisation et
présenta en 1759 son calcul des variations (nom donné par Euler avec lequel il entretenait une
riche correspondance). Lagrange offrait alors à la communauté scientifique un outil
mathématique permettant une formalisation algébrique de la mécanique. Choisissant le principe de moindre action comme principe fondateur de la mécanique, Lagrange déboucha sur des équations qui conditionnent le mouvement de systèmes matériels soumis à des efforts conservatifs. Cet événement majeur marqua un tournant dans l'étude du mouvement : la mécanique de Lagrange, algébriquement formalisée, s'éloignait des considérations géométriques. Pour reprendre une citation de Wilton Barroso Filho, « *Lagrange pense la mécanique dans le formalisme* » [Barroso Filho 1994]. La première preuve, issue d'une approche que je qualifie d'intégrale, s'appuie sur ces résultats.

<figure><img src="image_9.png" /><figcaption>Figure 1 Frise chronologique rassemblant les étapes majeures du développement du formalisme lagrangien.</figcaption></figure>

Il est important de mentionner que Lagrange n'était pas le seul scientifique cherchant à formaliser la mécanique. Dans son « Traité de la dynamique » rédigé en 1743, d'Alembert tentait de trouver un principe fondateur et unificateur pour la description du mouvement. Influencé par la loi fondamentale de la dynamique énoncée par Newton en 1687 dans ses « Principes Mathématiques », il dirigea ses recherches vers le **principe des puissances virtuelles** issu de la statique qui l'amena à énoncer un principe s'appliquant à la dynamique qui porte désormais son nom. Lagrange, prenant conscience de la qualité du travail de d'Alembert, utilisa son principe et le calcul des variations pour la rédaction d'un mémoire sur la libration de la lune : mémoire qui lui fit décrocher le prix de l'Académie des Sciences de Paris en 1764. Mais, c'est surtout en 1788, lors de la parution de la « Mécanique Analytique » que Lagrange rendit hommage au travail de d'Alembert. Dans cet ouvrage, synthèse d'une vie consacrée à la
formalisation de la mécanique, Lagrange choisit d'adopter le principe de d'Alembert comme principe fondateur. Il en déduisit ses équations, élargissant leur domaine d'application aux systèmes matériels soumis à des efforts dissipatifs tels que les forces de frottement. La deuxième preuve détaille cette démarche que je choisis de qualifier d'approche différentielle.

Ces travaux de Lagrange ont contribué à réduire et ordonner les principes connus de la mécanique : le formalisme lagrangien est le résultat de la démarche d'un mathématicien qui parcourt la physique, sans rapport direct avec l'observation. Les équations de Lagrange s'imposèrent alors naturellement par leur indépendance vis à vis de la représentation géométrique du problème mais aussi par leur puissance unificatrice qui avait réconcilié les partisans du principe de moindre action et ceux du principe des puissances virtuelles.

# 3 Définitions

Je rassemble dans ce paragraphe l'ensemble des définitions importantes communes aux deux preuves des équations de Lagrange. Dans un souci d'homogénéité, je choisis d'adopter les notations couramment utilisées au sein de l'équipe d'animation dynamique du LERI.

## 3.1 Configuration géométrique

Un **système matériel** est constitué d'un ensemble de corps matériels discrets ou compacts. La **configuration** d'un système matériel synthétise à la fois la position et la forme instantanée de celui-ci. Ainsi, le **mouvement** d'un système matériel correspond à la succession de ses configurations au fil du temps.

De toute évidence, la notion de mouvement est fortement dépendante des conditions d'observation, à savoir le déplacement de l'observateur et l'origine du temps. On nomme **référentiel d'observation** la donnée d'un repère éventuellement mobile lié à l'observateur et d'une origine sur la droite décrite par le temps. Dorénavant, je limiterai mon étude au cas des référentiels galiléens dans lesquels toute particule libre, ne subissant aucune force extérieure, possède un mouvement perçu comme rectiligne uniforme.

La **cinématique** d'un système matériel est la caractérisation mathématique du mouvement de ce système dans un référentiel d'observation donné. Historiquement, on distingue deux approches distinctes pour exprimer la cinématique d'un système matériel (Figure 2). La **description eulerienne** s'appuie sur la connaissance du champ des vitesses instantanées exprimant, en tout point **p** de l'espace, la vitesse de la particule occupant la position **p**. Il est important de noter que l'on ignore a priori de quelle particule il s'agit. Cette approche est à privilégier dans l'étude de l'écoulement des fluides dont les conditions d'observation ne permettent pas de suivre chaque particule mais plutôt de mesurer les vitesses instantanées en tout point de l'espace. La **description lagrangienne** s'appuie sur la connaissance des trajectoires particulaires. Cette approche suppose la possibilité d'identifier chacune des particules à partir d'un état de référence. L'étude du mouvement d'un système matériel cohésif (i.e. un système matériel dont deux particules voisines dans un état de référence demeurent proches au fil du temps) est naturellement envisageable dans une description lagrangienne.

Ainsi pour une description lagrangienne de la cinématique d'un système matériel, chaque particule constitutive d'un corps matériel est identifiée par un élément $\omega$ d'un ensemble $\Omega$.
Pour un corps compact, $\Omega$ est une partie compacte de $\mathbb{R}^d$ et une partie finie de $\mathbb{N}^d$ pour un corps discret. L'entier $d$ correspond à la dimension intrinsèque du corps matériel. On débouche naturellement sur la définition de **fonction particulière de position**.

$$ \chi : \begin{cases} \mathbb{R} \times \Omega \to \mathbb{R}^3 \\ (t, \omega) \mapsto \chi(t, \omega) \end{cases} \qquad (1.1) $$

$\chi(t, \omega)$ correspond à la position de la particule $\omega$ à l'instant $t$ lors de son mouvement.

<figure><img src="image_10.png" /><figcaption>Figure 2 Descriptions lagrangienne (a) et eulerienne (b) du mouvement.</figcaption></figure>

Le mouvement d'un corps matériel discret est entièrement déterminé par le nombre fini de fonctions de position de ses particules. Par contre, le mouvement d'un corps matériel compact constitué d'une infinité de particules ne peut être appréhendé de la même manière. Une démarche classique en mécanique consiste à réduire l'ensemble des configurations envisagées pour le système considéré en s'appuyant sur l'hypothèse d'une certaine cohésion interne.

Cette hypothèse m'amène à considérer que la configuration d'un corps matériel même compact dépend d'un nombre fini de paramètres. On appelle **coordonnées généralisées** (i.e. degrés de liberté) d'un système matériel les $n$ paramètres indépendants et indispensables permettant de déterminer exactement la **configuration** du système.

$$ \mathbf{u} = (u_1, \dots, u_n)^T \qquad (1.2) $$

Le vecteur $\mathbf{u}$ des coordonnées généralisées permet de déterminer la position de chaque particule du système dans cette configuration. Il n'est pas rare que la configuration d'un système matériel dépende explicitement du temps. Afin de prendre en compte cette dépendance explicite du temps de façon naturelle, j'introduis la notion de **configuration étendue** correspondant à un vecteur de dimension $n+1$ défini comme suit :

$$ \eta = (t \mathbf{u}^T)^T \qquad (1.3) $$
Pour la commodité des notations, je choisirai le système d'indexation suivant :

$$ 
\begin{align} 
\eta_0 &= t \tag{1.4} \\ 
\eta_i &= u_i, \quad \forall i:1 \le i \le n \nonumber 
\end{align} 
$$

L'espace décrit par les (n+1)-uplets représentant les configurations réellement accessibles par le système est une partie de $\mathbb{R}^{n+1}$. Cet espace que l'on note $\mathbb{U}$ est appelé **espace des configurations étendues**.

La configuration d'un système matériel doit permettre de déterminer la position de chaque particule constitutive du système. Ainsi La **fonction de positionnement** $\Phi$ exprime, à partir d'un descripteur $\eta$ de la configuration du système matériel et d'un identifiant de particule $\omega$, la position dans l'espace $\mathbb{R}^3$ de cette particule :

$$ \Phi: \begin{cases} \mathbb{U} \times \Omega \to \mathbb{R}^3 \\ (\eta, \omega) \mapsto \Phi(\eta, \omega) \end{cases} \quad (1.5) $$

## 3.2 Contraintes supplémentaires

Comme je l'ai détaillé dans le paragraphe 2.1, une configuration géométrique du système matériel est caractérisée par un vecteur $\eta$ de l'espace des configurations étendues $\mathbb{U}$.

Si à tout élément $\eta$ de $\mathbb{R}^{n+1}$ correspond une configuration géométrique accessible du système (c'est-à-dire $\mathbb{U} = \mathbb{R}^{n+1}$), le paramétrage est qualifié de **complet**. Dans le cas contraire (c'est à dire $\mathbb{U} \subsetneq \mathbb{R}^{n+1}$), le paramétrage est qualifié d'**incomplet** (Figure 3). L'ajout de contraintes supplémentaires constitue alors une démarche analytique permettant de caractériser l'espace des configurations étendues $\mathbb{U}$ comme un sous-ensemble strict de $\mathbb{R}^{n+1}$.

<figure><img src="image_11.png" /></figure>

### Exemple de paramétrage complet

**Pendule rigide simple** : système constitué d'une tige dont une extrémité est fixée en $O$, contrainte à se déplacer dans le plan $yOz$ autour de l'axe $Ox$.

1 coordonnée généralisée $\eta_1$ = angle de la tige avec l'axe $Oy$.

Chaque réel $\eta_1 \in \mathbb{R}$ correspond à une configuration accessible par le système.

<figure><img src="image_12.png" /></figure>

### Exemple de paramétrage incomplet

**Pendule rigide double** : système constitué de deux tiges dont une des extrémités est fixée en $O$, contraintes à se déplacer dans le plan $yOz$ autour de l'axe $Ox$.

2 coordonnées généralisées $\eta_1, \eta_2$ = angles respectifs de chaque tige avec l'axe $Oy$.

Tous les couples $(\eta_1, \eta_2) \in \mathbb{R}^2$ ne correspondent à une configuration accessible par le système puisque les tiges ne peuvent s'interpréner.

<caption>Figure 3 Exemples de paramétrages complet et incomplet.</caption>
Ces contraintes supplémentaires interviennent dans le formalisme lagrangien par le biais de conditions analytiques imposées aux coordonnées généralisées $\eta$. Les contraintes peuvent être répertoriées de plusieurs manières.

Les contraintes **holonomes** s'expriment sous la forme d'une équation implicite des coordonnées généralisées :

$$ \varphi(\eta) = 0 \qquad (1.6) $$

En présence d'une contrainte holonome, l'espace des configurations étendues $\mathcal{U}$ est restreint à une hyper surface de $\mathbb{R}^{n+1}$ (Figure 4).

<figure><img src="image_13.png" /><figcaption>Figure 4 Interprétation géométrique d'une contrainte holonome dans l'espace des configurations, projeté de l'espace des configurations étendues pour $t=cte$.</figcaption></figure>

Les **contraintes unilatérales** s'expriment sous la forme d'une inéquation implicite des coordonnées généralisées :

$$ \varphi(\eta) \ge 0 \qquad (1.7) $$

En présence d'une contrainte unilatérale, l'espace des configurations étendues $\mathcal{U}$ est restreint à un demi-espace fermé de $\mathbb{R}^{n+1}$ dont la frontière est l'hyper surface $\{\eta \in \mathbb{R}^{n+1} : \varphi(\eta) = 0\}$ (Figure 5).

<figure><img src="image_14.png" /><figcaption>Figure 5 Interprétation géométrique d'une contrainte unilatérale dans l'espace des configurations, projeté de l'espace des configurations étendues pour $t=cte$.</figcaption></figure>
Les contraintes scléronomes ne dépendent pas explicitement du temps :

$$ \frac{\partial \varphi}{\partial t}(\eta) = 0 \qquad (1.8) $$

Les contraintes rhéonomes dépendent explicitement du temps :

$$ \frac{\partial \varphi}{\partial t}(\eta) \neq 0 \qquad (1.9) $$

En imposant plusieurs contraintes simultanément, l'espace des configurations étendues $\mathcal{U}$ est déduit de l'intersection des sous-ensembles de $\mathbb{R}^{n+1}$ correspondant respectivement à chaque contrainte.

Enfin, Les contraintes **cinématiques** sont des équations qui dépendent à la fois des coordonnées généralisées $\eta$ et linéairement des taux de variation $\hat{\eta}$ des coordonnées, communément appelés vitesses généralisées, qui évoluent dans un espace de même dimension que $\mathcal{U}$ noté $\hat{\mathcal{U}}$ :

$$ \alpha(\eta) \cdot \hat{\eta} = 0 \qquad (1.10) $$

On déduit de l'équation (1.10) qu'en présence d'une contrainte cinématique, les vitesses généralisées $\hat{\eta}$ d'une configuration $\eta$ sont restreintes au plan tangent d'une hyper surface de vecteur normal $\alpha(\eta)$.

## 3.3 Cinématique

Puisque je m'intéresse à l'étude du mouvement, perçu comme une succession de configurations instantanées représentée par une suite de descripteurs $(\eta')_{t \in \mathbb{R}}$ indexée par le temps, je définis la **fonction des configurations étendues** comme une courbe paramétrée de l'espace des configurations étendues $\mathcal{U}$. Cette fonction notée $\varepsilon$ synthétise l'évolution temporelle des configurations instantanées.

$$ \varepsilon: \begin{cases} \mathbb{R} \to \mathcal{U} \\ t \mapsto \varepsilon(t) = \eta' \end{cases} \qquad (1.11) $$

Muni de cette fonction des configurations, je peux établir une relation simple entre la fonction particulière de position définie en (1.1) et la fonction de positionnement (1.5). En effet, on a :

$$ \chi(t, \omega) = \Phi(\varepsilon(t), \omega) \qquad (1.12) $$

S'appuyant sur la relation (1.12), je suis en mesure de donner une définition de la **vitesse particulière** à partir de la jacobienne de la fonction de positionnement. Pour cela, on écrit :

$$ v(t, \omega) = \dot{\chi}(t, \omega) = \frac{\partial \Phi}{\partial \eta}(\varepsilon(t), \omega) \cdot \dot{\varepsilon}(t) \qquad (1.13) $$
A l'image de la fonction de positionnement, j'introduis une **fonction de vitesse** qui, à partir d'une configuration $\eta$, du taux de variation $\hat{\eta}$ de cette configuration et d'un paramètre $\omega$ caractérisant une particule, donne la vitesse dans l'espace $\mathbb{R}^3$ de cette particule :

$$ \hat{\Phi} : \begin{cases} \mathbb{U} \times \hat{\mathbb{U}} \times \Omega \rightarrow \mathbb{R}^3 \\ (\eta, \hat{\eta}, \omega) \mapsto \hat{\Phi}(\eta, \hat{\eta}, \omega) = \frac{\partial \Phi}{\partial \eta}(\eta, \omega) \hat{\eta} \end{cases} \qquad (1.14) $$

Par conséquent, on obtient la relation suivante entre la vitesse particulaire et la fonction de vitesse :

$$ v(t, \omega) = \hat{\Phi}(\varepsilon(t), \dot{\varepsilon}(t), \omega) \qquad (1.15) $$

Cette démarche originale, qui permet de réaliser une séparation du modèle cinématique, est une des contributions importantes de Lagrange à l'essor d'une mécanique mathématiquement formalisée. Dans cette description éloignée des particularités géométriques, on trouve d'un coté les éléments structurels liés à la nature du système matériel tels que la fonction de positionnement et la fonction de vitesse. De l'autre coté, on considère le mouvement défini à partir des fonctions de configuration. Cette séparation me permet d'exprimer la cinématique du système au sein d'un mouvement quelconque $\varepsilon^*(t)$. Il suffit pour cela de substituer au **mouvement réel** $\varepsilon(t)$ l'expression du mouvement $\varepsilon^*(t)$, qualifié de **mouvement virtuel**.

Les mouvements virtuels, notamment les mouvements voisins du mouvement réel, jouent un rôle très important en mécanique lagrangienne : aussi bien dans l'approche intégrale où le principe de moindre action stipule que le mouvement réel réalise un minimum d'une quantité caractéristique du système parmi un ensemble donné de mouvements virtuels, que dans l'approche différentielle où le principe de d'Alembert repose sur l'expression des puissances virtuelles des efforts pour des mouvements virtuels proches du mouvement réel.

## 3.4 Dynamique

La seconde étape de la modélisation mécanique s'articule autour de la dynamique du système qui correspond à la caractérisation mathématique des causes de la modification du mouvement. Ces causes sont à la fois représentées par la masse inertielle du système et les efforts qui lui sont appliqués.

La masse inertielle exprime la façon dont est répartie la matière au sein du système matériel. On la définit par le biais d'une fonction de densité massique $\rho$ dépendant explicitement du temps $t$ et de l'identifiant particulière $\omega \in \Omega$.

$$ \rho : \begin{cases} \Omega \times \mathbb{R} \rightarrow \mathbb{R} \\ (\omega, t) \mapsto \rho(\omega, t) \end{cases} \qquad (1.16) $$
Le terme $\rho(\omega, t)$ équivaut, à l'instant $t$ et pour une particule identifiée par $\omega$, à :

* la masse de la particule dans le cas d'un système matériel discret ($\Omega \subset \mathbb{N}^d$).

* la densité de masse au voisinage de cette particule pour un système matériel compact ($\Omega \subset \mathbb{R}^d$). Ainsi la masse à l'instant $t$ d'un voisinage matériel infinitésimal $d\omega$ de la particule $\omega$ est égale à $\rho(\omega, t)d\omega$.

D'autre part les efforts, autres causes de modification du mouvement, peuvent être classés de différentes façons :

* On parle d'effort **intérieur** au système matériel lorsque cet effort est issu de l'interaction de deux corps appartenant au système.

* L'effort est **extérieur** au système matériel si le corps responsable de cet effort n'appartient pas au système.

* Les efforts donnés sont des efforts que l'on peut écrire explicitement en fonction des coordonnées et vitesses généralisées (i.e. configuration cinématique) du système matériel.

* Les efforts de liaison sont des efforts issus de contraintes liant les coordonnées généralisées et leurs taux de variation, dont on ne dispose pas d'une expression explicite.

L'ensemble des efforts donnés peut être décomposé en :

* un ensemble de forces ponctuelles $F^i(\varepsilon(t), \dot{\varepsilon}(t))$ appliquées à des particules $\omega^i$ ($1 \le i \le m$).

* des densités $F^i(\varepsilon(t), \dot{\varepsilon}(t), \omega)$ d'efforts appliqués à des sous-ensembles particulaires $\Omega^{F^i}$ du système matériel ($1 \le i \le m'$).

Un effort de densité $F^i(\varepsilon(t), \omega)$ (resp. un effort ponctuel $F^i(\varepsilon(t))$ appliqué en $\omega^i$) est **conservatif** lorsque son travail estimé entre deux points ne dépend pas du chemin parcouru pour joindre ces deux points. On démontre que cette hypothèse est équivalente à l'existence d'un potentiel réel $U: \mathbb{R}^3 \to \mathbb{R}$ dont l'effort $F(\varepsilon(t), \omega)$ dérive.

$$ \exists U: \mathbb{R}^3 \to \mathbb{R}: \begin{aligned} & F^i(\varepsilon(t), \omega) = -\operatorname{grad} U(\Phi(\varepsilon(t), \omega)) \\ & \quad (\text{resp. } F^i(\varepsilon(t)) = -\operatorname{grad} U(\Phi(\varepsilon(t), \omega^i))) \end{aligned} \quad (1.17) $$

Parmi ces efforts conservatifs, on distingue le cas particulier des efforts de gravitation. Ces efforts induits par la masse dérivent de l'énergie potentielle de gravitation qui s'exprime de la façon suivante :

$$ U^{\text{gravitation}}(\Phi(\varepsilon(t), \omega)) = \rho(\omega, t)\Phi(\varepsilon(t), \omega) \cdot \vec{g} \quad (1.18) $$
Le vecteur $\vec{g}$ est orienté selon la verticale gravitationnelle, sa norme $||\vec{g}||$ est la constante de gravitation locale. Cela suppose que le référentiel terrestre local choisi pour l'étude du mouvement est galiléen : hypothèse vérifiée pour des durées d'observations courtes.

La masse intervient aussi dans la définition de l'**énergie cinétique**, autre quantité caractéristique de la dynamique d'un système matériel dont l'expression est donnée par l'équation (1.19) :

$$K(t) = \frac{1}{2} \int_{\Omega} \rho(\omega, t) \dot{\chi}^2(t, \omega) d\omega = \frac{1}{2} \int_{\Omega} \rho(\omega, t) \hat{\Phi}(\varepsilon(t), \dot{\varepsilon}(t), \omega)^2 d\omega \quad (1.19)$$

Exploitant l'expression particulière de la fonction de vitesse (1.14), on modifie la formulation de l'énergie cinétique pour aboutir à une expression équivalente.

$$ 
\begin{gathered}
K(t) = \frac{1}{2} \dot{\varepsilon}(t)^T \mathbf{M}(\varepsilon(t)) \dot{\varepsilon}(t) \\
\text{avec } \mathbf{M}(\varepsilon(t)) = \int_{\Omega} \rho(\omega, t) \frac{\partial \Phi}{\partial \eta}(\varepsilon(t), \omega)^T \frac{\partial \Phi}{\partial \eta}(\varepsilon(t), \omega) d\omega
\end{gathered}
\quad (1.20)
$$

Cette nouvelle expression de l'énergie cinétique rejoint la préoccupation initiale de séparer au sein de la dynamique les éléments structurels liés à la nature du système et le mouvement. Les éléments structurels sont représentés par la **matrice des masses généralisées M** et le mouvement est caractérisé par la fonction de configuration $\varepsilon(t)$ ainsi que sa dérivée temporelle $\dot{\varepsilon}(t)$. Fidèle à cette démarche de séparation du modèle, on définit la **fonction d'énergie cinétique** qui, à partir d'une configuration $\eta$ et du taux de variation $\hat{\eta}$ de cette configuration exprime l'énergie cinétique du système matériel :

$$ \kappa : \begin{cases} \mathbb{U} \times \hat{\mathbb{U}} \to \mathbb{R} \\ (\eta, \hat{\eta}) \mapsto \kappa(\eta, \hat{\eta}) = \frac{1}{2} \hat{\eta}^T \mathbf{M}(\eta) \hat{\eta} \end{cases} \quad (1.21) $$

# 4 Approche intégrale

J'ai choisi d'utiliser le terme d'approche intégrale pour désigner cette preuve des équations de Lagrange, la première au sens chronologique, car elle s'intéresse à l'aspect global du mouvement, ce que Euler appelait les **causes finales**. On commence par introduire la notion de lagrangien avant de pouvoir énoncer le principe de moindre action. Les équations de Lagrange sont finalement une conséquence de ce principe fondateur de la mécanique. Je détaille alors le cheminement des calculs qui permettent d'arriver à cette conclusion.

## 4.1 Lagrangien

Je commence par rappeler brièvement quelques définitions d'analyse fonctionnelle. Soit $E$ un espace de dimension quelconque, on définit $\Gamma^{t_0, t_1}(E)$ l'ensemble des **courbes paramétrées** $\gamma$ définies sur l'intervalle $[t_0, t_1]$ à valeurs dans $E$. Par abus de notation, l'espace des courbes paramétrées sera noté $\Gamma(E)$, sans référence explicite aux bornes de l'intervalle de définition supposées figées dans l'étude.
On utilise la notation suivante pour désigner un élément de $\Gamma(E)$.

$$ \gamma : \begin{cases} [t_0, t_1] \to E \\ t \mapsto \gamma(t) \end{cases} \qquad (1.22) $$

Une **fonctionnelle** de $E$ est une application définie sur $\Gamma(E)$ à valeurs réelles. Par exemple, la longueur euclidienne d'une courbe paramétrée $\gamma \in \Gamma(\mathbb{R}^n)$ est une fonctionnelle de $\mathbb{R}^n$.

$$ I(\gamma) = \int_{t_0}^{t_1} \sqrt{1 + \dot{\gamma}(t)^2} dt \qquad (1.23) $$

Comme je l'ai mentionné dans le préambule historique, le principe de moindre action repose sur l'idée selon laquelle le mouvement réel d'un système matériel réalise, parmi l'infinité des mouvements virtuels, un minimum de l'intégrale vis-à-vis du temps d'une quantité caractéristique. Cette quantité scalaire $L(\eta, \hat{\eta})$, fonction des coordonnées et des vitesses généralisées, porte le nom de **lagrangien**, en hommage à Joseph Louis Lagrange pour sa contribution décisive dans l'énoncé définitif du principe moindre action. L'expression du Lagrangien peut être justifiée en comparant les équations de Lagrange et la loi fondamentale de la dynamique, qui sont finalement deux principes équivalents pour décrire le mouvement d'un point matériel. Pour cette raison, je donnerai la définition du Lagrangien après avoir présenté les équations de Lagrange.

L'action du lagrangien le long d'un mouvement $\boldsymbol{\varepsilon} \in \Gamma(\mathbb{U})$ est une fonctionnelle définie par :

$$ S(\boldsymbol{\varepsilon}) = \int_{t_0}^{t_1} L(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) dt \qquad (1.24) $$

Cette action, qui synthétise le comportement global du système matériel sur un intervalle de temps $[t_0, t_1]$, figure au centre du principe de moindre action.

## 4.2 Principe de moindre action

Ce principe, qualifié d'équation d'Euler-Lagrange dans certains ouvrages [Dacorogna 1992] a pour définition :

Le mouvement réel $\boldsymbol{\varepsilon} \in \Gamma(\mathbb{U})$ d'un système matériel est un minimum local pour l'action du lagrangien.

$$ S(\boldsymbol{\varepsilon}) = \inf_{\boldsymbol{\varepsilon}^* \in V(\boldsymbol{\varepsilon})} S(\boldsymbol{\varepsilon}^*) \qquad (1.25) $$

$V(\boldsymbol{\varepsilon})$ est un voisinage de l'espace $\Gamma(\mathbb{U})$ autour du mouvement réel $\boldsymbol{\varepsilon}$.

En analogie avec la recherche d'extrema dans $\mathbb{R}^n$ qui s'appuie sur l'étude des racines de la différentielle, on cherche à déterminer l'équivalent de la différentielle de l'action du lagrangien. Comme l'espace des mouvements $\Gamma(\mathbb{U})$ n'a pas une structure d'espace vectoriel normé, les règles classiques du calcul différentiel ne peuvent s'appliquer aux fonctionnelles. Le calcul des
variations introduit par Lagrange consiste en une forme dite faible du calcul différentiel [Fox 1987] [Giaquinta 1996]. L'analogue de la dérivée d'une fonctionnelle est nommée **variation**.

Afin de donner la définition de la variation d'une fonctionnelle, on construit des suites continues $(\varepsilon^*_{\tau})_{\tau \in \mathbb{R}^+}$ de mouvements virtuels convergeant vers le mouvement réel (Figure 6) :

$$ \varepsilon_{\tau}^{*} : \begin{cases} [t_0, t_1] \to \mathbb{U} \\ t \mapsto \varepsilon_{\tau}^{*}(t) \end{cases} \quad \forall \tau \in \mathbb{R}^{+} \quad \text{où} \quad \lim_{\tau \to 0} \varepsilon_{\tau}^{*}(t) = \varepsilon(t) \quad \forall t \in \mathbb{R} \qquad (1.26) $$

<figure><img src="image_15.png" /><figcaption>Figure 6 Famille de mouvements virtuels déduite du mouvement réel.</figcaption></figure>

De chacune de ces familles de mouvements virtuels, on déduit l'expression de la **variation du mouvement** [Bérest 1997] :

$$ \delta \varepsilon^*(t) = \lim_{\tau \to 0} \frac{\varepsilon_{\tau}^*(t) - \varepsilon(t)}{\tau} \qquad (1.27) $$

La variation du mouvement $\delta \varepsilon^*(t)$ déduite d'une suite de mouvements virtuels $(\varepsilon^*_{\tau})_{\tau \in \mathbb{R}^+}$ peut être interprétée comme un déplacement infinitésimal autour de $\varepsilon(t)$ à temps constant (Figure 7). L'usage de la notation $\delta$ proposée par Lagrange dans son calcul des variations témoigne de son souci de clarté en évitant ainsi de confondre la variation du mouvement $\delta \varepsilon^*(t)$ et sa différentielle $d\varepsilon(t)/dt$, tangente au mouvement à l'instant $t$, que l'on a déjà qualifiée de vitesse généralisée.

<figure><img src="image_16.png" /><figcaption>Figure 7 Vitesse généralisée du mouvement vs. variation du mouvement.</figcaption></figure>
La variation $\delta S$ d'une fonctionnelle $S$ est obtenue par passage à la limite des valeurs de $(S(\boldsymbol{\epsilon}_{\tau}^{*}(t)) - S(\boldsymbol{\epsilon}(t)))/\tau$ pour les mouvements virtuels construits en (1.26). De façon synthétique, on peut écrire :

$$ \delta S(\boldsymbol{\epsilon}^*) = \lim_{\tau \to 0} \frac{S(\boldsymbol{\epsilon}_{\tau}^*(t)) - S(\boldsymbol{\epsilon}(t))}{\tau} \qquad (1.28) $$

L'expression de la variation de l'action du lagrangien m'amène à donner un nouvel énoncé du principe de moindre action, communément appelé **forme faible de l'équation d'Euler-Lagrange** :

Le mouvement réel $\boldsymbol{\epsilon} \in \Gamma(U)$ d'un système matériel est tel que sa variation $\delta S(\boldsymbol{\epsilon}^*)$ est nulle pour toute famille de mouvements virtuels $(\boldsymbol{\epsilon}_{\tau}^*)_{\tau \in \mathbb{R}^+}$ satisfaisant la condition (1.26).

$$ \begin{gathered} \forall (\boldsymbol{\epsilon}_{\tau}^*)_{\tau \in \mathbb{R}^+} : \boldsymbol{\epsilon}_{\tau}^* \in V(\boldsymbol{\epsilon}) \quad \forall \tau \in \mathbb{R}^+ \\ \delta S(\boldsymbol{\epsilon}^*) = 0 \end{gathered} \qquad (1.29) $$

## 4.3 Equations de Lagrange

La forme faible du principe de moindre action constitue le point de départ de la construction des équations de Lagrange. Je vais présenter maintenant les différentes étapes du développement théorique qui permettent de déboucher sur la forme classique de ces équations. Pour cela, on commence par explicitier l'expression de la variation de l'action du lagrangien (1.24).

$$ \begin{align} \delta S(\boldsymbol{\epsilon}^*) &= \int_{t_0}^{t_1} \lim_{\tau \to 0} \frac{L(\boldsymbol{\epsilon}_{\tau}^*(t), \dot{\boldsymbol{\epsilon}}_{\tau}^*(t)) - L(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t))}{\tau} dt \tag{1.30} \\ &= \sum_{i=1}^{n} \int_{t_0}^{t_1} \frac{\partial L}{\partial \eta_i}(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) \lim_{\tau \to 0} \frac{\boldsymbol{\epsilon}_{\tau i}^*(t) - \boldsymbol{\epsilon}_i(t)}{\tau} + \frac{\partial L}{\partial \hat{\eta}_i}(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) \lim_{\tau \to 0} \frac{\dot{\boldsymbol{\epsilon}}_{\tau i}^*(t) - \dot{\boldsymbol{\epsilon}}_i(t)}{\tau} dt \nonumber \end{align} $$

Puisque les variables $\tau$ et $t$ sont indépendantes, on a la relation suivante :

$$ \lim_{\tau \to 0} \frac{\dot{\boldsymbol{\epsilon}}_{\tau i}^*(t) - \dot{\boldsymbol{\epsilon}}_i(t)}{\tau} = \lim_{\tau \to 0} \frac{d}{dt} \left( \frac{\boldsymbol{\epsilon}_{\tau i}^*(t) - \boldsymbol{\epsilon}_i(t)}{\tau} \right) = \frac{d}{dt} \lim_{\tau \to 0} \frac{\boldsymbol{\epsilon}_{\tau i}^*(t) - \boldsymbol{\epsilon}_i(t)}{\tau} = \frac{d}{dt} \delta \boldsymbol{\epsilon}_i^*(t) \qquad (1.31) $$

Ainsi, l'équation (1.30) devient

$$ \delta S(\boldsymbol{\epsilon}^*) = \sum_{i=1}^{n} \left\{ \int_{t_0}^{t_1} \frac{\partial L}{\partial \eta_i}(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) \delta \boldsymbol{\epsilon}_i^*(t) dt + \int_{t_0}^{t_1} \frac{\partial L}{\partial \hat{\eta}_i}(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) \frac{d}{dt} \delta \boldsymbol{\epsilon}_i^*(t) dt \right\} \qquad (1.32) $$
Après intégration par parties du deuxième terme apparaissant dans la somme, on obtient :

$$
\delta S(\boldsymbol{\epsilon}^*) = \sum_{i=1}^{n} \left\{ \begin{aligned} & \int_{t_0}^{t_1} \frac{\partial L}{\partial \eta_i}(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) \delta \epsilon_i^*(t) dt \\ & - \int_{t_0}^{t_1} \frac{d}{dt} \left( \frac{\partial L}{\partial \hat{\eta}_i}(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) \right) \delta \epsilon_i^*(t) dt \end{aligned} \right\} \tag{1.33} \\
& + \sum_{i=1}^{n} \left[ \frac{\partial L}{\partial \hat{\eta}_i}(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) \delta \epsilon_i^*(t) \right]_{t_0}^{t_1}
$$

On considère de plus que les mouvements de la famille $(\boldsymbol{\epsilon}_\tau^*)_{\tau \in \mathbb{R}^+}$ appartenant au voisinage $V(\boldsymbol{\epsilon})$
de $\boldsymbol{\epsilon}$ sont coïncidents en $t_0$ et $t_1$ :

$$
\begin{gather}
\forall \tau \in \mathbb{R}^+ : \boldsymbol{\epsilon}_{\tau}^{*}(t_0) = \boldsymbol{\epsilon}(t_0), \boldsymbol{\epsilon}_{\tau}^{*}(t_1) = \boldsymbol{\epsilon}(t_1) \Rightarrow \delta \boldsymbol{\epsilon}^{*}(t_0) = \delta \boldsymbol{\epsilon}^{*}(t_1) = 0 \notag \\
\text{Ainsi} \quad \left[ \frac{\partial L}{\partial \hat{\eta}_i} (\boldsymbol{\epsilon}^{*}(\tau,t), \dot{\boldsymbol{\epsilon}}^{*}(\tau,t)) \delta \epsilon_i^{*}(t) \right]_{t_0}^{t_1} = 0 \quad 1 \le i \le n \tag{1.34}
\end{gather}
$$

On obtient alors l’expression suivante de la variation de l’action du lagrangien :

$$
\delta S(\boldsymbol{\epsilon}^*) = \int_{t_0}^{t_1} \sum_{i=1}^{n} \left( \frac{\partial L}{\partial \eta_i}(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) - \frac{d}{dt} \frac{\partial L}{\partial \hat{\eta}_i}(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) \right) \delta \epsilon_i^*(t) dt \quad (1.35)
$$

La forme faible du principe de moindre action (1.29) me permet d’écrire l’égalité suivante :

$$
\int_{t_0}^{t_1} \sum_{i=1}^{n} \left( \frac{\partial L}{\partial \eta_i}(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) - \frac{d}{dt} \frac{\partial L}{\partial \hat{\eta}_i}(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) \right) \delta \epsilon_i^*(t) dt = 0 \quad (1.36)
$$

L'égalité (1.36) est valide pour une variation quelconque du mouvement (i.e. pour une famille quelconque de mouvements virtuels). Elle est donc aussi valide pour les *n* variations unitaires du mouvement, déduites des *n* familles $(\boldsymbol{\epsilon}_{\tau}^{*j})_{\tau \in \mathbb{R}^{+}}$ ($1 \le j \le n$) de mouvements, dont les coordonnées sont toutes nulles sauf une (1.37).

$$
\boldsymbol{\epsilon}_{\tau}^{*j} : \begin{cases} [t_0, t_1] \to \mathbb{U} \\ (t) \mapsto \boldsymbol{\epsilon}_{\tau}^{*j}(t) \end{cases} \quad \forall \tau \in \mathbb{R}^{+} \quad \text{pour } 1 \le j \le n \tag{1.37} \\
\text{avec } \forall t \in [t_0, t_1] : \delta \boldsymbol{\epsilon}_i^{*j}(t) = \delta_{ij}
$$

L'utilisation de ces *n* variations permet d'aboutir à un système de *n* équations que l'on
nomme **équations de Lagrange** :

$$
\frac{\partial L}{\partial \eta_i}(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) - \frac{d}{dt} \frac{\partial L}{\partial \hat{\eta}_i}(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) = 0 \quad 1 \le i \le n \quad (1.38)
$$
Dans la littérature contemporaine [Germain 1986] [Arnold 1989] [Calkin 1996], il est plus courant de trouver cette expression équivalente des équations de Lagrange :

$$ \frac{d}{dt} \frac{\partial L}{\partial \dot{\eta}_i}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) - \frac{\partial L}{\partial \eta_i}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) = 0 \quad 1 \le i \le n \quad (1.39) $$

## 4.4 Expression du lagrangien

Il est temps maintenant de donner la définition du lagrangien. Pour cela, il faut considérer le formalisme lagrangien comme une extension de la mécanique à l'étude du mouvement de systèmes matériels autres que des points matériels. Ainsi, les équations de Lagrange doivent être équivalentes à la loi fondamentale de la dynamique dans le cas d'un système matériel réduit à un point.

Je commence par rappeler la définition de la loi fondamentale de la dynamique qui conditionne l'évolution d'un point matériel de masse *m*, de vitesse $\dot{\mathbf{x}}(t)$ soumis à une force $\mathbf{F}(\mathbf{x})$ :

$$ \frac{d}{dt}(m\dot{\mathbf{x}}(t)) = \mathbf{F}(\mathbf{x}(t)) \quad (1.40) $$

L'équation (1.40) est équivalente à un système de 3 équations scalaires.

$$ \frac{d}{dt}(m\dot{\mathbf{x}}_i(t)) = F_i(\mathbf{x}(t)) \quad 1 \le i \le 3 \quad (1.41) $$

Si l'on considère que la force $\mathbf{F}(\mathbf{x})$ dérive d'un potentiel $U$ et que l'on se remémore la définition de la fonction d'énergie cinétique (1.21), on obtient une nouvelle formulation du système d'équations de (1.41).

$$ \frac{d}{dt} \frac{\partial \kappa}{\partial \dot{x}_i}(\mathbf{x}(t), \dot{\mathbf{x}}(t)) = - \frac{\partial U}{\partial x_i}(\mathbf{x}(t)) \quad 1 \le i \le 3 \quad (1.42) $$

Puisque l'énergie potentielle $U$ ne dépend pas explicitement de la vitesse de la particule et que l'énergie cinétique quant à elle est uniquement fonction de la vitesse,

$$ \frac{\partial U}{\partial \dot{x}_i}(\mathbf{x}(t)) = 0 \quad \text{et} \quad \frac{\partial \kappa}{\partial x_i}(\mathbf{x}(t), \dot{\mathbf{x}}(t)) = 0 \quad 1 \le i \le 3 \quad (1.43) $$

on débouche sur l'expression suivante:

$$ \frac{d}{dt} \frac{\partial}{\partial \dot{x}_i} \left( \kappa(\mathbf{x}(t), \dot{\mathbf{x}}(t)) - U(\mathbf{x}(t)) \right) - \frac{\partial}{\partial x_i} \left( \kappa(\mathbf{x}(t), \dot{\mathbf{x}}(t)) - U(\mathbf{x}(t)) \right) = 0 \quad 1 \le i \le 3 \quad (1.44) $$

Conformément à l'équation (1.39), le lagrangien d'un point matériel soumis à des efforts conservatifs dérivant d'un potentiel $U$ correspond à la différence de l'énergie cinétique $\kappa(\mathbf{x}(t), \dot{\mathbf{x}}(t))$ et de l'énergie potentielle $U(\mathbf{x}(t))$.
En exploitant la propriété d'additivité particulaire des fonctions d'énergie cinétique $\kappa(\mathbf{x}(t), \dot{\mathbf{x}}(t))$ et potentielle $U(\mathbf{x}(t))$ (c'est à dire que l'énergie totale d'un ensemble quelconque de particules est la somme des énergies de chaque particule), le lagrangien d'un système matériel quelconque est défini à partir de la différence de sa fonction d'énergie cinétique et de sa fonction d'énergie potentielle.

$$L(\eta, \dot{\eta}) = \kappa(\eta, \dot{\eta}) - U(\eta) = \frac{1}{2} \dot{\eta}^T M(\eta) \dot{\eta} - U(\eta) \quad (1.45)$$

Le lagrangien ne correspond pas à l'énergie totale du système puisque pour des efforts conservatifs cette quantité est constante. Le lagrangien doit plutôt être perçu comme une mesure des échanges entre énergie potentielle et énergie cinétique. De cette manière, le principe de moindre action postule que le mouvement réel est le mouvement qui tend à minimiser les transferts d'énergie au sein du système matériel.

En substituant l'expression du lagrangien proposée en (1.45) dans l'équation (1.39), on obtient une formulation des équations de Lagrange conforme à celle couramment utilisée [Germain 1986].

$$\frac{d}{dt} \frac{\partial \kappa}{\partial \dot{\eta}_i}(\varepsilon(t), \dot{\varepsilon}(t)) - \frac{\partial \kappa}{\partial \eta_i}(\varepsilon(t), \dot{\varepsilon}(t)) = - \frac{\partial U}{\partial \eta_i}(\varepsilon(t)) \quad 1 \le i \le n \quad (1.46)$$

Cette première justification des équations de Lagrange est restreinte au cas du mouvement de systèmes matériels soumis à des efforts conservatifs qui apparaissent dans le formalisme lagrangien par le biais du potentiel dont ils dérivent. Cette restriction sera levée grâce à la deuxième justification des équations de Lagrange présentée dans la section suivante.

# 5 Approche différentielle

Dans cette partie, je propose une deuxième justification des équations de Lagrange fidèle à celle présentée dans la **Mécanique Analytique** publiée en 1788. Je qualifie cette approche de différentielle car on privilégie ici une étude locale en temps du mouvement, cherchant à en déterminer les variations infinitésimales. La notion de vitesse virtuelle ayant atteint sa maturité, l'énoncé du **principe des puissances virtuelles** s'impose comme principe fondateur et constitue le nouveau point de départ de la construction des équations de Lagrange. Ce principe énoncé dans le cadre de la statique fut étendu au cas de la dynamique grâce à la contribution de d'Alembert.

## 5.1 Mouvements virtuels

On considère le mouvement $\varepsilon \in \Gamma(U)$ comme le mouvement réel du système matériel, décrivant l'évolution de la configuration géométrique de celui-ci. Au sein de l'espace des configurations étendues $U$, on peut construire une infinité de mouvements autres que $\varepsilon$, représentant tous des évolutions géométriques accessibles du système matériel compatibles avec le paramétrage choisi. Ces mouvements que j'ai déjà introduits pour les besoins de l'approche intégrale sont appelés **mouvements virtuels compatibles**.
Puisque l'on focalise l'étude du mouvement du système matériel à un instant donné $t'$, on s'intéresse plus particulièrement à une sous famille de mouvements virtuels $\boldsymbol{\varepsilon}^*$ définis comme suit :

$$ \boldsymbol{\varepsilon}^* : \begin{cases} \mathbb{R} \to \mathbb{U} \\ t \mapsto \boldsymbol{\varepsilon}^*(t) \end{cases} \quad \text{avec } \boldsymbol{\varepsilon}^*(t') = \boldsymbol{\varepsilon}(t') \qquad (1.47) $$

Ces mouvements sont appelés les **mouvements virtuels coïncidents en t'** (i.e. coïncidant à l'instant $t'$ avec le mouvement réel). De tels mouvements ont la particularité de conférer au système matériel une configuration identique à celle fournie par le mouvement réel à l'instant donné $t'$.

Tous ces mouvements étant confondus en $t'$, la seule caractéristique qui les distingue localement réside dans leur cinématique, c'est à dire l'allure de leurs dérivées temporelles $\dot{\boldsymbol{\varepsilon}}^*(t')$ à l'instant $t'$. Ainsi, on introduit de nouvelles familles de mouvements virtuels dont la cinématique vérifie certaines conditions. Parmi ces familles, on distingue $n$ classes d'équivalence $(\mathcal{E}^j)_{1\le j\le n}$ distinctes de mouvements virtuels dont les éléments ont des évolutions locales identiques remarquables (1.48). On choisit arbitrairement dans chaque famille $\mathcal{E}^j$ un représentant $\boldsymbol{\varepsilon}^*_{\mathcal{E}^j} (1\le j\le n)$ qui constituent les **n mouvements figés d'entraînement**.

$$ \forall \boldsymbol{\varepsilon}^* \in \mathcal{E}^j: \boldsymbol{\varepsilon}^*(t') = \boldsymbol{\varepsilon}(t') \quad \text{et} \quad \boldsymbol{\varepsilon}_i^*(t') = \delta_{ij} \quad \text{pour } 1\le j\le n \qquad (1.48) $$

Localement, les mouvements virtuels figés d'entraînement évoluent respectivement le long d'une seule coordonnée généralisée $\eta_i$, et ce à vitesse unitaire. Par la suite, bien que l'on s'intéresse à l'évolution du mouvement à l'instant $t'$, on utilisera la lettre $t$ pour désigner l'instant d'étude afin de ne pas alourdir les notations

## 5.2 Vitesse et puissance virtuelles

La **vitesse virtuelle** $v^*$ d'une particule $\omega$ dans un mouvement virtuel $\boldsymbol{\varepsilon}^*$ coïncident est déduite de la substitution de la vitesse généralisée $\dot{\boldsymbol{\varepsilon}}(t)$ par la vitesse virtuelle généralisée $\dot{\boldsymbol{\varepsilon}}^*(t)$ dans l'expression de la fonction de vitesse, et ce à temps fixe (i.e. $\dot{\boldsymbol{\varepsilon}}_0^*(t) = 0$).

$$ v^*(t, \omega) = \hat{\Phi}(\boldsymbol{\varepsilon}^*(t), \dot{\boldsymbol{\varepsilon}}^*(t), \omega) = \sum_{i=1}^{n} \frac{\partial \Phi}{\partial \eta_i}(\boldsymbol{\varepsilon}^*(t), \omega) \boldsymbol{\varepsilon}_i^*(t) = \sum_{i=1}^{n} \frac{\partial \Phi}{\partial \eta_i}(\boldsymbol{\varepsilon}(t), \omega) \dot{\boldsymbol{\varepsilon}}_i^*(t) \qquad (1.49) $$

La **puissance virtuelle** $P^{*F}$ d'un effort $F$ dans un mouvement virtuel $\boldsymbol{\varepsilon}^*$ est un nombre scalaire obtenu par produit scalaire du vecteur d'effort avec la vitesse virtuelle $v^*(t, \omega)$ du point d'application identifié par $\omega$.

Pour un effort ponctuel $F(\eta, \hat{\eta})$ dont le point d'application est identifié par $\omega$, on a :

$$ \begin{aligned} P^{*F} &= F(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}^*(t)) \cdot v^*(t, \omega) \\ &= \sum_{i=1}^{n} F(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}^*(t)) \cdot \left\{ \frac{\partial \Phi}{\partial \eta_i}(\boldsymbol{\varepsilon}(t), \omega) \boldsymbol{\varepsilon}_i^*(t) \right\} \end{aligned} \qquad (1.50) $$
On définit alors la puissance virtuelle $P^*F$ comme une fonction des coordonnées généralisées $\eta$ et de leurs taux de variation réels $\hat{\eta}$ et virtuels $\hat{\eta}^*$.

$$ P^*F(\eta, \hat{\eta}, \hat{\eta}^*) = \sum_{i=1}^{n} Q_i(\eta, \hat{\eta}) \hat{\eta}_i^* \qquad (1.51) \\ \text{avec } Q_i(\eta, \hat{\eta}) = F(\eta, \hat{\eta}) \cdot \frac{\partial \Phi}{\partial \eta_i}(\eta, \omega) $$

Pour un champ de forces $F(\eta, \hat{\eta}, \omega)$ de domaine d'application $\Omega^F$, on a :

$$ P^{*F} = \int_{\Omega^F} F(\varepsilon(t), \dot{\varepsilon}^*(t), \omega) \cdot v^*(t, \omega) d\omega \\ = \sum_{i=1}^{n} \int_{\Omega^F} F(\varepsilon(t), \dot{\varepsilon}^*(t), \omega) \cdot \left\{ \frac{\partial \Phi}{\partial \eta_i}(\varepsilon(t), \omega) \varepsilon_i^*(t) \right\} d\omega \qquad (1.52) $$

Tout comme en (1.51), on modifie l'expression de la puissance virtuelle de la manière suivante.

$$ P^{*F}(\eta, \hat{\eta}, \hat{\eta}^*) = \sum_{i=1}^{n} Q_i(\eta, \hat{\eta}) \hat{\eta}_i^* \qquad (1.53) \\ \text{avec } Q_i(\eta, \hat{\eta}) = \int_{\Omega^F} F(\eta, \hat{\eta}, \omega) \cdot \frac{\partial \Phi}{\partial \eta_i}(\eta, \omega) d\omega $$

Le vecteur $Q(\eta, \hat{\eta})$ composé des $n$ termes $Q_i(\eta, \hat{\eta})$ est communément appelé vecteur des **forces généralisées** appliquées au système matériel.

Dans le cas particulier d'efforts conservatifs, on a une expression particulière du vecteur des forces généralisées. On suppose l'existence d'un potentiel réel dont le champ de force dérive.

$$ \exists U : \mathbb{R}^3 \to \mathbb{R} : F(\eta, \omega) = -\operatorname{grad} U(\Phi(\eta, \omega)) \qquad (1.54) $$

En introduisant une quantité $\tilde{U}$ qui dépend explicitement de la configuration cinématique du système (coordonnées et vitesses généralisées),

$$ \tilde{U} : \begin{cases} \mathbb{U} \times \Omega \to \mathbb{R} \\ (\eta, \omega) \mapsto \tilde{U}(\eta, \omega) = U(\Phi(\eta, \omega)) \end{cases} \qquad (1.55) $$

l'expression de la puissance virtuelle pour cette densité d'efforts particulière devient :

$$ P^{*F} = \int_{\Omega^F} \frac{-\partial U}{\partial x}(\Phi(\varepsilon(t), \omega)) \cdot v^*(t, \omega) d\omega \\ = \sum_{i=1}^{n} \left\{ \int_{\Omega^F} \frac{-\partial U}{\partial x}(\Phi(\varepsilon(t), \omega)) \cdot \frac{\partial \Phi}{\partial \eta_i}(\varepsilon(t), \omega) d\omega \right\} \dot{\varepsilon}_i^*(t) \qquad (1.56) \\ = \sum_{i=1}^{n} \left\{ \int_{\Omega^F} \frac{-\partial \tilde{U}}{\partial \eta_i}(\varepsilon(t), \omega) d\omega \right\} \dot{\varepsilon}_i^*(t) $$
En définitive, on trouve :

$$P^{*F}(\eta, \hat{\eta}, \hat{\eta}^*) = \sum_{i=1}^{n} Q_i(\eta, \hat{\eta}) \hat{\eta}_i^* \quad \text{avec} \quad Q_i(\eta, \hat{\eta}) = \int_{\Omega^F} \frac{-\partial \tilde{U}}{\partial \eta_i}(\eta, \omega) d\omega \quad (1.57)$$

## 5.3 Principe des puissances virtuelles

Maintenant que j'ai donné la définition de la puissance virtuelle, je suis en mesure d'énoncer le premier **principe des puissances virtuelles**. Ce principe qui consiste en une équation bilan n'est valide que dans le cadre de la statique.

### Principe statique

Dans un référentiel galiléen, la somme des puissances virtuelles des efforts appliqués à un système matériel au repos est identiquement nulle. Ou de façon plus synthétique :

$$P^{*F} = 0 \quad (1.58)$$

où le vecteur F est la résultante de l'ensemble des forces appliquées au système.

## 5.4 Principe de d'Alembert

L'apport de d'Alembert consiste en l'extension du principe statique des puissances virtuelles à la dynamique. Il introduit pour cela la notion de quantité d'accélération qui intervient dans son principe par le biais de sa puissance virtuelle définie en (1.59).

$$ 
\begin{align} 
P^{*\gamma} &= \int_{\Omega} \frac{d}{dt} (\rho(\omega, t)v(t, \omega)) \cdot v^{*}(t, \omega) d\omega \tag{1.59} \\ 
&= \sum_{i=0}^{n} \int_{\Omega} \frac{d}{dt} (\rho(\omega, t) \hat{\Phi}(\varepsilon(t), \dot{\varepsilon}(t), \omega)) \cdot \left\{ \frac{\partial \Phi}{\partial \eta_i}(\varepsilon(t), \omega) \dot{\varepsilon}_i^*(t) \right\} d\omega \nonumber 
\end{align} 
$$

En utilisant les notations introduites pour exprimer les puissances virtuelles des efforts, on a :

$$ 
\begin{gather} 
P^{*\gamma}(\eta, \hat{\eta}, \hat{\eta}^*) = \sum_{i=1}^{n} P_i^{*\gamma}(\eta, \hat{\eta}) \hat{\eta}_i^* \tag{1.60} \\ 
\text{avec } P_i^{*\gamma}(\eta, \hat{\eta}) = \int_{\Omega} \frac{d}{dt} (\rho(\omega, \eta_0) \hat{\Phi}(\eta, \hat{\eta}, \omega)) \cdot \frac{\partial \Phi}{\partial \eta_i}(\eta, \omega) d\omega \nonumber 
\end{gather} 
$$

le principe de d'Alembert devient alors :

### Principe dynamique

Dans un référentiel galiléen, la puissance virtuelle de l'ensemble des efforts appliqués au système est égale à la puissance virtuelle des quantités d'accélération de celui-ci. Cet énoncé est résumé par l'équation (1.61).

$$P^{*\gamma} = P^{*F} \quad (1.61)$$
## 5.5 Equations de Lagrange

Le principe de d'Alembert constitue un autre principe fondateur dont la nature est très différente du principe de moindre action. Mais il apparaît clairement que le champ d'application du principe de d'Alembert est plus vaste que celui du principe de moindre action dans la mesure où les efforts, caractérisés par leurs puissances virtuelles, sont tout à fait quelconques (contrairement au principe de moindre action, qui ne prend en compte que des efforts conservatifs). Ce point fut sûrement crucial dans la décision de Lagrange d'adopter le principe de d'Alembert comme principe fondateur de la mécanique analytique. Je vais maintenant justifier l'équivalence du principe de d'Alembert avec les équations de Lagrange. Je commence par rappeler la définition de la fonction énergie cinétique.

$$ \kappa(\varepsilon(t), \dot{\varepsilon}(t)) = \frac{1}{2} \int_{\Omega} \rho(\omega, t) \hat{\Phi}(\varepsilon(t), \dot{\varepsilon}(t), \omega)^2 d\omega \quad (1.62) $$

Ainsi,

$$ \begin{aligned} \frac{\partial \kappa}{\partial \hat{\eta}_i}(\varepsilon(t), \dot{\varepsilon}(t)) &= \frac{\partial}{\partial \hat{\eta}_i} \left( \frac{1}{2} \int_{\Omega} \rho(\omega, t) \left( \hat{\Phi}(\varepsilon(t), \dot{\varepsilon}(t), \omega) \right)^2 d\omega \right) \\ &= \frac{1}{2} \int_{\Omega} \rho(\omega, t) \left( 2 \hat{\Phi}(\varepsilon(t), \dot{\varepsilon}(t), \omega) \cdot \frac{\partial \hat{\Phi}}{\partial \hat{\eta}_i}(\varepsilon(t), \dot{\varepsilon}(t), \omega) \right) d\omega \end{aligned} \quad (1.63) $$

En exploitant la définition de la fonction de vitesse (1.14), on obtient :

$$ \frac{\partial \kappa}{\partial \hat{\eta}_i}(\varepsilon(t), \dot{\varepsilon}(t)) = \int_{\Omega} \rho(\omega, t) \hat{\Phi}(\varepsilon(t), \dot{\varepsilon}(t), \omega) \cdot \frac{\partial \Phi}{\partial \eta_i}(\varepsilon(t), \omega) d\omega \quad (1.64) $$

Puis, en dérivant l'équation (1.64) par rapport au temps, on débouche sur l'expression suivante :

$$ \begin{aligned} \frac{d}{dt} \frac{\partial \kappa}{\partial \hat{\eta}_i}(\varepsilon(t), \dot{\varepsilon}(t)) = & \int_{\Omega} \frac{d}{dt} \left( \rho(\omega, t) \hat{\Phi}(\varepsilon(t), \dot{\varepsilon}(t), \omega) \right) \cdot \frac{\partial \Phi}{\partial \eta_i}(\varepsilon(t), \omega) d\omega \\ & + \int_{\Omega} \rho(\omega, t) \hat{\Phi}(\varepsilon(t), \dot{\varepsilon}(t), \omega) \cdot \frac{d}{dt} \frac{\partial \Phi}{\partial \eta_i}(\varepsilon(t), \omega) d\omega \end{aligned} \quad (1.65) $$

Le premier terme intégral du membre de droite est constitué de la force généralisée $P_i^*(\varepsilon(t), \dot{\varepsilon}(t))$ issue de la puissance virtuelle des quantités d'accélération (1.59).

On s'intéresse maintenant au deuxième terme intégral de l'équation (1.65). Pour en modifier l'expression, j'ai besoin d'utiliser le lemme suivant (démontré en bas de page<fnref n="3" />).

$$ \frac{d}{dt} \frac{\partial \Phi}{\partial \eta_i}(\varepsilon(t), \omega) = \frac{\partial \hat{\Phi}}{\partial \eta_i}(\varepsilon(t), \dot{\varepsilon}(t), \omega) \quad (1.66) $$

<footnote n="3">³ $\frac{d}{dt} \frac{\partial \Phi}{\partial \eta_i}(\varepsilon(t), \omega) = \sum_{j=1}^{n} \frac{\partial^2 \Phi}{\partial \eta_j \partial \eta_i} (\varepsilon(t), \omega) \dot{\varepsilon}_j(t) = \frac{\partial}{\partial \eta_i} \sum_{j=1}^{n} \frac{\partial \Phi}{\partial \eta_j} (\varepsilon(t), \omega) \dot{\varepsilon}_j(t) = \frac{\partial \hat{\Phi}}{\partial \eta_i} (\varepsilon(t), \dot{\varepsilon}(t), \omega)$</footnote>
L'équation devient alors, grâce à l'égalité (1.65) et le lemme (1.66) :

$$
\begin{equation}
\begin{split}
\frac{d}{dt} \frac{\partial \kappa}{\partial \hat{\eta}_i} (\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) ={}& P_i^{*\gamma}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) \\
& + \int_{\Omega} \rho(\omega, t) \hat{\Phi}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t), \omega) \cdot \frac{\partial \hat{\Phi}}{\partial \eta_i}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t), \omega) d\omega
\end{split}
\tag{1.67}
\end{equation}
$$

Mais le terme intégral de l'équation (1.67) satisfait une autre égalité. En effet,

$$
\frac{\partial \kappa}{\partial \eta_i}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) = \int_{\Omega} \rho(\omega, t) \hat{\Phi}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t), \omega) \cdot \frac{\partial \hat{\Phi}}{\partial \eta_i}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t), \omega) d\omega \quad (1.68)
$$

Les équations (1.67) et (1.68) m'amènent à écrire :

$$
P_i^{*\gamma}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) = \frac{d}{dt} \frac{\partial \kappa}{\partial \hat{\eta}_i}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) - \frac{\partial \kappa}{\partial \eta_i}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) \quad (1.69)
$$

La puissance virtuelle $P_i^{*\gamma}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t))$ des quantités d'accélération devient alors :

$$
P_i^{*\gamma}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) = \sum_{i=1}^{n} \left\{ \frac{d}{dt} \frac{\partial \kappa}{\partial \hat{\eta}_i}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) - \frac{\partial \kappa}{\partial \eta_i}(\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) \right\} \dot{\boldsymbol{\varepsilon}}_i^*(t) \quad (1.70)
$$

Conformément au principe de d'Alembert (1.61) et en utilisant les vitesses virtuelles
généralisées des *n* mouvements figés d'entraînement, on obtient les *n* équations de Lagrange.

$$
\left\{ \frac{d}{dt} \frac{\partial \kappa}{\partial \hat{\eta}_i} (\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) - \frac{\partial \kappa}{\partial \eta_i} (\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) \right\} Q_i (\boldsymbol{\varepsilon}(t), \dot{\boldsymbol{\varepsilon}}(t)) \quad 1 \le i \le n \quad (1.71)
$$

5.6 Prise en compte des contraintes

Dans le formalisme lagrangien, les contraintes peuvent être incorporées naturellement aux équations (1.71) par le biais de la puissance virtuelle des efforts nécessaires à leur maintien. Je limite délibérément mon étude à la réalisation des contraintes holonomes, cinématiques et unilatérales. En effet, la prise en compte des contraintes constitue un vaste sujet qui requiert sûrement autant de pages que celles déjà consacrées aux équations de Lagrange. Pour rester fidèle à la démarche adoptée dans cette étude bibliographique, ce paragraphe a pour vocation de donner une interprétation géométrique et si possible intuitive de l'expression des contraintes en formalisme lagrangien.

5.6.1 Contraintes holonomes

Soit une contrainte holonome de la forme :

$$
\varphi(\eta) = 0 \tag{1.72}
$$
Tout mouvement virtuel compatible $\boldsymbol{\varepsilon}^*(t)$ coïncident avec le mouvement réel $\boldsymbol{\varepsilon}(t)$ vérifie les égalités suivantes :

$$ \varphi(\boldsymbol{\varepsilon}^*(t)) = 0 \Rightarrow \frac{d\varphi}{dt}(\boldsymbol{\varepsilon}^*(t)) = \sum_{i=1}^{n} \frac{\partial \varphi}{\partial \eta_i}(\boldsymbol{\varepsilon}(t)) \dot{\boldsymbol{\varepsilon}}_i^*(t) = 0 \quad (1.73) $$

L'équation (1.73) est une relation d'orthogonalité entre le vecteur $\partial_{\eta} \varphi$ et les vitesses virtuelles généralisées $\dot{\boldsymbol{\varepsilon}}^*(t)$ des mouvements compatibles avec la contrainte holonome (1.72). Par conséquent, ces vitesses virtuelles généralisées appartiennent toutes à l'hyperplan de normale $\partial_{\eta} \varphi$ (Figure 8).

<figure><img src="image_17.png" /><figcaption>Figure 8 Interprétation géométrique de l'orthogonalité au plan tangent de l'hyper surface issue d'une contrainte holonome dans l'espace des configurations, projeté de l'espace des configurations étendues pour $t=cte$.</figcaption></figure>

D'autre part puisque les liaisons sont considérées comme parfaites, elles respectent la **loi de comportement des liaisons parfaites** [Germain 1986] qui stipule que la puissance virtuelle des efforts issus de ces contraintes est positive ou nulle.

$$ \sum_{i=1}^{n} Q_{i}^{\varphi}(\boldsymbol{\varepsilon}(t)) \dot{\boldsymbol{\varepsilon}}_{i}^{*}(t) \geq 0 \quad (1.74) $$

Le vecteur $Q^\varphi$ est le vecteur des forces généralisées garantissant la réalisation de la contrainte holonome $\varphi$. L'inégalité (1.74) est vérifiée pour tout mouvement virtuel. Ainsi, si l'on considère deux mouvements virtuels compatibles dont les vitesses virtuelles généralisées sont opposées dans l'hyperplan tangent, on arrive à la conclusion suivante :

$$ \sum_{i=1}^{n} Q_{i}^{\varphi}(\boldsymbol{\varepsilon}(t)) \dot{\boldsymbol{\varepsilon}}_{i}^{*}(t) = 0 \quad (1.75) $$
Les équations (1.73) et (1.75) permettent de déduire :

$$ \frac{\partial \varphi}{\partial \boldsymbol{\eta}} \perp \dot{\boldsymbol{\epsilon}}^* \text{ et } \mathbf{Q}^\varphi \perp \dot{\boldsymbol{\epsilon}}^* \qquad (1.76) \\ \forall \dot{\boldsymbol{\epsilon}}^* \in T_{\varepsilon(t)}(\mathbb{U}^\varphi) \quad \text{où} \quad \mathbb{U}^\varphi = \{\boldsymbol{\eta} \in \mathbb{U} : \varphi(\boldsymbol{\eta}) = 0\} $$

où l'ensemble $T_{\varepsilon(t)}(\mathbb{U}^\varphi)$ correspond à l'**espace tangent** au sous-espace $\mathbb{U}^\varphi$ en $\varepsilon(t)$.

D'un point de vue géométrique (Figure 8), une contrainte scalaire portant sur les coordonnées généralisées $\boldsymbol{\eta}$ est une équation implicite, au sens de la géométrie différentielle. Elle restreint donc l'évolution des configurations à une hyper surface de dimension $(n+1)-1=n$ de l'espace des configurations étendues $\mathbb{U}$. L'espace tangent $T_{\varepsilon(t)}(\mathbb{U}^\varphi)$ à l'hyper surface $\mathbb{U}^\varphi$ en $\varepsilon(t)$ est de dimension $n$ aussi et a naturellement pour vecteur normal, le vecteur $\partial_{\boldsymbol{\eta}}\varphi$.

l'espace orthogonal à $T_{\varepsilon(t)}(\mathbb{U}^\varphi)$ est appelé l'**espace cotangent** au sous-espace $\mathbb{U}^\varphi$ en $\varepsilon(t)$. On le note $\mathrm{coT}_{\varepsilon(t)}(\mathbb{U}^\varphi)$. Les espaces tangent et cotangent en $\varepsilon(t)$ étant supplémentaires,

$$ \mathbb{R}^{n+1} = T_{\varepsilon(t)}(\mathbb{U}^\varphi) \oplus \mathrm{coT}_{\varepsilon(t)}(\mathbb{U}^\varphi) \qquad (1.77) $$

on a la relation suivante entre les dimensions de chaque espace :

$$ \dim \mathrm{coT}_{\varepsilon(t)}(\mathbb{U}^\varphi) = \dim \mathbb{R}^{n+1} - \dim T_{\varepsilon(t)}(\mathbb{U}^\varphi) = (n+1) - n = 1 \qquad (1.78) $$

L'espace cotangent $\mathrm{coT}_{\varepsilon(t)}(\mathbb{U}^\varphi)$ en $\varepsilon(t)$ associé à l'hyper surface $\mathbb{U}^\varphi$ est donc de dimension 1, réduit à une droite vectorielle. Par conséquent, tous les vecteurs de $\mathrm{coT}_{\varepsilon(t)}(\mathbb{U}^\varphi)$ sont colinéaires. En s'appuyant sur (1.76), on a :

$$ \frac{\partial \varphi}{\partial \boldsymbol{\eta}}(\boldsymbol{\varepsilon}(t)), \mathbf{Q}^{\varphi}(\boldsymbol{\varepsilon}(t)) \in \mathrm{coT}_{\varepsilon(t)}(\mathbb{U}^{\varphi}) \Rightarrow \exists \lambda(t) \in \mathbb{R}: \mathbf{Q}^{\varphi}(\boldsymbol{\varepsilon}(t)) = \lambda(t) \frac{\partial \varphi}{\partial \boldsymbol{\eta}}(\boldsymbol{\varepsilon}(t)) \qquad (1.79) $$

### 5.6.2 Contraintes cinématiques

Soit une contrainte cinématique de la forme :

$$ \psi(\boldsymbol{\eta}, \hat{\boldsymbol{\eta}}) = \mathbf{a}(\boldsymbol{\eta}) \cdot \hat{\boldsymbol{\eta}} = 0 \qquad (1.80) $$

Tout mouvement virtuel compatible $\boldsymbol{\varepsilon}^*(t)$ coïncident avec le mouvement réel $\boldsymbol{\varepsilon}(t)$ vérifie aussi l'égalité (1.80).

$$ \psi(\boldsymbol{\varepsilon}^*(t), \dot{\boldsymbol{\varepsilon}}^*(t)) = \mathbf{a}(\boldsymbol{\varepsilon}^*(t)) \cdot \dot{\boldsymbol{\varepsilon}}^*(t) = \mathbf{a}(\boldsymbol{\varepsilon}(t)) \cdot \dot{\boldsymbol{\varepsilon}}^*(t) = 0 \qquad (1.81) $$

L'équation (1.81) est, tout comme l'équation (1.73), une relation d'orthogonalité contraignant les vitesses virtuelles généralisées $\dot{\boldsymbol{\varepsilon}}^*(t)$ à appartenir à un hyperplan $\Pi(\boldsymbol{\varepsilon}(t))$ de normale $\mathbf{a}(\boldsymbol{\varepsilon}(t))$.
En supposant les liaisons parfaites et en considérant des vitesses virtuelles généralisées opposées dans $\Pi(\boldsymbol{\varepsilon}(t))$, on obtient la deuxième égalité suivante,

$$ \sum_{i=1}^{n} Q_{i}^{\varphi}(\boldsymbol{\varepsilon}(t)) \dot{\boldsymbol{\varepsilon}}_{i}^{*}(t) = 0 \qquad (1.82) $$

Les équations (1.81) et (1.82) permettent de déduire :

$$ \mathbf{a} \perp \dot{\boldsymbol{\varepsilon}}^* \quad \text{et} \quad \mathbf{Q}^\varphi \perp \dot{\boldsymbol{\varepsilon}}^* \quad \forall \dot{\boldsymbol{\varepsilon}}^* \in \Pi(\boldsymbol{\varepsilon}(t)) \qquad (1.83) $$

En reproduisant le raisonnement utilisé pour les contraintes holonomes, on aboutit à une conclusion similaire, à savoir :

$$ \mathbf{a}(\boldsymbol{\varepsilon}(t)), \mathbf{Q}^{\varphi}(\boldsymbol{\varepsilon}(t)) \in \mathrm{coT}_{\varepsilon(t)}(\mathbb{U}^{\varphi}) \Rightarrow \exists \lambda(t) \in \mathbb{R} : \mathbf{Q}^{\varphi}(\boldsymbol{\varepsilon}(t)) = \lambda(t)\mathbf{a}(\boldsymbol{\varepsilon}(t)) \qquad (1.84) $$

### 5.6.3 Contraintes unilatérales

Soit une contrainte unilatérale de la forme

$$ \varphi(\eta) \ge 0 \qquad (1.85) $$

En suivante la démarche utilisée dans (1.73), l'équation (1.85) permet de déduire la relation suivante :

$$ \varphi(\boldsymbol{\varepsilon}^*(t)) \ge 0 \quad \Rightarrow \quad \frac{d\varphi}{dt}(\boldsymbol{\varepsilon}^*(t)) = \sum_{i=1}^{n} \frac{\partial \varphi}{\partial \eta_i}(\boldsymbol{\varepsilon}(t)) \dot{\boldsymbol{\varepsilon}}_i^*(t) \ge 0 \qquad (1.86) $$

On restreint alors l'ensemble des mouvements virtuels considérés aux mouvements qui garantissent la positivité de l'inéquation (1.86). A partir de ce moment, deux cas de figure se présentent :

D'une part, la contrainte n'est pas réalisée (i.e. $\varphi(\boldsymbol{\varepsilon}^*(t)) > 0$). Toutes les évolution cinématiques sont donc envisageables et conformément à la loi de comportement des liaisons parfaites (1.75), on a la relation suivante :

$$ \sum_{i=1}^{n} Q_{i}^{\varphi}(\boldsymbol{\varepsilon}(t)) \dot{\boldsymbol{\varepsilon}}_{i}^{*}(t) = 0 \quad \forall \dot{\boldsymbol{\varepsilon}}^{*} \in \mathbb{R}^{n+1} \quad \Rightarrow \quad \mathbf{Q}^{\varphi} \perp \dot{\boldsymbol{\varepsilon}}^{*} \quad \forall \dot{\boldsymbol{\varepsilon}}^{*} \in \mathbb{R}^{n+1} \qquad (1.87) $$

Les efforts $Q_i^\varphi(\boldsymbol{\varepsilon}(t))$ relatifs à cette contrainte sont donc identiquement nuls. D'autre part, la contrainte est réalisée (i.e. $\varphi(\boldsymbol{\varepsilon}^*(t)) = 0$). Pour prévenir la violation de la contrainte, les vitesses virtuelles généralisées évoluent dans un des deux demi-espaces de $\hat{\mathbb{U}}$ délimités par l'espace $T_{\varepsilon(t)}(\mathbb{U}^\varphi)$ tangent au sous-espace $\mathbb{U}^\varphi$ en $\boldsymbol{\varepsilon}(t)$. Le vecteur $\partial_\eta \varphi(\boldsymbol{\varepsilon}(t))$ est un vecteur normal à $T_{\varepsilon(t)}(\mathbb{U}^\varphi)$. Ainsi, si l'on considère la décomposition d'une vitesse virtuelle généralisée $\dot{\boldsymbol{\varepsilon}}^*(t)$ en une composante tangentielle $\boldsymbol{\varepsilon}^*(t)''$ et une composante normale $\dot{\boldsymbol{\varepsilon}}^*(t)\perp$ (Figure 9), on a :

$$ \exists \gamma(t) \in \mathbb{R} : \quad \dot{\boldsymbol{\varepsilon}}^*(t)\perp = \gamma(t) \frac{\partial \varphi}{\partial \eta}(\boldsymbol{\varepsilon}(t)) \qquad (1.88) $$
<figure><img src="image_18.png" /><figcaption>Figure 9 Décomposition du vecteur de vitesse généralisée.</figcaption></figure>

Conformément à l'équation (1.86), on a :

$$ \frac{\partial \varphi}{\partial \boldsymbol{\eta}}(\boldsymbol{\varepsilon}(t)) \cdot \dot{\boldsymbol{\varepsilon}}^*(t) = \frac{\partial \varphi}{\partial \boldsymbol{\eta}}(\boldsymbol{\varepsilon}(t)) \cdot (\dot{\boldsymbol{\varepsilon}}^*(t)^\| + \dot{\boldsymbol{\varepsilon}}^*(t)^\perp) = \frac{\partial \varphi}{\partial \boldsymbol{\eta}}(\boldsymbol{\varepsilon}(t)) \cdot \left( \gamma(t) \frac{\partial \varphi}{\partial \boldsymbol{\eta}}(\boldsymbol{\varepsilon}(t)) \right) \geq 0 \quad (1.89) $$

Par conséquent, le scalaire $\gamma(t)$ doit être positif (Figure 10).

D'autre part, puisque les liaisons sont considérées comme parfaites, la puissance virtuelle des efforts issus des contraintes est positive ou nulle.

$$ \sum_{i=1}^{n} Q_{i}^{\varphi}(\boldsymbol{\varepsilon}(t)) \dot{\boldsymbol{\varepsilon}}_{i}^{*}(t) \geq 0 \quad (1.90) $$

Si l'on considère des vitesses virtuelles généralisées appartenant à l'espace tangent $T_{\varepsilon(t)}(\mathbb{U}^\varphi)$ (i.e. telles que $\gamma(t) = 0$) de directions opposées, on arrive à la conclusion suivante :

$$ \sum_{i=1}^{n} Q_{i}^{\varphi}(\boldsymbol{\varepsilon}(t)) \dot{\boldsymbol{\varepsilon}}_{i}^{*}(t)^{\|} = 0 \quad (1.91) $$

Comme pour les contraintes holonomes, on a finalement :

$$ \frac{\partial \varphi}{\partial \boldsymbol{\eta}} \perp \dot{\boldsymbol{\varepsilon}}^{*\|} \quad \text{et} \quad \mathbf{Q}^{\varphi} \perp \dot{\boldsymbol{\varepsilon}}^{*\|} \quad (1.92) \\ \forall \boldsymbol{\varepsilon}^{*\|} \in T_{\varepsilon(t)}(\mathbb{U}^{\varphi}) \quad \text{où} \quad \mathbb{U}^{\varphi} = \{\boldsymbol{\eta} \in \mathbb{U} : \varphi(\boldsymbol{\eta}) = 0\} $$

où l'ensemble $T_{\varepsilon(t)}(\mathbb{U}^\varphi)$ correspond à l'espace tangent au sous-espace $\mathbb{U}^\varphi$ en $\boldsymbol{\varepsilon}(t)$.
<figure><img src="image_19.png" /><figcaption>Figure 10 Interprétation géométrique de l'orientation des vitesses généralisées vis à vis du plan tangent de l'hyper surface issue d'une contrainte holonome.</figcaption></figure>

Puisque l'espace cotangent $\mathrm{coT}_{\varepsilon(t)}(\mathbb{U}^{\varphi})$ est de dimension 1, on déduit que :

$$ \frac{\partial \varphi}{\partial \eta}(\boldsymbol{\varepsilon}(t)), \mathbf{Q}^{\varphi}(\boldsymbol{\varepsilon}(t)) \in \mathrm{coT}_{\varepsilon(t)}(\mathbb{U}^{\varphi}) \Rightarrow \exists \lambda(t) \in \mathbb{R} : \mathbf{Q}^{\varphi}(\boldsymbol{\varepsilon}(t)) = \lambda(t) \frac{\partial \varphi}{\partial \eta}(\boldsymbol{\varepsilon}(t)) \quad (1.93) $$

En exploitant la relation de colinéarité de (1.93), on donne une nouvelle expression à la loi de comportement des liaisons parfaites pour une vitesse virtuelle généralisée $\dot{\boldsymbol{\varepsilon}}^*(t)^{\perp}$ définie en (1.88) :

$$ \mathbf{Q}^{\varphi}(\boldsymbol{\varepsilon}(t)) \cdot \dot{\boldsymbol{\varepsilon}}^*(t)^{\perp} = \left( \lambda(t) \frac{\partial \varphi}{\partial \eta}(\boldsymbol{\varepsilon}(t)) \right) \cdot \left( \gamma(t) \frac{\partial \varphi}{\partial \eta}(\boldsymbol{\varepsilon}(t)) \right) \geq 0 \quad \Rightarrow \quad \lambda(t) \geq 0 \quad (1.94) $$

Imposant au coefficient scalaire $\lambda(t)$ d'être positif pour une contrainte unilatérale réalisée.

### 5.6.4 Multiplicateurs de Lagrange

En conclusion, pour un ensemble de $c_1$ contraintes holonomes, $c_2$ contraintes cinématiques et $c_3$ contraintes unilatérales, la puissance virtuelle des efforts visant à garantir le maintien de ces contraintes est de la forme :

$$ \begin{gathered} P^{\text{liaisons}}(t) = \sum_{i=1}^{n} Q_i^{\text{liaisons}}(\boldsymbol{\varepsilon}(t)) \dot{\varepsilon}_i^*(t) \quad (1.95) \\ \text{avec } Q_i^{\text{liaisons}}(\boldsymbol{\varepsilon}(t)) = \sum_{j=1}^{c_1+c_2+c_3} \lambda_j(t) L_{ji}(\boldsymbol{\varepsilon}(t)) \end{gathered} $$

Chaque $\lambda_j(t)$ est une nouvelle inconnue scalaire associée à la contrainte $j$, communément appelé **multiplicateur de Lagrange**.

$L_{ji}(\boldsymbol{\varepsilon}(t))$ est la puissance virtuelle unitaire (i.e. $\lambda_j(t)=1$) induite par la contrainte $j$ dans un mouvement figé d'entraînement $\boldsymbol{\varepsilon}^i(t)$. Le vecteur $L_j$ peut être perçu comme la force
généralisée issue de la contrainte $j$ pour une réponse unitaire. A noter que le terme $L_{ji}(\varepsilon(t))$, malgré sa notation ambiguë, n'a rien à voir avec le lagrangien du système matériel.

$$L_{ji}(\varepsilon(t)) = \frac{\partial \varphi^j}{\partial \eta_i}(\varepsilon(t)) \quad \text{pour une contrainte holonome (ou unilatérale réalisée)} \quad \varphi^j(\eta).$$

$$L_{ji}(\varepsilon(t)) = a_i^j(\varepsilon(t)) \quad \text{pour une contrainte cinématique} \quad \psi^j(\eta, \hat{\eta}) = \mathbf{a}^j(\eta) \cdot \hat{\eta}.$$

Ainsi, le mouvement $\varepsilon(t)$ d'un système matériel conditionné par $j$ contraintes et soumis à des efforts caractérisés par leur vecteur $\mathbf{Q}$ des forces généralisées est régi par les équations de Lagrange ci-dessous :

$$\begin{array}{l@{\quad}l@{\quad}l}
\displaystyle \frac{d}{dt} \frac{\partial \kappa}{\partial \hat{\eta}_i} (\varepsilon(t), \dot{\varepsilon}(t)) - \frac{\partial \kappa}{\partial \eta_i} (\varepsilon(t), \dot{\varepsilon}(t)) = Q_i(\varepsilon(t), \dot{\varepsilon}(t)) + \sum_{j=1}^{c} \lambda_j(t) L_{ji}(\varepsilon(t)) & 1 \le i \le n & \\
\varphi^j(\varepsilon(t)) = 0 & 1 \le j \le c_1 & (1.96) \\
\mathbf{a}^j(\varepsilon(t)) \cdot \dot{\varepsilon}(t) = 0 & c_1 + 1 \le j \le c_1 + c_2 & \\
\varphi^j(\varepsilon(t)) \ge 0 & c_1 + c_2 + 1 \le j \le c_1 + c_2 + c_3 &
\end{array}$$

# 6 Petit guide à l'attention de l'informaticien

En guise de conclusion, je souhaite d'abord justifier l'intérêt de l'utilisation des équations de Lagrange pour l'animation dynamique. Puis, avec la même volonté de « vulgarisation », je trace les grandes lignes de la modélisation lagrangienne avec une optique d'implantation informatique. Je termine ce chapitre en introduisant une formulation équivalente des équations de Lagrange plus couramment utilisée au sein de la communauté de l'animation.

## 6.1 Quelques bonnes raisons pour privilégier le formalisme lagrangien

### 6.1.1 Généricité

Un des avantages du formalisme lagrangien réside dans sa généricité. En effet, l'équation de Newton (i.e. la loi fondamentale de la dynamique) se limite à décrire le mouvement d'un point matériel. Par extension, l'équation de Newton-Euler conditionne l'évolution des solides rigides seulement. Les équations de Lagrange permettent de décrire le mouvement d'une grande variété de systèmes englobant les points matériels ainsi que les corps rigides. Comme je le détaillerai après, il suffit d'être en mesure d'exprimer la fonction d'énergie cinétique d'un système matériel pour déduire les équations qui régissent son évolution, et ce grâce à la démarche lagrangienne. L'impact informatique d'une telle particularité est loin d'être négligeable. En effet, le choix de ce formalisme pour la réalisation d'un moteur d'animation dynamique permet d'envisager une plate-forme logicielle commune et extrêmement modulaire [Baraff 1996] pour l'animation dynamique d'objets de natures très diverses.
### 6.1.2 Prise en compte naturelle des liaisons

Un autre avantage remarquable du formalisme lagrangien vient de la prise en compte naturelle des liaisons au sein du système. Alors que le formalisme newtonien décrit les mouvements individuels des constituants d'un système en imposant des contraintes *a posteriori* pour représenter les liaisons, le formalisme lagrangien s'appuie sur le jeu de paramètres le plus pertinent (vecteur des coordonnées généralisées) pour décrire l'évolution globale du système. A noter que les paramètres choisis ne correspondent pas forcément à des coordonnées spatiales. Ainsi, les coordonnées généralisées peuvent représenter un angle de rotation ou bien un paramètre identifiant un point d'une courbe paramétrée sur laquelle le corps est contraint à se déplacer. Ici encore l'impact informatique est significatif. En optant pour le formalisme lagrangien, l'évolution du système matériel est caractérisé par le plus petit jeu de variables indépendantes que l'on puisse choisir. Par conséquent, le nombre d'inconnues ainsi que le nombre d'équations à résoudre se trouvent être minimaux permettant de réduire la quantité de mémoire nécessaire mais surtout d'accroître la rapidité de calcul des animations.

## 6.2 Mise en œuvre informatique

Une préoccupation majeure de l'animateur dans le passage de la théorie physique à l'implantation informatique est de cerner les données nécessaires et suffisantes pour décrire un système matériel. Dans le formalisme lagrangien, l'énergie cinétique est un élément central qui joue un rôle majeur dans l'expression des équations de Lagrange. Afin d'être en mesure d'exprimer l'énergie cinétique d'un système donné, on doit disposer d'un minimum d'informations sur ses propriétés structurelles qui sont rassemblées dans ce que je qualifie de **forme cinématique**. Une forme cinématique est composée des éléments suivants :

* La **dimension cinématique** $n$ correspond à la dimension de l'espace des configurations étendues $\mathcal{U}$ (i.e. au nombre de degrés de liberté) du système.

* L'**espace maximal d'identification particulière** $\Omega$ ($\Omega \subset \mathbb{Z}^d$ pour un corps discret et $\Omega \subset \mathbb{R}^d$ pour un corps compact) renseigne sur la répartition interne des particules ainsi que sur la dimension intrinsèque $d$ du système.

* La **fonction de positionnement** $\Phi$ exprime, à partir d'un descripteur $\eta \in \mathcal{U}$ de la configuration du système matériel et d'un identifiant de particule $\omega \in \Omega$, donne la position dans l'espace $\mathbb{R}^3$ de cette particule :

$$ \Phi: \begin{cases} \mathcal{U} \times \Omega \to \mathbb{R}^3 \\ (\eta, \omega) \mapsto \Phi(\eta, \omega) \end{cases} \qquad (1.97) $$

Un implantation informatique des concepts lagrangiens doit et peut se contenter de l'explicitation de la forme cinématique des systèmes qu'elle est amenée à coder. En effet, à partir de la forme cinématique d'un système matériel, je suis en mesure de donner l'expression de son énergie cinétique en s'appuyant sur la définition suivante :

$$ \kappa(\eta, \hat{\eta}) = \frac{1}{2} \hat{\eta}^T \mathbf{M}(\eta) \hat{\eta} \quad \text{avec} \quad \mathbf{M}(\eta) = \int_{\Omega} \rho(\omega, t) \frac{\partial \Phi}{\partial \eta}(\eta, \omega)^T \frac{\partial \Phi}{\partial \eta}(\eta, \omega) d\omega \qquad (1.98) $$
De même, le coordonnées $Q_i(\eta, \hat{\eta})$ vecteur des forces généralisées se déduit naturellement de la forme cinématique et des champs des forces donnés grâce aux relations suivantes :

$$Q_i(\eta, \hat{\eta}) = \int_{\Omega^F} \mathbf{F}(\eta, \hat{\eta}, \omega) \cdot \frac{\partial \Phi}{\partial \eta_i}(\eta, \omega) d\omega \quad \text{pour une densité d'efforts } \mathbf{F}: \mathbb{U} \times \hat{\mathbb{U}} \times \Omega^F \to \mathbb{R}^3.$$

$$Q_i(\eta, \hat{\eta}) = \mathbf{F}(\eta, \hat{\eta}) \cdot \frac{\partial \Phi}{\partial \eta_i}(\eta, \tilde{\omega}) \quad \text{pour un effort } \mathbf{F}(\eta, \hat{\eta}) \text{ appliqué en } \tilde{\omega}.$$

$$Q_i(\eta, \hat{\eta}) = \int_{\Omega^F} \frac{-\partial \tilde{U}}{\partial \eta_i}(\eta, \omega) d\omega \quad \text{pour une énergie potentielle } \tilde{U}: \mathbb{U} \times \Omega^F \to \mathbb{R}.$$

## 6.3 Formulation équivalente des équations de Lagrange

A partir de l'expression particulière de l'énergie cinétique (1.98), on se propose de donner une formulation équivalente des équations de Lagrange ci-dessous

$$\frac{d}{dt} \frac{\partial \kappa}{\partial \hat{\eta}_i} (\varepsilon(t), \dot{\varepsilon}(t)) - \frac{\partial \kappa}{\partial \eta_i} (\varepsilon(t), \dot{\varepsilon}(t)) = Q_i(\varepsilon(t), \dot{\varepsilon}(t)) + \sum_{j=1}^{c} \lambda_j(t) L_{ji}(\varepsilon(t)) \quad 1 \le i \le n \quad (1.99)$$

### 6.3.1 Equations de Lagrange

On s'intéresse d'abord à modifier le membre de gauche de l'équation (1.99). Ainsi, on écrit :

$$\frac{\partial \kappa}{\partial \hat{\eta}_i}(\varepsilon(t), \dot{\varepsilon}(t)) = M_{i\circ}(\varepsilon(t))\dot{\varepsilon}(t) \quad (1.100)$$

Après dérivation par rapport au temps,

$$\frac{d}{dt} \frac{\partial \kappa}{\partial \hat{\eta}_i} (\varepsilon(t), \dot{\varepsilon}(t)) = M_{i\circ}^*(\varepsilon(t)) \ddot{\varepsilon}(t) + \left[ \sum_{j=0}^{n} \dot{\varepsilon}_j(t) \frac{\partial M_{i\circ}}{\partial \hat{\eta}_j}(\varepsilon(t)) \right] \dot{\varepsilon}(t) \quad (1.101)$$

D'autre part,

$$\frac{\partial \kappa}{\partial \eta_i}(\varepsilon(t), \dot{\varepsilon}(t)) = \frac{1}{2} \dot{\varepsilon}(t)^T \frac{\partial M}{\partial \eta_i}(\varepsilon(t)) \dot{\varepsilon}(t) \quad (1.102)$$

Pour alléger les notations, on choisit d'introduire deux nouveaux termes $\mathbf{A}^i(\varepsilon(t))$ et $\mathbf{B}^i(\varepsilon(t))$.

$$\mathbf{B}^i(\varepsilon(t)), \mathbf{A}^i(\varepsilon(t)) \in M_{n+1}(\mathbb{R}) \\ B_{jk}^i(\varepsilon(t)) = \frac{\partial M_{jk}}{\partial \eta_i}(\varepsilon(t)), \quad A_{jk}^i(\varepsilon(t)) = \frac{\partial M_{ik}}{\partial \eta_j}(\varepsilon(t)) = B_{ik}^j(\varepsilon(t)) \quad (1.103)$$

De cette manière, le terme de gauche de l'équation (1.99) devient, en utilisant les matrices définies en (1.103) :

$$\frac{d}{dt} \frac{\partial \kappa}{\partial \hat{\eta}_i} (\varepsilon(t), \dot{\varepsilon}(t)) - \frac{\partial \kappa}{\partial \eta_i} (\varepsilon(t), \dot{\varepsilon}(t)) = M_{i\circ}^{*\omega}(\varepsilon(t)) \ddot{\varepsilon}(t) + \dot{\varepsilon}(t)^T \left[ \mathbf{A}^i(\varepsilon(t)) - \frac{1}{2} \mathbf{B}^i(\varepsilon(t)) \right] \dot{\varepsilon}(t) \quad (1.104)$$
On définit le vecteur $\mathbf{In}(\eta, \hat{\eta})$ comme étant le vecteur des forces généralisées d'inertie,

$$ \mathbf{In}(\eta, \hat{\eta}) \in \mathbb{R}^{n+1} : In_i(\eta, \hat{\eta}) = \hat{\eta}^T \left[ A^i(\eta) - \frac{1}{2} B^i(\eta) \right] \hat{\eta} \quad 1 \le i \le n \quad (1.105) $$

le vecteur $\mathbf{E}(\eta, \hat{\eta})$ correspond au vecteur des forces généralisées totales appliquées au système matériel.

$$ \mathbf{E}(\eta, \hat{\eta}) = \mathbf{Q}(\eta, \hat{\eta}) - \mathbf{In}(\eta, \hat{\eta}) \quad (1.106) $$

Par conséquent, les équations de Lagrange sont équivalentes au système différentiel suivant :

$$ \mathbf{M}^{**}(\varepsilon(t))\ddot{\varepsilon}(t) = \mathbf{E}(\varepsilon(t), \dot{\varepsilon}(t)) + \mathbf{L}^{*}(\varepsilon(t))^{T} \lambda(t) \quad (1.107) $$

### 6.3.2 Equations des contraintes

Mais le système différentiel (1.107) contient plus d'inconnues que d'équations. En effet les multiplicateurs de Lagrange issus des contraintes imposées au système viennent s'ajouter aux inconnues que constituent les coordonnées généralisées. Afin d'injecter de nouvelles équations au système (1.107), on introduit les dérivées des équations de contrainte. Plus précisément, on dérive chaque équation afin de faire apparaître dans leur expression le vecteur des accélérations généralisées $\ddot{\varepsilon}(t)$.

$$ \mathbf{L}^{*}(\varepsilon(t))\ddot{\varepsilon}(t) + \mathbf{D}(\varepsilon(t), \dot{\varepsilon}(t)) = 0 \quad (1.108) $$

La matrice $\mathbf{L}(\varepsilon(t))$ est composée des termes $L_{ji}(\varepsilon(t))$, puissances virtuelles unitaires des efforts de maintien des contraintes. L'expression du terme $\mathbf{D}(\varepsilon(t), \dot{\varepsilon}(t))$ varie selon la façon dont sont gérées les contraintes. Je détaille par la suite plusieurs stratégies de gestion des contraintes ainsi que l'expression du terme $\mathbf{D}(\varepsilon(t), \dot{\varepsilon}(t))$ associé.

En rassemblant les équations (1.107) et (1.108), on débouche sur un système complet d'équations différentielles qui consiste en une forme alternative équivalente des équations de Lagrange et des équations de contraintes.

$$ \begin{bmatrix} \mathbf{M}^{**}(\varepsilon(t)) & -\mathbf{L}^{*}(\varepsilon(t))^{T} \\ -\mathbf{L}^{*}(\varepsilon(t)) & \mathbf{0} \end{bmatrix} \begin{bmatrix} \ddot{\varepsilon}(t) \\ \lambda(t) \end{bmatrix} = \begin{bmatrix} \mathbf{E}(\varepsilon(t), \dot{\varepsilon}(t)) \\ \mathbf{D}(\varepsilon(t), \dot{\varepsilon}(t)) \end{bmatrix} \quad (1.109) $$

Je reviens maintenant à l'expression du terme $\mathbf{D}(\varepsilon(t), \dot{\varepsilon}(t))$ figurant dans les équations de contraintes. Je signale au lecteur qu'une étude détaillée de ce problème a été réalisée par Yannick Remion pour les besoins de son habilitation à diriger des recherches [Remion 2000]. Je me contente donc de donner les grandes lignes du développement renvoyant au travail précédemment cité pour de plus amples détails.
Dans une logique théorique, le formalisme lagrangien garantit que les contraintes sont réalisées en tout instant. Par conséquent les équations de contraintes et leurs dérivées successives sont identiquement nulles.

$$ \frac{d^2\varphi^i}{dt^2}(\boldsymbol{\epsilon}(t)) = 0 \quad \text{et} \quad \frac{d}{dt}(\alpha^i(\boldsymbol{\epsilon}(t))\cdot\dot{\boldsymbol{\epsilon}}(t)) = 0 \qquad (1.110) $$

Dans ce cas, le terme $D(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t))$ équivaut à :

$$ D(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) = \dot{\boldsymbol{\epsilon}}(t)^T \frac{\partial L}{\partial \boldsymbol{\eta}} (\boldsymbol{\epsilon}(t))^T \dot{\boldsymbol{\epsilon}}(t) \qquad (1.111) $$

Mais le raisonnement précédent n'est en fait valide que dans une optique de résolution formelle des équations de Lagrange. Lors d'une résolution numérique avec un échantillonnage de l'échelle du temps, les paramètres cinématiques ont des évolutions qui approximent leurs évolutions réelles, induisant une dérive dans la résolution des contraintes. Baumgarte [Baumgarte 1972] propose une solution pour contrecarrer ce problème de dérive en ajoutant des termes correcteurs pondérées par des scalaires $\gamma_1, \gamma_2, \gamma_3$ (1.112).

$$ \begin{gathered} \frac{d^2\varphi^i}{dt^2}(\boldsymbol{\epsilon}(t)) + \gamma_1 \frac{d\varphi^i}{dt}(\boldsymbol{\epsilon}(t)) + \gamma_2\varphi^i(\boldsymbol{\epsilon}(t)) = 0 \qquad (1.112) \\ \frac{d}{dt}(\alpha^i(\boldsymbol{\epsilon}(t))\cdot\dot{\boldsymbol{\epsilon}}(t)) + \gamma_3\alpha^i(\boldsymbol{\epsilon}(t))\cdot\dot{\boldsymbol{\epsilon}}(t) = 0 \end{gathered} $$

Le terme $D_i(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t))$ associé à la contrainte $i$ devient alors, pour une contrainte holonome et une contrainte cinématique :

$$ \begin{gathered} D_i(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) = \dot{\boldsymbol{\epsilon}}(t)^T \frac{\partial^2\varphi^i}{\partial\boldsymbol{\eta}^2} (\boldsymbol{\epsilon}(t))^T \dot{\boldsymbol{\epsilon}}(t) + \gamma_1 \frac{\partial\varphi^i}{\partial\boldsymbol{\eta}} (\boldsymbol{\epsilon}(t)) \dot{\boldsymbol{\epsilon}}(t) + \gamma_2\varphi^i (\boldsymbol{\epsilon}(t)) \qquad (1.113) \\ D_i(\boldsymbol{\epsilon}(t), \dot{\boldsymbol{\epsilon}}(t)) = \dot{\boldsymbol{\epsilon}}(t)^T \frac{\partial\alpha^i}{\partial\boldsymbol{\eta}} (\boldsymbol{\epsilon}(t))^T \dot{\boldsymbol{\epsilon}}(t) + \gamma_3\alpha^i (\boldsymbol{\epsilon}(t)) \cdot \dot{\boldsymbol{\epsilon}}(t) \end{gathered} $$
## 6.4 Architecture logicielle

En résumé, l'implantation d'un moteur d'animation dynamique lagrangien doit s'articuler autour des grandes étapes suivantes :

### Choix d'une forme cinématique

① Explicitation de la matrice des masses généralisées $M(\eta)$ et du vecteur des forces généralisées $E(\eta, \hat{\eta})$

### Choix d'une modélisation des contraintes

② Explicitation de la force généralisée $L(\eta)$ des contraintes et du terme correcteur $D(\eta, \hat{\eta})$

### Boucle sur l'intervalle de temps

③ Evaluation des termes $M(\varepsilon(t))$, $E(\varepsilon(t), \dot{\varepsilon}(t))$, $L(\varepsilon(t))$, $D(\varepsilon(t), \dot{\varepsilon}(t))$ de l'équation (1.109) à l'instant $t$.

④ Résolution numérique du système linéaire (1.109) pour déterminer les valeurs des vecteurs $\ddot{\varepsilon}(t)$ et $\dot{\lambda}(t)$ à l'instant $t$.

⑤ Résolution numérique itérative de l'équation différentielle du second ordre déduite de l'étape ④ pour évaluer les vecteurs $\varepsilon(t+h)$ et $\dot{\varepsilon}(t+h)$.

⑥ $t \leftarrow t+h$, retour à l'étape ③.

L'étape ⑤ exige d'être capable de trouver une solution approchée d'une équation différentielle du second ordre. Le chapitre suivant, consacré aux méthodes d'intégration numérique, décrit dans le détail la mise en œuvre d'algorithmes de résolution itérative de ces équations différentielles.
<figure><img src="image_20.png" /></figure>
## 1 Introduction

La particularité de l'animation dynamique, en comparaison avec l'animation par ordinateur qualifiée d'animation traditionnelle, réside dans le fait que les évolutions des objets composant la scène sont déduites de la résolution des équations du mouvement.

Les équations de Lagrange, qui correspondent aux équations du mouvement pour lesquelles nous avons opté, sont équivalentes après simplifications formelles à un système d'équations différentielles (cf Chapitre 1 : Justification des équations de Lagrange). Dans une optique de conception informatique d'un moteur lagrangien généraliste, les équations différentielles déduites d'une modélisation mécanique donnée peuvent être tout à fait quelconques. Ainsi, une résolution exacte du système différentiel par le biais d'une méthode analytique n'est pas envisageable dans notre cadre de travail.

Dans ces conditions, l'intégration numérique semble être la seule issue : elle permet d'obtenir une approximation de la solution exacte du problème différentiel. Les méthodes d'intégration numérique sont nombreuses et variées, adaptées à des problèmes de natures diverses. Mais ces méthodes reposent toutes sur le même schéma : l'échelle du temps, paramètre des évolutions, est discrétisée afin de calculer de façon itérative les valeurs successives d'une solution approchée.

Ce chapitre a pour vocation, après avoir introduit la terminologie liée aux méthodes d'intégration numérique, d'exposer des critères objectifs de comparaison. Par la suite, je donne une classification, selon leur champ d'application, des méthodes les plus courantes. En effet, il semble illusoire de vouloir dresser une liste exhaustive des méthodes d'intégration, un domaine vaste qui constitue encore à ce jour un thème de recherche actif.
# 2 Définitions

## 2.1 Discrétisation du problème

On cherche à déterminer une solution approchée de l'équation différentielle (2.1) sur un intervalle de temps borné $[t_0, t_0 + T]$ avec une condition sur la valeur de la solution en $t_0$.

$$ egin{cases} \frac{dy}{dt}(t) = f(t, y(t)) \\ y(t_0) = y_0 \end{cases} \quad \text{où} \quad f: [t_0, t_0 + T] \times \mathbb{R}^n \to \mathbb{R}^n \qquad (2.1) $$

Une équation différentielle telle que (2.1) munie d'un condition initiale est appelée **problème de Cauchy**. Si la fonction $f$ est $\Lambda$ lipschitzienne par rapport à $y$ sur l'ensemble $[t_0, t_0 + T] \times \mathbb{R}^n$ (propriété synthétisée par l'équation (2.2)), le problème de Cauchy admet alors une solution unique $y(t)$ définie sur l'intervalle $[t_0, t_0 + T]$.

$$ \forall t \in [t_0, t_0 + T], \forall y_1, y_2 \in \mathbb{R}^n : \|f(t, y_1) - f(t, y_2)\| \le \Lambda \|y_1 - y_2\| \qquad (2.2) $$

Soit $(t_n)_{0 \le n \le N}$ une subdivision croissante de l'intervalle $[t_0, t_0 + T]$.

$$ t_0 < t_1 < \dots < t_{N-1} < t_N = t_0 + T $$

On pose alors,

$$ \begin{cases} h_n = t_{n+1} - t_n \\ h_{\max} = \max_{0 \le n < N} (h_n) \end{cases} \qquad 0 \le n < N \qquad (2.3) $$

On note $y_n$ la valeur approchée calculée de la solution exacte $y(t_n)$ en $t_n$. Par suite, on note $f_n$ la valeur de $f(t_n, y_n)$. Ainsi,

$$ y_n \approx y(t_n) \quad \text{et} \quad f_n = f(t_n, y_n) \qquad (2.4) $$

## 2.2 Méthodes d'intégration

Une méthode d'intégration est un procédé permettant de déterminer de façon itérative les valeurs approchées successives de la solution du problème (2.1).

* Une méthode d'intégration est dite **explicite** si le calcul de $y_{n+1}$ ne dépend que de $y_n, y_{n-1}, \dots, y_{n-r}$, les valeurs de la solution approchée calculées antérieurement.

* Une méthode d'intégration est dite **implicite** si le calcul de $y_{n+1}$ dépend de $y_n, y_n, y_{n-1}, \dots, y_{n-r}$, mais aussi de $f_{n+1}$, la dérivée de la solution en $t_{n+1}$ approchée par la valeur $y_{n+1}$.
* Une **méthode d'intégration à un pas** est une méthode qui ne s'appuie que sur la connaissance de $y_n$ (éventuellement de $f_{n+1}$ si la méthode est implicite) pour le calcul de $y_{n+1}$.

* Une **méthode d'intégration à r+1 pas** (i.e. à pas multiples) est une méthode qui s'appuie sur la connaissance de $y_n, y_{n-1}, \dots, y_{n-r}$ (ou $f_{n+1}, f_n, y_{n-1}, \dots, y_{n-r-1}$ si la méthode est implicite) pour le calcul de $y_{n+1}$.

## 2.3 Consistance, stabilité et convergence des méthodes d'intégration

J'expose dans ce paragraphe des propriétés associées aux méthodes d'intégration, permettant de les caractériser et de les comparer.

### 2.3.1 Consistance

L'erreur de consistance $e_n$ à l'étape $n$ correspond à la différence de la solution exacte $y(t_{n+1})$ et de sa valeur approchée $y_{n+1}$ obtenue par intégration numérique en supposant que $y(t_n) = y_n$.

$$e_n = y(t_{n+1}) - y_{n+1} \qquad (2.5)$$

<figure><img src="image_21.png" /><figcaption>Figure 11 Représentation géométrique de l'erreur de consistance.</figcaption></figure>

L'**erreur de consistance** $e_n$ correspond à une mesure de l'erreur locale générée par la méthode d'intégration lors d'une itération. Cette erreur est aussi qualifiée d'erreur de troncature dans certains ouvrages [Cartwright 1992] en référence à la troncature du développement de Taylor de la fonction $f(t, y)$ pour la construction de certaines méthodes d'intégration numérique.

Une méthode d'intégration est **consistante** si :

$$\sum_{n=0}^{N} |e_n| \xrightarrow{h_{\max} \to 0} 0 \qquad (2.6)$$

Une méthode d'intégration est d'**ordre** $p$ si :

$$\begin{gathered} \forall n, 0 \le n \le N : e_n = o(h_{\max}^{p+1}) \quad h_{\max} \to 0 \qquad (2.7) \\ (\text{i.e. } \exists C > 0 : |e_n| \le C.h_{\max}^{p+1}) \end{gathered}$$

L'ordre d'une méthode d'intégration renseigne sur la précision de la méthode d'intégration. En effet, l'ordre conditionne la rapidité de convergence du terme (2.6) qui représente la somme des
erreurs locales lors du calcul de la solution approchée. Il est clair que plus l'ordre d'une méthode d'intégration est élevé, plus cette méthode est consistante.

### 2.3.2 Stabilité

La notion de stabilité est quelque peu ambiguë dans ce contexte car on peut à la fois parler de la stabilité de la solution exacte de l'équation différentielle (inhérente à la nature du problème de Cauchy) et de la stabilité de la méthode d'intégration (liée à l'algorithme de la méthode).

Une solution $y(t)$ d'un problème de Cauchy (2.1) est stable (au sens de Liapunov [Dieudonné 1968]) si, pour une valeur initiale $\tilde{y}_0$ proche de $y_0$, la solution exacte $\tilde{y}(t)$ correspondante demeure voisine de $y(t)$ au fil du temps :

$$ 
\begin{array}{l}
\forall \varepsilon > 0, \exists \delta > 0 : \forall y_0, \tilde{y}_0 \in \mathbb{R}^n \text{ tels que } \|y_0 - \tilde{y}_0\| \le \delta, \text{ on a:} \\
\qquad \|y(t) - \tilde{y}(t)\| \le \varepsilon \quad \forall t \ge t_0
\end{array}
\qquad (2.8)
$$

Une méthode d'intégration est stable si :

$$ 
\exists S > 0 : \forall (y_n)_{0 \le n \le N}, (\tilde{y}_n)_{0 \le n \le N} \text{ deux suites de valeurs approchées} \\
\max_{0 \le n \le N} \|\tilde{y}_n - y_n\| \le S (\|\tilde{y}_0 - y_0\|) \qquad (2.9)
$$

On peut illustrer la notion de stabilité de la méthode d'intégration ainsi : pour une erreur initiale $\|\tilde{y}_0 - y_0\|$, le calcul itératif des valeurs successives $y_n$ et $\tilde{y}_n$ génère un écart $\max\|\tilde{y}_n - y_n\|$ majorable, donc contrôlable. L'équation (2.9) montre que ce critère demeure fortement dépendant de la stabilité de l'équation différentielle que la méthode d'intégration cherche à résoudre. Il semble en effet difficile de formuler un critère objectif de stabilité, complètement indépendant du problème de Cauchy considéré.

Malheureusement dans la pratique, le critère (2.9) est assez difficile à mettre en oeuvre. On trouve dans des ouvrages d'analyse numérique [Hairer 1991] [Cartwright 1992] des définitions alternatives de la stabilité.

La **stabilité absolue** (ou $\lambda$ h-stabilité) permet de quantifier la stabilité de la méthode d'intégration pour un problème modèle correspondant à l'équation différentielle linéaire suivante (Figure 12) :

$$ 
\left\{ 
\begin{array}{l}
\displaystyle \frac{dy}{dt}(t) = \lambda y(t) \\
\\
\displaystyle y(t_0) = \alpha
\end{array}
\right.
\quad \text{pour } \alpha, \lambda \in \mathbb{C} \qquad (2.11)
$$

Le problème (2.11) peut être résolu de façon analytique. L'expression de la solution exacte est donnée par l'équation (2.12).

$$ 
\begin{array}{l}
y(t) = \alpha e^{\lambda(t-t_0)} \\
\\
\text{avec } y(t) \underset{t \to +\infty}{\longrightarrow} 0 \quad \text{pour } \text{Re}(\lambda) < 0
\end{array}
\qquad (2.12)
$$
<figure><img src="image_22.png" /><figcaption>Figure 12 solutions de l'équation différentielle linéaire $y' = λy$, stable pour $λ < 0$ et instable pour $λ > 0$.</figcaption></figure>

Le **domaine de λ h-stabilité** est la région $\{\lambda h : \lambda \in \mathbb{C}, h \in \mathbb{R}^+\}$ du plan complexe dans laquelle la méthode d'intégration donne, pour une constante $\lambda$ et un pas de temps $h$, une solution approchée du problème (2.11) qui converge aussi vers 0 quand $t$ tend vers $+\infty$. De par leur définition, les domaines de $\lambda h$-stabilité sont tous homothétiques à un rapport $1/h$ pour un pas de temps $h$ donné (Figure 13).

<figure><img src="image_23.png" /><figcaption>Figure 13 Domaines homothétiques de $\lambda h$-stabilité de la méthode d'Euler explicite pour $h=1$, $h=2/3$ et $h=1/2$.</figcaption></figure>

Ainsi, pour un pas de temps $h$ de plus en plus petit, les valeurs « acceptables » de $\lambda$, pour que la solution approchée converge, sont de plus en plus grandes. Quand le pas de temps $h$ tend vers zéro, les domaines de $\lambda h$-stabilité convergent plus ou moins rapidement vers le demi plan complexe $\{z \in \mathbb{C} : \mathrm{Re}(z) < 0\}$.

Naturellement, une méthode d'intégration est d'autant plus stable (au sens de la stabilité absolue) que son domaine de $\lambda h$-stabilité est grand.

Une méthode d'intégration est **$\lambda h$-stable** si son domaine de $\lambda h$-stabilité est le demi plan complexe $\{z \in \mathbb{C} : \mathrm{Re}(z) < 0\}$.

Les domaines de $\lambda h$-stabilité, même s'ils constituent un critère de comparaison des méthodes d'intégration, ne sont vraiment significatifs que dans le cadre d'un problème linéaire, c'est à
dire des problèmes où $f(t, y) = \lambda y$. Des approches moins répandues donnent une définition de la **stabilité absolue non-linéaire** à partir du problème modèle suivant :

$$ egin{cases} \frac{dy}{dt}(t) = \lambda y(t)(1 - y(t)) \\ y(a) = \alpha \end{cases} \quad \text{pour } \alpha, \lambda \in \mathbb{C} \qquad (2.13) $$

Bien que non-linéaire, ce critère s'applique à une équation différentielle particulière et ne permet pas d'extrapoler de façon certaine sur la stabilité de la méthode pour une équation quelconque. Finalement, de par la difficulté de quantifier la notion de stabilité, ces démarches doivent se cantonner à donner une simple indication de stabilité pour une méthode donnée.

### 2.3.3 Convergence

Comme nous l'avons vu dans le paragraphe 2.3.1, la consistance de la méthode mesure la précision de la méthode, en supposant qu'à chaque itération la valeur initiale coïncide avec la solution exacte. En pratique, ce cas de figure n'est pas envisageable. Au mieux, on peut espérer qu'au fur et à mesure de la progression du calcul, la valeur approchée reste voisine de la solution exacte. Cette idée est mathématiquement formalisée par la stabilité de la méthode d'intégration. Ainsi, la qualité de la solution approchée dépend simultanément de la consistance et de la stabilité de la méthode d'intégration. Cette notion est synthétisée par la **convergence** de la méthode [Demailly 1991].

On commence par donner la définition de l'**erreur globale** :

$$ \max_{0 \le n \le N} \| y(t_n) - y_n \| \qquad (2.14) $$

L'erreur globale correspond à l'écart entre la solution exacte et la valeur issue des itérations successives lors de l'intégration, englobant l'accumulation des erreurs locales et l'erreur sur la valeur initiale. Une méthode d'intégration est convergente si :

$$ \max_{0 \le n \le N} \| y(t_n) - y_n \| \xrightarrow{h_{\max} \to 0, y_0 \to y(t_0)} 0 \qquad (2.15) $$

#### Théorème de la convergence

Si une méthode d'intégration est à la fois consistante (2.6) et stable (2.9), alors cette méthode est convergente.

A noter que l'erreur globale est une erreur théorique, c'est à dire qu'elle ne prend pas en compte les erreurs d'arrondis qui se produisent lors des calculs successifs. En pratique, lors des évaluations intermédiaires de la fonction $f$ requises par l'ordre de la méthode d'intégration, l'accumulation des erreurs d'arrondis peut nuire à la précision numérique de la méthode.

## 2.4 Classification des problèmes d'intégration numérique

A l'aide des définitions précédentes, je suis en mesure de classer les problèmes d'intégration numérique en sous catégories. J'en cite ici deux d'importance majeure.
Les **problèmes bien posés** sont des problèmes dont l'équation différentielle associée est stable au sens de la définition (2.8). En d'autres termes, la fonction $f(t, y)$ est suffisamment lisse (i.e. plusieurs fois différentiable par rapport à la variable $y$) pour que la solution ne soit pas trop perturbée par une petite variation de la donnée initiale. Ce genre de problèmes ne requiert pas impérativement de méthodes d'intégration stables.

Les **problèmes raides** sont des problèmes qui imposent d'utiliser un pas de temps beaucoup plus petit que celui requis pour satisfaire la consistance [Lambert 1991]. Par exemple, un problème est qualifié de raide lorsque le sens de variation de la fonction $f$ est beaucoup plus important dans une direction que dans d'autres. Plus précisément, si l'on considère les valeurs propres $(\lambda_i)_{1\le i\le n}$ de la jacobienne de la fonction $f$, le problème est raide si le **coefficient de raideur** $r$ défini en (2.16) est grand.

$$r = \frac{\mathrm{Re}\lambda_n}{\mathrm{Re}\lambda_1} \quad \text{avec} \quad \mathrm{Re}\lambda_n < \dots < \mathrm{Re}\lambda_1 < 0 \qquad (2.16)$$

Les problèmes raides nécessitent des méthodes d'intégration stables pour garantir la convergence de la solution approchée.

# 3 Méthodes d'intégration à un pas

Ce chapitre se focalise sur les méthodes à un pas. Je commence par en donner la forme générale et énoncer les critères de convergence. Puis, j'énumère dans l'ordre chronologique d'apparition quelques méthodes classiques pour finalement aboutir sur l'algorithme de construction de méthodes d'intégration à un pas proposé par Runge et Kutta [Runge 1895] [Kutta 1901] [Crouzeix 1989] [Demailly 1991].

## 3.1 Méthodes explicites

Toute méthode d'intégration explicite à un pas s'écrit sous la forme suivante :

$$y_{n+1} = y_n + h_n \varphi(t_n, y_n, h_n) \quad \text{où} \quad \varphi: [t_0, t_0+T] \times \mathbb{R}^n \times \mathbb{R} \to \mathbb{R}^n \qquad (2.17)$$

### 3.1.1 Critères de convergence

J'expose dans cette section deux résultats classiques en analyse numérique permettant de quantifier à la fois la consistance et la stabilité des méthodes d'intégration à un pas [Crouzeix 1989] [Demailly 1991].

Le premier critère est une condition nécessaire et suffisante de consistance, aisément démontrable à partir du développement de Taylor de la fonction $f(t, y)$ :

Si $\varphi \in C^p(\mathbb{R} \times \mathbb{R}^n \times \mathbb{R}, \mathbb{R}^n)$ (i.e. $p$ fois continûment dérivable) et,

$$\begin{gathered} \forall (t, y) \in [t_0, t_0+T] \times \mathbb{R}^n: \\ \frac{\partial^i \varphi}{\partial h^i}(t, y, 0) = \frac{1}{i+1} f^{(i)}(t, y) \quad (0 \le i \le p) \quad \Leftrightarrow \quad \begin{array}{l} \text{la méthode est} \\ \text{d'ordre } p+1 \end{array} \end{gathered} \qquad (2.18)$$
Cas particulier ($p = 0$) :

$$ 
\begin{array}{l}
\forall (t, y) \in [t_0, t_0 + T] \times \mathbb{R}^n : \\
\qquad \varphi(t, y, 0) = f(t, y) \Leftrightarrow \text{la méthode est consistante}
\end{array}
\qquad (2.19)
$$

Le deuxième critère est une condition suffisante de stabilité :

$$ 
\begin{array}{l}
\forall (t, y_1, h), (t, y_2, h) \in [t_0, t_0 + T] \times \mathbb{R}^n \times \mathbb{R} : \\
\qquad \|\varphi(t, y_1, h) - \varphi(t, y_2, h)\| \le \Lambda \|y_1 - y_2\| \\
\text{(i.e.} \varphi \text{ est } \Lambda\text{-lipschitzienne par rapport à y)
\end{array}
\qquad (2.20)
$$

### 3.1.2 La méthode d'Euler (1768)

Premièrement, on rappelle que :

$$ \frac{dy}{dt}(t) = \lim_{h \to 0} \frac{y(t+h) - y(t)}{h} \qquad (2.21) $$

La méthode repose sur l'idée que pour un réel $h$ suffisamment petit, le terme de droite de (2.21) constitue une approximation satisfaisante de la dérivée (Figure 14). Par substitution dans (2.1), on obtient le système suivant :

$$ \frac{y(t+h) - y(t)}{h} \approx f(y(t), t) \Leftrightarrow y(t+h) \approx y(t) + h f(y(t), t) \qquad (2.22) $$

L'équation (2.22) amena Euler à formuler la méthode d'intégration suivante :

$$ 
\begin{cases} 
t_{n+1} = t_n + h_n \\ 
y_{n+1} = y_n + h_n f(t_n, y_n) 
\end{cases} 
\qquad (2.23)
$$

<figure><img src="image_24.png" /><figcaption>Figure 14 Représentation géométrique de la méthode d'Euler.</figcaption></figure>

La méthode d'Euler vérifie le critère de consistance (2.19) car $\varphi(t, y, h) = f(t, y)$. Très intuitive mais quelque peu grossière, elle est seulement d'ordre 1 puisque $\partial_t \varphi(t, y, h) = 0$ (cf. critère de consistance (2.18)). L'erreur de consistance est donc mal contrôlée. Cette méthode est
une méthode d'intégration faiblement convergente, même pour une fonction $f(t, y)$ satisfaisant les conditions de stabilité.

### 3.1.3 La méthode de Runge (1895)

Géométriquement, la méthode d'Euler s'appuie sur la pente de la solution de (2.1) en $t_n$, estimée grâce à la valeur $f(t_n, y_n)$, pour déterminer la valeur approchée $y_{n+1}$. L'idée de Runge a été d'utiliser la pente estimée de la solution en $t_{n+1/2} = t_n + h_n/2$, au milieu du sous intervalle $[t_n, t_{n+1}]$, afin de déterminer $y_{n+1}$ de façon plus précise (Figure 15). C'est pourquoi cette méthode d'intégration porte aussi le nom de **méthode du point milieu**.

$$ egin{cases} t_{n+1} = t_n + h_n \\ y_{n+1/2} = y_n + \frac{h_n}{2} f(t_n, y_n) \\ f_{n+1/2} = f\left(t_n + \frac{h_n}{2}, y_{n+1/2}\right) \\ y_{n+1} = y_n + h_n f_{n+1/2} \end{cases} \quad (2.24) $$

<figure><img src="image_25.png" /><figcaption>Figure 15 Représentation géométrique de la méthode de Runge.</figcaption></figure>

La méthode de Runge vérifie aussi le critère de consistance (2.19). Conformément à mon intuition, cette méthode est plus précise que la méthode d'Euler puisque son ordre est égal à 2. En s'appuyant sur l'algorithme (2.24), la fonction $\varphi: [t_0, t_0+T] \times \mathbb{R}^n \times \mathbb{R} \to \mathbb{R}^n$ associée à la méthode du point milieu s'exprime ainsi :

$$ \varphi(t, y, h) = f\left(t + h/2, y + h/2f\left(t, y\right)\right) \quad (2.25) $$

On justifie alors l'ordre de la méthode en vérifiant les égalités du critère de consistance (2.18),

$$ egin{align*} \varphi(t, y, 0) &= f(t, y), \\ \frac{\partial \varphi}{\partial h}(t, y, h) &= \frac{1}{2} \frac{\partial f}{\partial t}(t + h/2, y + h/2f(t, y)) + \frac{1}{2} f(t, y) \frac{\partial f}{\partial y}(t + h/2, y + h/2f(t, y)) \\ &= \frac{1}{2} \left( \frac{\partial f}{\partial t}(t + h/2, y + h/2f(t, y)) \frac{dt}{dt} + \frac{dy(t)}{dt} \frac{\partial f}{\partial y}(t + h/2, y + h/2f(t, y)) \right) \\ \frac{\partial \varphi}{\partial h}(t, y, 0) &= \frac{1}{2} \frac{df}{dt}(t, y) \end{align*} $$
$$
\begin{align*}
\frac{\partial^2 \varphi}{\partial h^2}(t, y, h) &= \frac{\partial}{\partial h} \frac{\partial \varphi}{\partial h}(t, y, h) = \frac{\partial}{\partial h} \left( \frac{1}{2} \frac{df}{dt} \left( t + h/2, y + h/2f(t,y) \right) \right) \\
&= \frac{1}{2} \frac{d}{dt} \frac{\partial f}{\partial h} \left( t + h/2, y + h/2f(t,y) \right) = \frac{1}{2} \frac{d}{dt} \left( \frac{1}{2} \frac{df}{dt} \left( t + h/2, y + h/2f(t,y) \right) \right) \\
\frac{\partial^2 \varphi}{\partial h^2}(t, y, h) &= \frac{1}{4} \frac{d^2 f}{dt^2}(t, y) \neq \frac{1}{3} \frac{d^2 f}{dt^2}(t, y)
\end{align*}
$$

3.1.4 Construction des méthodes explicites à un pas (1901)

C'est avec le concours de Kutta, que Runge mit au point un procédé de construction de méthodes d'intégration à un pas, englobant et unifiant l'élaboration d'un bon nombre de méthodes existantes. Ils partirent de la constatation suivante :

$$
y(t) = y(t_0) + \int_{t_0}^{t} \frac{dy}{dt}(u)du \quad (2.26)
$$

Comme $y(t)$ est solution du problème (2.1), on obtient par substitution :

$$
y(t) = y(t_0) + \int_{t_0}^{t} f(u, y(u)) du \quad (2.27)
$$

Ainsi, sur le sous intervalle $[t_n, t_{n+1}]$, on peut écrire :

$$
y(t_{n+1}) = y(t_n) + \int_{t_n}^{t_{n+1}} f(u, y(u)) du \quad (2.28)
$$

Avec le changement de variable $u = t_n + vh_n$, on obtient :

$$
y(t_{n+1}) = y(t_n) + h_n \int_0^1 f(t_n + vh_n, y(t_n + vh_n)) dv \quad (2.29)
$$

En s'appuyant sur (2.29), on débouche alors sur une écriture générique d'une classe de
méthodes d'intégration à un pas :

$$
\begin{align}
y_{n+1} &= y_n + h_n \varphi(t_n, y_n, h_n) \notag \\
\text{où } \varphi(t_n, y_n, h_n) &\approx \int_0^1 f(t_n + v h_n, y(t_n + v h_n)) dv \tag{2.30}
\end{align}
$$

En définitive, les différentes méthodes d'intégration à un pas de cette nature se distinguent selon
la méthode d'approximation du terme intégral utilisée, d'où le choix du terme d'intégration
numérique pour désigner les méthodes de résolution numérique des équations différentielles. On
présente maintenant en détail l'algorithme de construction de Runge-Kutta.
Soit $(c_n)_{1\le n\le s}$ une subdivision croissante de l'intervalle d'intégration telle que :

$$0 \le c_1 \le c_2 \le \dots \le c_{s-1} \le c_s \le 1 \qquad (2.31)$$

On pose :

$$t_{n,i} = t_n + c_i h_n$$

$$y(t_{n,i}) = y(t_n) + h_n \int_0^{c_i} f(t_n + v h_n, y(t_n + v h_n)) dv \quad 1 \le i \le s \qquad (2.32)$$

On choisit $s$ méthodes de calcul d'intégrale de coefficients respectifs $a_{ij}$ afin d'estimer les $s$ termes intégraux de (2.32).

$$y_{n,i} = y_n + h_n \sum_{j=0}^{i} a_{ij} f_{n,j} \quad 1 \le i \le s \qquad (2.33)$$

$$f_{n,i} = f(t_{n,i}, y_{n,i})$$

Le calcul de $y_{n,i}$ ne dépend que des valeurs de $f_{n,j}$ pour $0 \le j < i$ de sorte que les valeurs successives de $f_{n,i}$ peuvent être obtenues de façon itérative sans avoir à recourir à la résolution d'un système linéaire.

<figure><img src="image_26.png" /><figcaption>Figure 16 Approximation polygonale du terme intégral de l'équation (2.28) pour la méthode d'Euler (a), la méthode de Runge (b) et la méthode de Runge-Kutta 4 (c).</figcaption></figure>

Finalement, l'expression de (2.29) peut être approchée de la manière suivante en choisissant une méthode d'intégration numérique de coefficients $b_i$, approximation polygonale de l'aire sous la courbe $f(t, y)$ estimée à partir des valeurs $f_{n,i}$ (Figure 16).

$$y_{n+1} = y_n + h_n \sum_{i=1}^{s} b_i f_{n,i} \qquad (2.34)$$
En résumé,

$$
\left[
\begin{array}{l}
\left[ \begin{array}{l} t_{n,i} = t_n + c_i h_n \\ y_{n,i} = y_n + h_n \sum_{1 \le j < i} a_{ij} f_{n,j} \\ f_{n,i} = f(t_{n,i}, y_{n,i}) \end{array} \right] \quad 1 \le i \le s \\[1em]
\left\{
\begin{array}{l}
t_{n+1} = t_n + h_n \\[1em]
y_{n+1} = y_n + h_n \sum_{i=1}^{s} b_i f_{n,i}
\end{array}
\right.
\end{array}
\right.
\qquad
\begin{array}{c|ccccccc}
c_1 & & & & & & & \\
c_2 & a_{2,1} & & & & & & \\
c_3 & a_{3,1} & a_{3,2} & & & & & \\
\dots & \dots & \dots & \dots & & & & \\
c_s & a_{s,1} & a_{s,2} & \dots & \dots & a_{s,s-1} & \\
\hline
& b_1 & b_2 & \dots & \dots & b_{s-1} & b_s
\end{array}
$$

Figure 17 Algorithme de la méthode de Runge Kutta explicite à s étages et son tableau de Butcher associé.

Une telle méthode d'intégration est appelée méthode explicite de Runge-Kutta à s étages (i.e. méthode ERK à s étages). Le tableau de Butcher [Butcher 1964] ci-dessus est une représentation symbolique de la méthode, regroupant tous les paramètres intervenant dans sa construction.

* Les termes $(c_i)_{1 \le i \le s}$ sont appelés les nœuds.

* Les termes $(a_{i,j})_{1 \le j \le i \le s}$ sont appelés les poids intérieurs.

* Les termes $(b_i)_{1 \le i \le s}$ sont appelés les poids extérieurs.

### 3.1.5 Classification des méthodes explicites à un pas

Toutes les méthodes d'intégration présentées auparavant peuvent être représentées à l'aide d'un tableau de Butcher. Dans cette section, on donne une liste non exhaustive de méthodes couramment utilisées, classées par ordre de consistance croissant.

#### 3.1.5.1 Méthode d'ordre 1

$$ \begin{array}{c|c} 0 & 0 \\ & 1 \end{array} $$

<caption>Tableau 1 Méthode d'Euler.</caption>

#### 3.1.5.2 Méthodes d'ordre 2

$$ \begin{array}{c|c} 0 & 1/2 \\ & 1/2 \\ \hline & 0 \quad 1 \end{array} $$

(a)

$$ \begin{array}{c|c} 0 & a \\ & a \\ \hline & 1 - \frac{1}{2a} \quad \frac{1}{2a} \end{array} $$

(b)

<caption>Tableau 2 (a) Méthode de Runge (i.e. point milieu)</caption>

(b) Méthode de Runge Kutta 2 avec $0 < a < 1$.
### 3.1.5.3 Méthodes d'ordre 3

<table><caption>Tableau 3 Méthode de Heun.</caption><tr><td>0</td><td></td><td></td></tr><tr><td>1/3</td><td>1/3</td><td></td></tr><tr><td>2/3</td><td>0</td><td>2/3</td></tr><tr><td></td><td>1/4</td><td>0</td><td>3/4</td></tr></table>

### 3.1.5.4 Méthodes d'ordre 4

<table><tr><td>0</td><td></td><td></td></tr><tr><td>1/2</td><td>1/2</td><td></td></tr><tr><td>1/2</td><td>0</td><td>1/2</td></tr><tr><td>1</td><td>0</td><td>0</td><td>1</td></tr><tr><td></td><td>1/6</td><td>2/6</td><td>2/6</td><td>1/6</td></tr></table>

(a)

<table><tr><td>0</td><td></td><td></td></tr><tr><td>1/3</td><td>1/3</td><td></td></tr><tr><td>2/3</td><td>-1/3</td><td>1</td></tr><tr><td>1</td><td>1</td><td>-1</td><td>1</td></tr><tr><td></td><td>1/8</td><td>3/8</td><td>3/8</td><td>1/8</td></tr></table>

(b)

<caption>Tableau 4 (a) Méthode de Runge Kutta 4 (b) Méthode de Runge Kutta.</caption>

### 3.1.5.5 Méthodes d'ordre 5

<table><caption>Tableau 5 Méthode de Lawson [Lawson 1966].</caption><tr><td>0</td><td></td><td></td></tr><tr><td>1/2</td><td>1/2</td><td></td></tr><tr><td>1/4</td><td>1/8</td><td>1/8</td></tr><tr><td>1/2</td><td>0</td><td>0</td><td>1/2</td></tr><tr><td>3/4</td><td>0</td><td>-3/16</td><td>6/16</td><td>9/16</td></tr><tr><td>1</td><td>1/7</td><td>4/7</td><td>67</td><td>-12/7</td><td>8/7</td></tr><tr><td></td><td>7/90</td><td>0</td><td>32/90</td><td>12/90</td><td>32/90</td><td>7/90</td></tr></table>

<table><caption>Tableau 6 Méthode de Fehlberg [Fehlberg 1974].</caption><tr><td>0</td><td></td><td></td></tr><tr><td>1/5</td><td>1/5</td><td></td></tr><tr><td>3/10</td><td>3/40</td><td>9/40</td></tr><tr><td>3/5</td><td>3/10</td><td>-9/10</td><td>6/5</td></tr><tr><td>9/10</td><td>3/40</td><td>27/40</td><td>-3/5</td><td>3/4</td></tr><tr><td>1</td><td>107/162</td><td>5/2</td><td>-140/27</td><td>35/9</td><td>-70/81</td></tr><tr><td></td><td>8/81</td><td>0</td><td>25/63</td><td>25/108</td><td>25/81</td><td>-1/28</td></tr></table>
### 3.1.5.6 Méthodes d'ordre 6

<table><caption>Tableau 7 Méthodes de Butcher pour $a = \pm\sqrt{5}$ [Butcher 1964].</caption><tr><td>1/12</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0</td><td>(5-a)/10</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0</td><td>-a/10</td><td>(5+2a)/10</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>0</td><td>(-15+7a)/2</td><td>(-1+a)/4</td><td>(15-7a)/10</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>5/12</td><td>(5-a)/60</td><td>0</td><td>1/6</td><td>(15+7a)/60</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>5/12</td><td>(5+a)/60</td><td>0</td><td>(9-5a)/12</td><td>1/6</td><td>(-5+3a)/10</td><td></td><td></td><td></td><td></td></tr><tr><td>1/12</td><td>1/6</td><td>0</td><td>(-55+25a)/12</td><td>(-25-7a)/12</td><td>5-2a</td><td>(5+a)/2</td><td></td><td></td><td></td></tr><tr><td></td><td>0</td><td>(5-a)/10</td><td>(5+a)/10</td><td>(5-a)/10</td><td>(5+a)/10</td><td>(5-a)/10</td><td>1</td><td></td><td></td></tr></table>

### 3.1.6 Récapitulatif sur les méthodes explicites à un pas

L'algorithme de Runge-Kutta permet de construire une méthode d'intégration d'un ordre *p* donné. Par contre, plus l'ordre *p* est élevé, plus le nombre *s* d'étages de la méthode est grand. Le tableau 8 synthétise les nombres *s* d'étages correspondants à plusieurs ordres de consistance *p*. On remarque que l'ordre 4 constitue un ordre charnière après lequel le nombre d'étages devient plus grand que l'ordre, induisant plus d'évaluations intermédiaires de la fonction *f*(*t*, *y*). Cette propriété contribue à la popularité de la méthode d'intégration explicite de Runge-Kutta d'ordre 4 (i.e. méthode ERK4).

<table><caption>Tableau 8 Récapitulatif des ordres de consistance *p* et des nombres *s* d'étages correspondants.</caption><tr><td><i>p</i></td><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td><td>7</td><td>8</td><td>9</td><td>10</td></tr><tr><td><i>S</i></td><td>1</td><td>2</td><td>3</td><td>4</td><td>6</td><td>7</td><td>9</td><td>11</td><td>12-17</td><td>13-17</td></tr></table>

Les méthodes explicites de Runge-Kutta ne sont pas $\lambda$-h-stables. Ces méthodes ne sont donc pas adaptées à la résolution de problèmes raides. La Figure 18 présente les domaines de $\lambda$-h-stabilité des méthodes d'Euler, de Runge et de Runge-Kutta 4. Il apparaît que la méthode de Runge-Kutta 4 a le plus grand domaine de stabilité, un autre argument positif en sa faveur.

<figure><img src="image_27.png" /><figcaption>Figure 18 Domaines de $\lambda$-h-stabilité des méthodes d'Euler, de Runge (i.e. point milieu) et de Runge-Kutta 4.</figcaption></figure>
### 3.1.7 Contrôle de l'erreur, sélection du pas

Le choix du pas de temps lors de la résolution numérique est déterminant sur l'ordre de grandeur de l'erreur globale. En utilisant un pas de temps très faible, on s'assure en général une bonne convergence de la méthode bien que celui-ci, à certains temps de l'intégration, pourrait être plus grand sans pour autant altérer la précision de la solution. Il apparaît donc primordial de mettre en place un mécanisme d'ajustement automatique du pas de temps, garantissant une majoration de l'erreur par une tolérance donnée TOL.

L'extrapolation de Richardson (2.35) est une stratégie classique et générique d'ajustement automatique du pas de temps. Ce mécanisme consiste en un premier calcul de la solution approchée $y_{n+1}$ pour un pas de temps $h$. Puis on calcule une solution approchée $\tilde{y}_{n+1/2}$ pour un pas de temps $h/2$ qui sert de donnée initiale pour le calcul de $\tilde{y}_{n+1}$ pour un pas de temps $h/2$. La solution $\tilde{y}_{n+1}$ issue d'un calcul avec un pas de temps plus petit est plus précise que $y_{n+1}$.

$$ 
\begin{aligned}
y_{n+1} &= y_n + h_n \varphi(t_n, y_n, h_n) \\
\tilde{y}_{n+1/2} &= y_n + \frac{h_n}{2} \varphi\left(t_n, y_n, \frac{h_n}{2}\right) \rightarrow \tilde{y}_{n+1} = \tilde{y}_{n+1/2} + \frac{h_n}{2} \varphi\left(t_n + \frac{h_n}{2}, \tilde{y}_{n+1/2}, \frac{h_n}{2}\right)
\end{aligned} 
\quad (2.35) 
$$

la différence $ERR = |y_{n+1} - \tilde{y}_{n+1}|$ fournit une mesure de l'erreur globale qui permet, si besoin est, de diminuer le pas de temps de façon efficace grâce à la formule suivante [Press 1988] :

$$ 
\tilde{h}_n = \begin{cases} \beta h_n \left( \frac{TOL}{ERR} \right)^{1/p} & \text{si } ERR > TOL \\ \beta h_n \left( \frac{TOL}{ERR} \right)^{1/(p+1)} & \text{sinon} \end{cases} \quad \text{avec } 0 < \beta < 1 \quad (2.36) 
$$

$p$ est l'ordre de la méthode d'intégration.

L'extrapolation de Richardson est un mécanisme d'ajustement automatique du pas de temps efficace mais nécessite 3s évaluations de la fonction $f$ pour une méthode explicite de Runge-Kutta à $s$ étages. L'évaluation de la fonction $f$, correspondant aux efforts appliqués au système dans le cadre de la résolution des équations du mouvement, est en règle générale coûteuse. Cet état de fait a motivé la recherche de méthodes d'ajustement du pas de temps plus optimales.

J. A. Zonneveld [Zonneveld 1964], exploitant les particularités des méthodes de Runge-Kutta, proposa un mécanisme d'ajustement du pas de temps beaucoup plus performant. Il choisit d'utiliser en fait deux méthodes ERK emboîtées d'ordres respectifs $p$ et $q$ à $s$ et $s+1$ étages, partageant les mêmes poids intérieurs (Tableau 9). La différence des valeurs $y_{n+1}$ et $\tilde{y}_{n+1}$ déduites des deux intégrations simultanées avec le même pas de temps $h$ fournit aussi une mesure de l'erreur globale. Les deux méthodes partageant les mêmes poids intérieurs, le calcul des valeurs intermédiaires de la fonction $f$ n'est effectué que $s$ fois, contrairement aux 3s évaluations requises par le mécanisme proposé par Richardson pour une méthode de Runge-Kutta à $s$ étages.
<table><caption>Tableau 9 Couple de méthodes ERK emboîtées à s et s+1 étages.</caption><tbody><tr><td>c<sub>1</sub></td><td colspan="4">a<sub>2,1</sub></td><td></td><td></td></tr><tr><td>c<sub>3</sub></td><td>a<sub>3,1</sub></td><td>a<sub>3,2</sub></td><td></td><td></td><td></td><td></td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td></td><td></td><td></td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td></td><td></td></tr><tr><td>c<sub>s</sub></td><td>a<sub>s,1</sub></td><td>a<sub>s,2</sub></td><td>...</td><td>...</td><td>a<sub>s,s-1</sub></td><td></td></tr><tr><td>1</td><td>b<sub>1</sub></td><td>b<sub>2</sub></td><td>...</td><td>...</td><td>b<sub>s-1</sub></td><td>b<sub>s</sub></td></tr><tr><td></td><td>b'<sub>1</sub></td><td>b'<sub>2</sub></td><td>...</td><td>...</td><td>b'<sub>s-1</sub></td><td>b'<sub>s</sub></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>b'<sub>s+1</sub></td></tr></tbody></table>

Je présente maintenant plusieurs couples de méthodes ERK emboîtées appelées **paires (p,q)** où *p* et *q* sont les ordres respectifs de chaque méthode, classés par ordre de consistance croissant [Stroud 1974] [Crouzeix 1989] [Enright 1995].

<table><caption>Tableau 10 Une paire (1,2) de méthodes ERK.</caption><tbody><tr><td>0</td><td></td><td></td><td></td></tr><tr><td>1</td><td>1</td><td></td><td></td></tr><tr><td></td><td>1/2</td><td>1/2</td><td></td></tr></tbody></table>

<table><caption>Tableau 11 Une paire (2,4) de méthodes ERK.</caption><tbody><tr><td>0</td><td></td><td></td><td></td></tr><tr><td>1/2</td><td>1/2</td><td></td><td></td></tr><tr><td>1/2</td><td>0</td><td>1/2</td><td></td></tr><tr><td>1</td><td>0</td><td>0</td><td>1</td></tr><tr><td></td><td>1/6</td><td>2/6</td><td>2/6</td></tr><tr><td></td><td></td><td></td><td>1/6</td></tr></tbody></table>

<table><caption>Tableau 12 Une paire (3,4) de méthodes ERK.</caption><tbody><tr><td>0</td><td></td><td></td><td></td></tr><tr><td>2/7</td><td>2/7</td><td></td><td></td></tr><tr><td>4/7</td><td>-8/35</td><td>4/5</td><td></td></tr><tr><td>6/7</td><td>29/42</td><td>-2/3</td><td>5/6</td></tr><tr><td>1</td><td>1/6</td><td>1/6</td><td>5/12</td></tr><tr><td></td><td>11/96</td><td>7/24</td><td>35/96</td></tr><tr><td></td><td></td><td></td><td>7/48</td></tr><tr><td></td><td></td><td></td><td>1/12</td></tr></tbody></table>

## 3.2 Méthodes implicites

Toute méthode d'intégration implicite à un pas est de la forme suivante :

$$y_{n+1} = y_n + h_n \varphi(t_n, y_n, f_{n+1}, h_n) \quad \text{où} \quad \varphi: [t_0, t_0+T] \times \mathbb{R}^n \times \mathbb{R}^n \times \mathbb{R} \to \mathbb{R}^n \quad (2.37)$$

Contrairement à une méthode d'intégration explicite (2.17), la méthode implicite à un pas requiert l'évaluation du terme $f_{n+1}$ qui équivaut à $f(t_{n+1}, y_{n+1})$. Cette particularité confère en général aux méthodes d'intégration implicites une plus grande stabilité, même dans le cadre des problèmes raides. Alors que les méthodes d'intégration explicites extrapolent la valeur $y_{n+1}$ de
la solution approchée en $t_{n+1}$ en fonction de la tendance en $t_n$, les méthodes implicites déterminent la valeur $y_{n+1}$ afin qu'elle soit compatible avec les efforts exprimés en $t_{n+1}$. Mais cette stabilité relative offerte par les méthodes implicites a un prix. En effet, l'équation (2.37) est une équation de nature souvent non linéaire dont la résolution nécessite des méthodes bien plus coûteuses que celles mises en œuvre dans le cadre des méthodes explicites. Les différentes méthodes d'intégration implicites à un pas se distinguent par la stratégie choisie pour résoudre cette équation non linéaire. J'énumère maintenant les schémas de résolution les plus courants.

### 3.2.1 Résolution par itérations

En considérant que $f_{n+1} = f(t_{n+1}, y_{n+1})$, l'équation (2.37) devient :

$$y_{n+1} = y_n + h_n \varphi(t_n, y_n, f(t_{n+1}, y_{n+1}), h_n) \quad (2.38)$$

Si on définit la fonction $\Psi(u) = y_n + h_n \varphi(t_n, y_n, f(t_{n+1}, u), h_n) - u$, la résolution de l'équation non linéaire (2.38) se ramène au problème suivant :

$$\text{Trouver u tel que: } \Psi(u) = 0 \quad (2.39)$$

On peut alors envisager de trouver les racines de cette équation non linéaire grâce à des méthodes classiques telles que la méthode de bissection, la méthode de la sécante ou bien la méthode de Newton-Raphson [Press 1988]. Malheureusement, toutes ces méthodes nécessitent la réalisation de critères difficiles à garantir pour un problème quelconque (unicité de la solution dans l'intervalle de recherche par exemple). Et ces méthodes nécessitent un grand nombre d'itérations pour aboutir à une solution convenable. Pour toutes ces raisons, la résolution par itérations est très rarement utilisée dans un cadre généraliste. Pourtant, très récemment, Hauth et al [Hauth 2001] proposent d'utiliser une méthode qualifiée de « méthode de Newton simplifiée inexacte » pour la résolution de l'équation linéaire (2.39). En fait, la méthode de Newton repose sur plusieurs résolutions successives de systèmes linéaires faisant intervenir la jacobienne de la fonction $\Psi$. Pour la résolution de ces systèmes linéaires successifs, ils privilégient la méthode du gradient conjugué, méthode elle aussi itérative qui fournit une solution approchée justifiant ainsi le terme de méthode de Newton inexacte. De plus, pour des raisons de performance, ils imposent que la jacobienne intervenant dans les itérations de la méthode de Newton soit constante. Cette hypothèse nécessaire pour garantir la rapidité de la méthode n'est pas sans conséquences sur la précision de la solution. En fait, ils démontrent l'efficacité de leur démarche pour des systèmes matériels dont les efforts appliqués sont linéaires, par conséquent avec une jacobienne des efforts constante.

### 3.2.2 Linéarisation

Une autre façon d'entreprendre la résolution de l'équation (2.37) consiste à remplacer $f(t_{n+1}, y_{n+1})$ par son développement de Taylor à l'ordre 1. Cette opération revient à ne traiter que la partie linéaire du problème en négligeant les termes d'ordre plus élevé. On débouche alors sur une expression linéaire, plus simple à résoudre. Cette technique a été employée par Baraff [Baraff 1998] dans le cadre de l'animation de textiles qui est un domaine réputé pour engendrer des équations du mouvement raides, dues en l'occurrence aux efforts élastiques internes.
On considère la méthode d’intégration d’Euler implicite (2.40) :

$$y_{n+1} = y_n + h_n f(t_{n+1}, y_{n+1}) \quad (2.40)$$

On réalise alors un développement de Taylor à l'ordre 1 de la fonction $f(t, y)$ par rapport à la variable y.

$$y_{n+1} = y_n + h_n \left[ f(t_{n+1}, y_n) + \frac{\partial f}{\partial y}(t_{n+1}, y_n)(y_{n+1} - y_n) \right] \quad (2.41)$$

Cette hypothèse permet de déboucher sur le système linéaire suivant :

$$\left[ I - h_n \frac{\partial f}{\partial y}(t_{n+1}, y_n) \right] (y_{n+1} - y_n) = h_n f(t_{n+1}, y_n) \quad (2.42)$$

En utilisant la méthode du gradient conjugué pour chaque pas de temps, Baraff détermine en peu d'itérations les valeurs successives $y_{n+1}$ de la solution approchée. Il est important de mentionner que dans cette étude il avait délibérément fait l'hypothèse que les efforts représentés par la fonction $f(t, y)$ dérivaient d'un potentiel quadratique. Ainsi, les efforts étaient de nature linéaire. Dans ce cas particulier, le mécanisme de linéarisation n'a pas d'effet « destructeur » sur la fonction d'efforts $f(t, y)$. La particularité de la méthode d'intégration d'Euler implicite, due à sa faible consistance et sa grande stabilité, réside dans le fait que la solution approchée calculée subit un certain amortissement. Cette propriété permettait en l'occurrence d'atténuer les oscillations caractéristiques de masses-ressorts et de produire des animations de textiles très convaincantes en utilisant un pas de temps d'intégration assez grand.

Dans le même ordre d'idée, Desbrun [Desbrun 1999] propose une optimisation de la méthode de Baraff en supposant la matrice du système linéaire (2.42) constante. En considérant le champ de forces stationnaire uniforme (i.e. indépendant du temps et de la position), Desbrun s'éloigne plus encore d'une démarche rigoureuse de modélisation du phénomène physique mais arrive à générer des animations tout à fait convaincantes et ce dans le cadre d'une application de réalité virtuelle dont les contraintes temps réel sont très difficiles à réaliser en général.

Plus récemment, Volino [Volino 2000] s'appuyant sur la précision assez médiocre de la méthode d'Euler, propose une méthode d'intégration implicite déduite de la méthode du point milieu pour résoudre les équations du mouvement d'un système masses-ressorts. Cette méthode d'intégration d'ordre 2 a les mêmes propriétés de stabilité que la méthode d'Euler implicite tout en étant plus précise.

L'engouement récent pour les méthodes d'intégration implicites à un pas par troncature est tout de même à tempérer. Bien que ces méthodes soient stables (en effet les deux méthodes citées ci-dessus sont $\lambda h$-stables), leur efficacité n'est flagrante que pour le cas particulier d'une fonction d'efforts $f(t, y)$ linéaire. Dans un cadre général, ces méthodes peuvent s'avérer inadaptées à cause de l'élimination des termes d'ordres supérieurs à 2 dans l'expression de $f(t, y)$. De plus, elles reposent sur l'expression explicite de la jacobienne de $f(t, y)$ dont on ne dispose pas a priori
### 3.2.3 Méthode Prédiction-Evaluation-Correction

Cette dernière approche est sûrement la plus générique puisqu'elle peut être appliquée à une méthode d'intégration implicite quelconque à un pas d'ordre $q$, de fonction $\varphi^q$. En effet, on commence par calculer une valeur $\tilde{y}_{n+1}$ de la solution à l'aide d'une méthode d'intégration explicite d'ordre $p$ ($p<q$) de fonction $\varphi^p$: cette étape est qualifiée de **Prédiction**. De manière classique, la valeur $\tilde{y}_{n+1}$ consiste en une extrapolation de la solution en fonction des tendances à l'instant $t_n$. Mais cette valeur intermédiaire $\tilde{y}_{n+1}$, après **Evaluation** (i.e. calcul de $f(t_{n+1}, \tilde{y}_{n+1})$), est ensuite injectée dans le terme de droite de la méthode d'intégration implicite (2.37) : cette étape appelée **Correction** permet d'aboutir à la valeur de $y_{n+1}$.

Ce mécanisme « casse » en quelque sorte la nature implicite de la méthode en recourant à une méthode d'intégration explicite. Ainsi, pour garantir la stabilité de la méthode on peut avoir recours à plusieurs étapes de Prédiction-Correction afin de converger vers une valeur approchée satisfaisante (2.43). Du point de vue méthodologique, ces itérations peuvent être comparées aux itérations requises par la méthode du gradient conjugué dans le cadre des problèmes linéarisés.

$$ P \overbrace{\overbrace{\overbrace{P}^{n \text{ fois}} ECEC \cdots EC}^{\cdots}}^{\cdots E} $$

$P$: prédiction

$$ \tilde{y}_{n+1}^1 = y_n + h_n \varphi^p(t_n, y_n, h_n) \qquad (2.43) $$

$E$: évaluation

$$ \tilde{f}_{n+1}^i = f(t_n, \tilde{y}_{n+1}^i) \qquad 1 \le i \le n $$

$C$: correction

$$ \tilde{y}_{n+1}^{i+1} = y_n + h_n \varphi^q(t_n, y_n, \tilde{f}_{n+1}^i, h_n) $$

## 4 Méthodes d'intégration à pas multiples

Parallèlement aux travaux menés pour la construction de méthodes d'intégration à un pas, des scientifiques (notamment J.C. Adams) ont mis en place des méthodes de résolution s'appuyant sur la connaissance des valeurs de la solution approchée calculées antérieurement. Les méthodes présentées dans cette section reposent sur la substitution de la fonction $f(t, y)$ par un polynôme interpolant les valeurs $(f_{n-i})_{0 \le i \le r}$ [Crouzeix 1989] [Demailly 1991] [Hairer 1993]. Contrairement aux méthodes d'intégration à un pas qui ont recours à des évaluations intermédiaires de la fonction $f(t, y)$ pour calculer la valeur $y_{n+1}$ de la solution approchée, les méthodes d'intégration à pas multiples s'appuient sur la connaissance des valeurs antérieures de $f(t, y)$. Ainsi, à ordre de consistance équivalent, les méthodes à pas multiples sont beaucoup plus rapides que les méthodes à un pas. Cependant, l'usage d'un nombre important de valeurs antérieures pour garantir la consistance peut nuire à la stabilité de la solution déduite d'un polynôme de degré élevé.
Je traite ici des méthodes d'intégration (explicites ou implicites) à r+1 pas de la forme :

$$ \sum_{i=0}^{r+1} \alpha_i y_{n+1-i} = h_n \sum_{i=0}^{r+1} \beta_i f_{n+1-i} \quad \text{où } \alpha_0 \neq 0 \text{ et } |\alpha_{r+1}| + |\beta_{r+1}| > 0 \quad (2.44) $$

Si $\beta_0 = 0$, la méthode d'intégration est explicite.

# 4.1 Méthodes explicites

## 4.1.1 Méthode d'Adams-Bashforth (1855)

Adams choisit d'aborder le problème de la même façon que Runge et Kutta, choisissant comme point de départ la relation (2.27). L'originalité de sa méthode réside dans le fait de remplacer la fonction $f$ par son polynôme d'interpolation aux points $(t_{n-i}, f_{n-i})_{0 \le i \le r}$.

On rappelle la définition du **polynôme d'interpolation de Lagrange** :

$$ p_{n,r}(t) = \sum_{i=0}^{r} f_{n-i} L_{n,i,r}(t) \quad \text{où} \quad L_{n,i,r}(t) = \prod_{\substack{0 \le j \le r \\ j \ne i}} \frac{t - t_{n-j}}{t_{n-i} - t_{n-j}} \quad (2.45) $$

L'expression (2.28) devient alors :

$$ \begin{aligned} y(t_{n+1}) &\approx y(t_n) + \int_{t_n}^{t_{n+1}} p_{n,r}(u) du \\ &\approx y(t_n) + h_n \sum_{i=0}^{r} b_{n,i,r} f_{n-i} \quad \text{où} \quad b_{n,i,r} = \frac{1}{h_n} \int_{t_n}^{t_{n+1}} L_{n,i,r}(u) du \end{aligned} \quad (2.46) $$

D'où la méthode suivante :

$$ \left\{ \begin{array}{l} t_{n+1} = t_n + h_n \\ y_{n+1} = y_n + h_n \sum_{i=0}^{r} b_{n,i,r} f_{n-i} \end{array} \right. \quad (2.47) $$

La méthode énoncée en (2.47) est appelée méthode d'Adams Bashforth à r+1 pas (i.e. méthode $AB_{r+1}$).

Dans le cadre d'une résolution à pas de temps adaptatif, les termes $b_{n,i,r}$ qui dépendent de la longueur de l'intervalle $[t_n, t_{n+1}]$ doivent être évalués de nouveau à chaque itération. Par contre, pour une résolution à pas fixe, ces mêmes termes deviennent constants permettant ainsi de déboucher sur une expression simplifiée des méthodes $AB_{r+1}$.
Je cite en exemple quelques méthodes d'Adams-Bashforth de 1 à 6 pas [Stroud 1974] :

$$r=0: y_{n+1} = y_n + hf_n$$

$$r=1: y_{n+1} = y_n + h\left(\frac{3}{2}f_n - \frac{1}{2}f_{n-1}\right)$$

$$r=2: y_{n+1} = y_n + h\left(\frac{23}{12}f_n - \frac{16}{12}f_{n-1} + \frac{5}{12}f_{n-2}\right)$$

$$r=3: y_{n+1} = y_n + h\left(\frac{55}{24}f_n - \frac{59}{24}f_{n-1} + \frac{37}{24}f_{n-2} - \frac{9}{24}f_{n-3}\right)$$

$$r=4: y_{n+1} = y_n + h\left(\frac{1901}{720}f_n - \frac{2774}{720}f_{n-1} + \frac{2616}{720}f_{n-2} - \frac{1274}{720}f_{n-3} + \frac{251}{720}f_{n-4}\right)$$

$$r=5: y_{n+1} = y_n + h\left(\frac{4277}{1440}f_n - \frac{7923}{1440}f_{n-1} + \frac{9982}{1440}f_{n-2} - \frac{7298}{1440}f_{n-3} + \frac{2877}{1440}f_{n-4} - \frac{475}{1440}f_{n-5}\right)$$

Une méthode $AB_{r+1}$ est d'ordre $r+1$. A condition d'être en mesure de calculer les $r$ premières valeurs de la solution approchée, ces méthodes semblent performantes dans la mesure où elles ne requierent pas d'évaluations intermédiaires de la fonction $f(t, y)$ comme les méthodes ERK tout en offrant une bonne consistance. Seule ombre au tableau : les méthodes $AB_{r+1}$ sont peu stables. En effet, les domaines de $\lambda$ h-stabilité de ces méthodes sont plus petits que ceux des méthodes ERK de même ordre de consistance (Figure 19). Ce phénomène peut être en partie expliqué par le caractère oscillant des polynômes d'interpolation pour des valeurs de $r$ de plus en plus grandes.

<figure><img src="image_28.png" /><figcaption>Figure 19 Domaines de $\lambda$ h-stabilité des méthodes d'Adams-Bashforth d'ordre $k=1, 2, 3, 4$.</figcaption></figure>

## 4.2 Méthodes implicites

### 4.2.1 Méthode d'Adams-Moulton (1855)

Cette fois-ci, on remplace $f(t, y)$ par son polynôme d'interpolation en $(t_{n+1-i}, f_{n+1-i})_{0 \le i \le r}$

$$p_{n,r}^*(t) = \sum_{i=-1}^{r} f_{n-i} L_{n,i,r}^*(t) \quad \text{où} \quad L_{n,i,r}^*(t) = \prod_{-1 \le j \le r \atop j \ne i} \frac{t - t_{n-j}}{t_{n-i} - t_{n-j}} = \frac{t - t_{n+1}}{t_{n-i} - t_{n+1}} L_{n,i,r}(t) \quad (2.48)$$
En utilisant la même démarche que précédemment, on débouche sur la méthode suivante :

$$
\begin{equation} \tag{2.49}
\begin{cases}
t_{n+1} = t_n + h_n \\
y_{n+1} = y_n + h_n \sum_{i=-1}^{r} b_{n,i,r}^* f_{n-i} \quad \text{où} \quad b_{n,i,r}^* = \frac{1}{h_n} \int_{t_n}^{t_{n+1}} L_{n,i,r}^*(u) du
\end{cases}
\end{equation}
$$

La méthode énoncée en (2.49) est appelée méthode d'Adams Moulton à r+1 pas (i.e. méthode AM<sub>r+1</sub>). Lors d'une résolution à pas fixe, les termes b<sup>*</sup><sub>n,i,r</sub> sont constants permettant ainsi de donner une formulation simplifiée des méthodes AM<sub>r+1</sub>. Afin d'illustrer cette démarche, j'expose quelques exemples de méthodes d'Adams-Moulton de 1 à 6 pas [Stroud 1974] :

$$
r = 0: \quad y_{n+1} = y_n + hf_{n+1}
$$

$$
r = 1: \quad y_{n+1} = y_n + h \left( \frac{1}{2} f_{n+1} + \frac{1}{2} f_n \right)
$$

$$
r = 2: \quad y_{n+1} = y_n + h \left( \frac{5}{12} f_{n+1} + \frac{8}{12} f_n - \frac{1}{12} f_{n-1} \right)
$$

$$
r = 3: \quad y_{n+1} = y_n + h \left( \frac{9}{24} f_{n+1} + \frac{19}{24} f_n - \frac{5}{24} f_{n-1} + \frac{1}{24} f_{n-2} \right)
$$

$$
r = 4: \quad y_{n+1} = y_n + h \left( \frac{251}{720} f_{n+1} + \frac{646}{720} f_n - \frac{264}{720} f_{n-1} + \frac{106}{720} f_{n-2} - \frac{19}{720} f_{n-3} \right)
$$

$$
r = 5: \quad y_{n+1} = y_n + h \left( \frac{475}{1440} f_{n+1} + \frac{1427}{1440} f_n - \frac{798}{1440} f_{n-1} + \frac{482}{1440} f_{n-2} - \frac{173}{1440} f_{n-3} + \frac{27}{1440} f_{n-4} \right)
$$

Une méthode d'Adams-Moulton à *r*+1 pas est d'ordre *r*+2. De plus, grâce à leur nature implicite, les méthodes AM<sub>*r*+1</sub> sont plus stables que les méthodes AB<sub>*r*+1</sub> (Figure 20). Par contre l'équation non linéaire induite par les méthodes implicites AM<sub>*r*+1</sub> nécessite la mise en œuvre d'une des stratégies mentionnées dans la section 3.2. Classiquement, on utilise un schéma PEC...EC en choisissant une méthode AB<sub>*r*+1</sub> pour l'étape de prédiction.

<figure><img src="image_29.png" /><figcaption>Figure 20 Domaines de $\lambda$. h-stabilité des méthodes d'Adams-Moulton d'ordre $k=1, 2, 3, 4$.</figcaption></figure>
### 4.2.2 La méthode BDF (Backward Differentiation Formulas)

Au lieu de construire un polynôme d'interpolation pour la fonction $f(t, y)$, la méthode BDF interpole directement les valeurs $(t_{n+1-i}, y_{n+1-i})_{0 \le i \le r}$ de la solution approchée.

$$ q_{n,r}^*(t) = \sum_{i=-1}^{r} y_{n-i} L_{n,i,r}^*(t) \quad \text{où} \quad L_{n,i,r}^*(t) = \prod_{-1 \le j \le r, j \ne i} \frac{t - t_{n-j}}{t_{n-i} - t_{n-j}} \qquad (2.50) $$

$$ \frac{dq_{n,r}^*}{dt}(t) = \sum_{i=-1}^{r} y_{n-i} \frac{dL_{n,i,r}^*}{dt}(t) \qquad (2.51) $$

Puisque le polynôme $q_{n,r}^*(t)$ interpole les valeurs de la solution approchée du problème de Cauchy (2.1), on considère qu'il est presque solution du même problème. Ainsi,

$$ \frac{dq_{n,r}^*}{dt}(t) \approx f(t, q_{n,r}^*(t)) \Rightarrow \frac{dq_{n,r}^*}{dt}(t_{n+1}) \approx f(t_{n+1}, y_{n+1}) \text{ puisque } q_{n,r}^*(t_{n+1}) = y_{n+1} \quad (2.52) $$

On aboutit alors à la méthode d'intégration suivante :

$$ \left\{ \begin{array}{l} t_{n+1} = t_n + h_n \\ \displaystyle \sum_{i=-1}^{r} y_{n-i} \frac{dL_{n,i,r}^*}{dt}(t_{n+1}) = f(t_{n+1}, y_{n+1}) \end{array} \right. \qquad (2.53) $$

Lorsque l'intégration numérique se fait à pas de temps fixe, on peut procéder au précalcul des dérivées des termes $L_{n,i,r}^*(t_{n+1})$. Je présente quelques exemples de méthodes BDF à 1,2,3 et 4 pas [Hairer 1993] :

$$ r = 0: y_{n+1} - y_n = hf_{n+1} $$

$$ r = 1: \frac{3}{2} y_{n+1} - 2y_n + \frac{1}{2} y_{n-1} = hf_{n+1} $$

$$ r = 2: \frac{11}{6} y_{n+1} - 3y_n + \frac{3}{2} y_{n-1} - \frac{1}{3} y_{n-2} = hf_{n+1} $$

$$ r = 3: \frac{25}{12} y_{n+1} - 4y_n + 3y_{n-1} - \frac{4}{3} y_{n-2} + \frac{1}{4} y_{n-3} = hf_{n+1} $$

Une méthode BDF à $r+1$ pas est d'ordre $r+1$. De plus, en raison de son schéma de construction, elle est particulièrement stable. Cette méthode est donc à privilégier pour la résolution de problèmes raides.

## 5 Synthèse

Comme on peut facilement le constater, il semble illusoire de trouver « la » méthode d'intégration numérique générale et efficace au vu de la grande diversité des équations à résoudre. On doit alors se contenter de choisir la méthode la plus appropriée au problème que l'on souhaite traiter.
<table><caption>Tableau 13. Récapitulatif des propriétés des méthodes d’intégration classiques étudiées.</caption><thead><tr><th>Méthode</th><th>Ordre</th><th>Stabilité</th><th>Evaluations</th><th>Généricité</th></tr></thead><tbody><tr><td>Euler</td><td>1</td><td>---</td><td>1</td><td>oui</td></tr><tr><td>Point milieu</td><td>2</td><td>--</td><td>2</td><td>oui</td></tr><tr><td>Heun</td><td>3</td><td>--</td><td>3</td><td>oui</td></tr><tr><td>Runge Kutta 4</td><td>4</td><td>-</td><td>4</td><td>oui</td></tr><tr><td>Lawson</td><td>5</td><td>--</td><td>6</td><td>oui</td></tr><tr><td>Fehlberg</td><td>5</td><td>--</td><td>6</td><td>oui</td></tr><tr><td>Butcher</td><td>6</td><td>---</td><td>7</td><td>oui</td></tr><tr><td>Euler (linéarisation)</td><td>1</td><td>+++</td><td>1+</td><td>non</td></tr><tr><td>Point milieu (linéarisation)</td><td>2</td><td>+++</td><td>2+</td><td>non</td></tr><tr><td>AB<sub>r+1</sub></td><td>r+1</td><td>--</td><td>1</td><td>oui</td></tr><tr><td>AM<sub>r+1</sub></td><td>P(EC)<sup>n</sup>E</td><td>-</td><td>n+1</td><td>oui</td></tr><tr><td>BDF<sub>r+1</sub></td><td>P(EC)<sup>n</sup>E</td><td>++</td><td>n+1</td><td>oui</td></tr></tbody></table>

A cette fin, je récapitule dans le tableau 13 les méthodes exposées dans ce chapitre en mentionnant leurs caractéristiques importantes :

* L'ordre de la méthode renseigne sur la précision de la méthode d'intégration

* La **stabilité** est une indication sur la capacité de la méthode à converger même face à un problème raide.

* Le nombre d'évaluations de la fonction $f(t, y)$ pour une itération.

* La **généricité** indique si la méthode peut s'appliquer à tout type de problèmes ou bien si elle procède à une simplification des termes de l'équation, négligeant certaines quantités.

Les lignes grisées correspondent aux méthodes d'intégration explicites et les lignes blanches aux méthodes d'intégration implicites.

Le nombre d'évaluations pour les méthodes implicites d'Euler et du point milieu sont marqués du signe + car aux évaluations de la fonction $f(t, y)$ il faut ajouter les évaluations de sa jacobienne.
3

CORPS DEFORMABLES EN
ANIMATION DYNAMIQUE.
ETAT DE L'ART

# 1 Introduction

L'animation d'objets déformables constitue en synthèse d'images un véritable défi. En effet, aux difficultés de produire un mouvement global cohérent de ces objets viennent s'ajouter les difficultés de représenter leurs changements complexes de forme. Dès le début des années quatre vingt, les travaux de recherche en animation par ordinateur ont eu recours aux lois de la physique pour représenter les déformations subies par un objet de synthèse. L'utilisation des lois de la dynamique permet de concevoir des animations convaincantes mais surtout assure l'automatisation de la production de celles ci. En effet si l'on confère aux objets de synthèse des propriétés mécaniques, ils évoluent en fonction des contraintes physiques qui leurs sont imposées et ce conformément aux équations du mouvement.

En physique, l'évolution des corps déformables est décrite par la mécanique des milieux continus. Dans ce formalisme, les déformations d'un corps supposé continu dépendent d'une infinité de paramètres. En description lagrangienne, ces paramètres correspondent aux positions dans l'espace de l'infinité des particules qui constituent le corps matériel alors qu'en description eulerienne, ils correspondent à l'infinité des valeurs prises par le champ de vitesses instantanées. Afin d'être en mesure d'exprimer et surtout de résoudre les équations qui régissent le mouvement d'un corps déformable, on doit recourir à une approximation méthodologique qui consiste à réduire le nombre des paramètres qui conditionnent l'évolution du corps. Cette démarche classique en mécanique revient à réduire l'ensemble des configurations d'intérêt garantissant ainsi la finitude des paramètres caractéristiques. Les différents modèles utilisés en synthèse d'images se distinguent par la démarche de réduction de l'ensemble des configurations acceptables mais aussi, en conséquence, par le choix des équations du mouvement. Dans ce chapitre, j'expose les modèles mécaniques les plus couramment utilisés en synthèse d'images, classés en fonction des lois physiques sur lesquelles ils reposent. Puis, j'introduis le modèle innovant que je propose en le comparant aux modèles précédents.
# 2 Modèles classiques

## 2.1 Mécanique du point matériel

Les modèles considérés dans cette section s'appuient sur l'hypothèse suivante : le corps déformable étudié est représenté par un **ensemble discret de particules** noté $(\mathbf{x}^i(t))_{i \le i \le n}$. La masse totale du corps est répartie de façon discontinue parmi les particules qui le constituent. Et le mouvement du corps est déduit de façon émergente des mouvements individuels de chaque particule.

La loi fondamentale de la dynamique régit le mouvement d'une particule élémentaire représentée par sa position $\mathbf{x}(t)$, de masse $m$, soumise à des forces de résultante $\mathbf{F}(\mathbf{x}(t), \dot{\mathbf{x}}(t), t)$.

$$ \frac{d}{dt}(m\dot{\mathbf{x}}(t)) = \mathbf{F}(\mathbf{x}(t), \dot{\mathbf{x}}(t), t) \qquad (3.1) $$

Les modèles reposant sur cette approche particulière se distinguent en fait par la représentation des interactions entre particules.

### 2.1.1 Modèle masses-ressorts

Afin de conférer à l'objet déformable une certaine cohésion interne, les particules matérielles (réparties en un réseau uni, bi ou tridimensionnel généralement régulier) sont connectées par des ressorts à leurs particules « voisines » (Figure 21). Pour une particule de position $\mathbf{x}^i(t)$, les ressorts connectés aux particules voisines $\mathbf{x}^j(t)$ pour $j \in V^i$ génèrent des efforts de la forme :

$$ \sum_{j \in V^i} \mathbf{F}^{ij}(t) \qquad (3.2) $$

$\mathbf{F}^{ij}(t)$ correspond à la réponse du ressort d'extrémités $\mathbf{x}^i(t)$, $\mathbf{x}^j(t)$.

Afin de réduire le comportement oscillatoire du réseau masses-ressorts, on ajoute une force d'amortissement de coefficient $\gamma^i$, proportionnelle à la vitesse $\dot{\mathbf{x}}^i(t)$ de la particule $i$. L'équation fondamentale de la dynamique devient alors, dans le cadre du modèle masses-ressorts,

$$ \frac{d}{dt}(m\dot{\mathbf{x}}^i(t)) = \sum_{j \in V^i} \mathbf{F}^{ij}(t) - \gamma^i \dot{\mathbf{x}}^i(t) + \mathbf{F}(\mathbf{x}(t), \dot{\mathbf{x}}(t), t) \qquad (3.3) $$

Ce modèle fut utilisé très tôt par Platt et Badler [Platt 1981] pour reproduire des expressions du visage. Des points caractéristiques du visage sont reliés par plusieurs ressorts, de sorte que, lorsque plusieurs points sont déplacés, l'ensemble du réseau s'agence en une nouvelle configuration conforme au déplacement subi. Luciani [Luciani 1985] utilisa ce même modèle pour mettre en œuvre un outil d'animation de marionnettes, tandis que Miller [Miller 1988] s'en servit pour animer des reptiles en ajoutant des contraintes de préservation de volume et un modèle de propagation des ondes sinusoïdales au travers du réseau. Plus récemment, le modèle masses-ressorts fut utilisé pour la modélisation de textiles tissés. Provot [Provot 1995], en imposant une distance maximale d'élongation, garantit une réponse non linéaire des ressorts.
Boux de Casson [B. de Casson 2000], autorise l'élimination de certaines connections dans le réseau pour simuler des déchirures.

<figure><img src="image_30.png" /><figcaption>Figure 21 Réseau masses-ressorts bidimensionnel.</figcaption></figure>

En raison de sa simplicité, le modèle masses-ressorts a connu un vif succès au sein de la communauté de l'animation par ordinateur. D'implémentation rapide et potentiellement parallélisable, ce type de modèle permet d'obtenir des animations de bonne facture avec une certaine économie de moyens. Par contre, ce type de modélisation, tout comme les autres d'ailleurs, devient problématique lorsque l'on cherche à représenter des matériaux assez rigides. En effet, à mesure que l'on augmente la raideur *k* des ressorts liant des particules de masses *m*, le pas de temps *h* décroît conformément à la règle de Coulomb.

$$h = \sqrt{\frac{m}{k}} \qquad (3.4)$$

Baraff [Baraff 1998] propose alors d'utiliser une méthode d'intégration d'Euler implicite pour résoudre le système d'équations différentielles (3.3) réputé pour être un problème « raide » (vous pouvez vous reporter au chapitre consacré à l'intégration numérique pour une définition de ce terme). Bien que la méthode de Baraff garantisse la stabilité de la solution même pour des pas de temps élevés, son domaine d'application est restreint au cas d'efforts linéaires. Dans le même état d'esprit, Volino étend la démarche au cas des efforts quadratiques en utilisant une méthode d'intégration implicite dérivée de la méthode du point milieu [Volino 2000].

Un autre inconvénient du modèle masses-ressorts réside dans l'identification des paramètres. Dans une optique de simulation où l'on dispose des caractéristiques physiques du matériau, quelle structure de réseau privilégier ? Quelles raideurs choisir ? En guise de réponse, Van Gelder [Van Gelder 1998] a démontré l'impossibilité de régler les raideurs d'un modèle masses-ressorts afin de le faire coïncider avec un modèle continu de la matière.

### 2.1.2 Modèle à base de particules

A l'image du comportement microscopique, les interactions entre particules interviennent ici sous la forme de forces d'attraction répulsion. Le modèle à base de particules, à l'inverse du modèle masses-ressorts, n'est pas un modèle structuré : Cozot [Cozot 1996] utilise le terme de modèle lagrangien libre-libre pour désigner les modèles à base de particules pour les distinguer des modèles lagrangiens libres définis à partir d'un maillage discret évolutif. Les particules sont donc libres et peuvent potentiellement interagir avec toutes les autres particules, alors que le modèle masses-ressorts réduit le nombre d'interactions aux « voisins » de la particule au sein du réseau structuré. Dans la pratique, on a aussi recours à une notion de voisinage lors du calcul
des interactions entre particules afin de réduire l'importante combinatoire résultant d'un grand nombre de particules. Ici le voisinage n'est pas issu d'un graphe de connexion mais d'une distance inter-particules qui correspond à un seuillage des forces d'attraction répulsion (Figure 22). Le grand intérêt des systèmes de particules réside dans la nature a priori quelconque de l'énergie potentielle dont dérivent les forces d'attraction répulsion. Terzopoulos et al. [Terzopoulos 1989] simulent des matériaux en fusion en faisant varier l'intensité de la force d'interaction entre deux particules en fonction de la température, jusqu'à l'annuler pour rendre les particules libres de tout mouvement. Breen et al. [Breen 1992] mettent à profit cette flexibilité pour simuler le drapé d'un tissu à l'aide d'un système de particules. L'énergie potentielle du système est composée d'une énergie de répulsion pour empêcher le chevauchement de matière et d'énergies pour résister à la traction, à la flexion et au cisaillement. Bien que les résultats soient satisfaisants, la méthode se limite au cas de la statique et surtout, les coefficients de ces énergies sont complètement empiriques. Deux ans après, ils proposent une amélioration de leur modèle [Breen 1994] en utilisant des données expérimentales issues de tests de Kawabata [Kawabata 1980]. Volino et al. [Volino 1995] optent aussi pour un système de particules mais pour modéliser le mouvement des textiles. Ils déduisent leurs énergies de la théorie de l'élasticité des coques minces. Enfin, Eberhardt et al. [Eberhardt 1996] s'inspirent aussi des données expérimentales de Kawabata pour proposer un modèle dynamique de tissu à base de particules.

<figure><img src="image_31.png" /><figcaption>Figure 22 Système de particules et l'énergie potentielle d'attraction répulsion associée. La limite de la zone d'influence entre particules est symbolisée par le cercle en pointillés.</figcaption></figure>

Le modèle à base de particules qui ne souffre pas du problème crucial d'identification des paramètres pose tout de même un problème. En effet, quelle est la validité d'une telle approche, issue d'une observation des phénomènes microscopiques d'attraction répulsion, à l'échelle macroscopique ? De manière équivalente, quelle taille de particules choisir pour obtenir un modèle cohérent ?

### 2.1.3 Modèle SPH

Le modèle SPH [Benz 1988] qui signifie Smoothed Particle Hydrodynamics est une variante du modèle à base de particules. Il a été introduit par les physiciens pour simuler de façon plus
exacte la dynamique des fluides. Ici, les particules sont désormais considérées comme des
positions d'échantillonnage de quantités continues, caractéristiques du milieu simulé. Une
expression continue de ces quantités est obtenue par lissage des valeurs échantillonnées à l'aide
d'un noyau de lissage $W_h$. Ce noyau de lissage fournit la distribution massique sur une plage de
longueur $h$. Ainsi un champ scalaire continu $f: \mathbb{R}^3 \rightarrow \mathbb{R}$, dont on connaît seulement les valeurs
$f_i$ estimées aux points d'échantillonnage représentés par les particules, peut être approché par
la fonction suivante :

$$
\langle f(\mathbf{r}) \rangle = \sum_j m_j \frac{f_j}{\rho(\mathbf{r}_j)} W_h(\mathbf{r} - \mathbf{r}_j) \quad (3.5)
$$

La fonction $\rho: \mathbb{R}^3 \rightarrow \mathbb{R}$ représente la densité massique du milieu, $m_j$ est la masse de la particule $j$.

Plus encore, des valeurs échantillonnées $f_i$ du champ $f$, on peut déduire une expression
continue approchée de son gradient :

$$
\langle \nabla f(\mathbf{r}) \rangle = \sum_j m_j \frac{f_j}{\rho(\mathbf{r}_j)} \nabla W_h(\mathbf{r} - \mathbf{r}_j) \quad (3.6)
$$

Cette démarche, qui offre toute la souplesse du modèle à base de particules, permet de formuler
des équations d'état à partir de quantités continues. Le modèle SPH est donc plus proche de la
nature continue du milieu qu'il représente conduisant à des simulations de meilleure qualité.
Desbrun [Desbrun 1997] utilise ce modèle pour simuler des corps hautement déformables.

2.2 Mécanique des milieux continus

<figure><img src="image_32.png" /><figcaption>Figure 23 Représentation des configurations de référence et instantanée d'un milieu déformable.</figcaption></figure>

Cette deuxième section recense les méthodes qui s’appuient sur l’hypothèse de continuité de la matière constituant le corps déformable. Cette fois-ci le corps n’est plus représenté par un ensemble discret de particules $(\mathbf{x}^i(t))_{i \le i \le n}$. On définit une application continue $\mathbf{r}$ caractérisant la forme instantanée de l’objet. En fait, le terme $\mathbf{r}(\mathbf{r}^0, t)$ correspond à la position à l’instant $t$ d’une particule identifiée par sa position $\mathbf{r}^0$ dans une configuration de référence (Figure 23).
L'évolution du corps déformable est alors régie par l'équation suivante :

$$ \frac{\partial}{\partial t} \left( \mu(\mathbf{r}^0, t) \frac{\partial \mathbf{r}}{\partial t}(\mathbf{r}^0, t) \right) + \gamma(\mathbf{r}^0, t) \frac{\partial \mathbf{r}}{\partial t}(\mathbf{r}^0, t) + \frac{\partial \varepsilon}{\partial \mathbf{r}}(\mathbf{r}^0, t) = \mathbf{f}(\mathbf{r}(\mathbf{r}^0, t), t) \quad (3.7) $$

$\mu(\mathbf{r}^0, t)$ est la densité massique du corps déformable dans la configuration de référence, $\gamma$ un coefficient d'amortissement, $\varepsilon(\mathbf{r}^0, t)$ une énergie potentielle élastique responsable des efforts de cohésion interne et $\mathbf{f}(\mathbf{r}(\mathbf{r}^0, t), t)$ le vecteur résultant des efforts extérieurs.

### 2.2.1 Modèle de Terzopoulos

Terzopoulos et al. [Terzopoulos 1987] choisissent une formulation proche de l'équation (3.7) pour décrire l'évolution d'un corps déformable.

$$ \frac{\partial}{\partial t} \left( \mu(\mathbf{a}, t) \frac{\partial \mathbf{r}}{\partial t}(\mathbf{a}, t) \right) + \gamma(\mathbf{a}, t) \frac{\partial \mathbf{r}}{\partial t}(\mathbf{a}, t) + \frac{\partial \varepsilon}{\partial \mathbf{r}}(\mathbf{a}, t) = \mathbf{f}(\mathbf{r}(\mathbf{a}, t), t) \quad (3.8) $$

Dans leur modèle, le paramètre **a** de l'équation (3.8) ne représente pas obligatoirement la position dans l'espace de la particule pour une configuration de référence. En l'occurrence, ils considèrent des objets déformables définis à partir de courbes, surfaces ou volumes paramétriques. Ce choix est motivé par la volonté de pouvoir animer des objets issus d'une modélisation géométrique à base de surfaces de Bézier ou de surfaces NURBS. Ainsi le paramètre **a** correspond à l'indice paramétrique sur une courbe, surface ou volume identifiant une particule. Sous ces conditions, le paramètre **a** n'est pas rattaché à une grandeur physique. Afin de caractériser le changement de forme de l'objet, ils choisissent d'utiliser le **tenseur métrique** $G(\mathbf{a}, t)$.

$$ G(\mathbf{a}, t) = \frac{\partial \mathbf{r}}{\partial \mathbf{a}}(\mathbf{a}, t)^T \frac{\partial \mathbf{r}}{\partial \mathbf{a}}(\mathbf{a}, t) \qquad G^0(\mathbf{a}) = \frac{\partial \mathbf{r}^0}{\partial \mathbf{a}}(\mathbf{a})^T \frac{\partial \mathbf{r}^0}{\partial \mathbf{a}}(\mathbf{a}) \quad (3.9) $$

Ce tenseur permet de mesurer la dilatation de deux fibres matérielles à l'instant $t$ autour de la particule identifiée par le paramètre **a**, en fonction de deux directions paramétriques (Figure 24).

Intuitivement, on peut rapprocher l'expression du tenseur métrique $G(\mathbf{a}, t)$ avec celle du tenseur de Cauchy-Green $C(\mathbf{r}^0, t)$ de la théorie de l'élasticité (le lecteur trouvera la définition du tenseur de Cauchy-Green dans le chapitre consacré à l'élasticité non linéaire continue). Mais dans la mesure où le paramètre **a** ne correspond pas systématiquement à une position de l'espace, le tenseur métrique est dépourvu de toute grandeur physique.

La déformation locale instantanée du corps autour d'une particule identifiée par **a** est alors obtenue par la différence des tenseurs métriques $G(\mathbf{a}, t)$ et $G^0(\mathbf{a})$ correspondant respectivement à la configuration instantanée et à la configuration de référence. L'énergie totale de déformation est déduite par sommation :

$$ \mathcal{E}(t) = \int_{\Omega} \| G(\mathbf{a}, t) - G^0(\mathbf{a}) \|_{\alpha}^{2} d\mathbf{a} \quad (3.10) $$
Comme les tenseurs métriques ne correspondent pas à des grandeurs physiques, les coefficients de la norme matricielle $||\cdot||_α$ ne peuvent être identifiés à des quantités caractéristiques du matériau et doivent être par conséquent fixés de manière empirique.

<figure><img src="image_33.png" /><figcaption>Figure 24 Illustration du tenseur métrique, comparaison avec le tenseur de dilatation de Cauchy-Green.</figcaption></figure>

Au vu de l'expression complexe des équations à résoudre, Terzopoulos procède à plusieurs simplifications. Il effectue une troncature dans l'expression de la variation $\delta\epsilon/\partial\mathbf{r}$ de l'énergie de déformation pour ne conserver que les termes d'ordre 2. Puis, s'appuyant sur une grille discrète régulière, il débouche sur le système différentiel suivant :

$$ \mathbf{M} \frac{\partial^2 \mathbf{r}}{\partial t^2} + \mathbf{D} \frac{\partial \mathbf{r}}{\partial t} + \mathbf{K}(\mathbf{r}) \mathbf{r} = \mathbf{f} \qquad (3.11) $$

La matrice **M** est la **matrice des masses** (i.e. mass matrix), la matrice **D** est la **matrice d'amortissement** (i.e. damping matrix) et la matrice **K**(r) est la **matrice des raideurs** (i.e. stiffness matrix). Enfin, les vecteurs **r** et **f** correspondent respectivement aux valeurs des fonction r et f, échantillonnées aux points de la grille régulière. Il substitue alors l'ensemble des dérivées partielles par des différences finies. Bien qu'il considère des équations régissant le mouvement d'un milieu continu, Terzopoulos a recours à une discrétisation du milieu en une grille régulière pour la résolution de ces équations : résolution qui demeure coûteuse car elle nécessite à chaque itération l'inversion de la matrice des raideurs K(r).

En 1988, Terzopoulos et Fleischer [Terzopoulos 1988] proposent une modification du modèle initial en décomposant le mouvement d'un corps déformable en un mouvement rigide (translation et rotation) et une déformation. Ainsi, en considérant de petits déformations, il suppose la matrice des raideurs **K** constante, accélérant ainsi les temps de calcul des animations.
### 2.2.2 Méthode des éléments finis

La méthode des éléments finis est une méthode d'analyse numérique pour la recherche de solutions approchées d'équations d'évolution. Les équations d'évolution englobent à la fois les équations différentielles ordinaires (EDO) ainsi que les équations aux dérivées partielles (EDP). Cette méthode très populaire dans de nombreux domaines comme la mécanique, la biochimie mais aussi l'économie a fait l'objet de nombreux ouvrages dont ceux de Zienkiewicz [Zienkiewicz 1977], de Bathe [Bathe 1982], de Hughes [Hughes 1987]. Je conseille aussi au lecteur le cours très didactique de Garrigues [Garrigues 2001].

Je détaille le déroulement de la méthode des éléments finis pour le problème d'évolution suivant :

$$ \text{Trouver } u : \Omega \to \mathbb{R} \text{ telle que } \begin{cases} A(u(x)) = f(x) \quad \forall x \in \Omega \quad \text{(problème d'évolution)} \\ u(x) = g(x) \quad \forall x \in \partial\Omega \quad \text{(conditions aux limites)} \end{cases} \qquad (3.12) $$

A est un opérateur différentiel.

Le but de la méthode est de déterminer une solution approchée $\tilde{u}(x)$ de l'équation (3.12), définie par morceaux sur des sous-domaines $\Omega_i$ de l'ensemble $\Omega$. Chaque restriction $\tilde{u}|_{\Omega_i}(x)$ de la solution approchée $\tilde{u}(x)$ au sous-domaine $\Omega_i$ est en fait un polynôme interpolant les valeurs de la solution en des points caractéristiques du sous-domaine. De façon plus précise, la méthode des éléments finis se décompose en plusieurs étapes :

(1) Le **maillage** consiste à découper le domaine $\Omega$ en sous-domaines $\Omega_i$ appelés **mailles**. Au sein de chaque maille, on choisit un ensemble de **nœuds** correspondant aux points d'échantillonnage de la solution approchée de l'équation (3.12) et d'une fonction d'interpolation permettant d'exprimer les valeurs de la solution approchée en tout point de la maille $\Omega_i$. Une maille dotée à la fois de ses nœuds et de sa fonction d'interpolation est appelée **élément**.

(2) La **formulation variationnelle** consiste à transformer le problème d'évolution (3.12) en un problème intégral de la forme :

$$ \text{Trouver } u : \Omega \to \mathbb{R} \text{ telle que } \int A(u(x))\psi(x)dx = \int f(x)\psi(x)dx \quad \forall \psi \qquad (3.13) $$

En exploitant l'intégration par parties (ou la formule de Green-Stokes en dimension plus élevée), on peut procéder à des simplifications visant à diminuer l'ordre de dérivation de la solution du problème initial.

(3) La **discrétisation** permet, en choisissant un ensemble fini de $n$ fonctions $\psi$ particulières, de générer un nombre fini de $n$ équations linéaires dont les inconnues sont les valeurs de la solution approchée et de ses dérivées successives estimées en chaque nœud du maillage.

La méthode des éléments finis, réputée pour être coûteuse en temps de calculs, a été utilisée en synthèse d'images dans sa forme la plus simple. Les éléments considérés sont des tétraèdres, privilégiés pour la forme triangulaire de ses cotés qui constitue une forme de base pour la