MAO: CALCUL SCIENTIFIQUE

FILIPA CAETANO, ANTOINE LEVITT

VERSION PRÉLIMINAIRE, JUIN 2024

Table des matières

1. Programme 2
2. Introduction 2
2.1. Modélisation 3
2.2. Analyse numérique 3
2.3. Calcul scientifique 4
3. Arithmétique machine 4
4. Algèbre linéaire 5
4.1. Normes matricielles 5
4.2. Rayon spectral 5
4.3. Conditionnement et stabilité 6
5. Équations linéaires 7
5.1. Systèmes denses 7
5.2. La méthode de Richardson 8
5.3. Préconditionnement 9
5.4. Méthodes de Krylov 9
5.5. En pratique 9
TD/TP 9
6. Équations non-linéaires 10
6.1. Méthode de Newton 10
TD/TP 11
7. Optimisation 11
7.1. Optimisation sans contraintes 11
7.2. Optimisation sous contrainte 12
8. EDO 13
8.1. Méthode d’Euler 14
8.2. Commentaires et extensions 15
TD/TP 15
9. Introduction aux EDP 16
9.1. Modélisation 16
9.2. Classification 16
9.3. Solution par Fourier 17
TD/TP 17
10. EDP elliptiques 18
10.1. Caractère bien posé 18
10.2. Discrétisation 18
10.3. Analyse d’erreur 19
11. EDP d’évolution 20
11.1. Discrétisation en espace 21
1
 11.2. Discrétisation en temps 21

1. Programme
Le présent document forme des notes de cours volontairement concises. On se réfèrera aux ou-
vrages de référence pour plus de détails.
Équipe enseignante :
— Filipa Caetano
— Antoine Levitt
Contact : pré[email protected].
Syllabus (très prévisionnel...) :
— AL Lundi 29 avril : CM : introduction, normes matricielles, conditionnement
— AL Mardi 30 avril : TD/TP : arithmétique machine, algèbre linéaire
— AL Jeudi 2 mai : TD/TP : algèbre linéaire
— AL Lundi 6 mai : CM/TD : systèmes linéaires
— AL Mardi 7 mai : TD/TP : systèmes linéaires
— FC Lundi 13 mai : CM : équations non-linéaires et optimisation
— FC Mardi 14 mai : TD/TP : équations non-linéaires et optimisation
— FC Jeudi 16 mai : TD/TP : équations non-linéaires et optimisation
— FC Mardi 21 mai : EDO
— AL Jeudi 23 mai : Partiel
— FC Lundi 27 mai : EDO
— FC Mardi 28 mai : EDO
— FC Jeudi 30 mai : EDP
— FC Lundi 3 juin : Transport
— FC Mardi 4 juin : Transport
— FC Jeudi 6 juin : Transport
— AL Lundi 10 juin : EDP
— AL Mardi 11 juin : EDP
— AL Jeudi 13 juin : EDP
— AL Lundi 17 juin : EDP
— AL Mardi 18 juin : EDP
— AL Jeudi 20 juin : EDP
— FC Lundi 24 juin : examen
Bibliographie (disponible à l’agrégation)
— Programme de l’agrégation : Francis Filbet, Analyse numérique
— EDO : Jean-Pierre Demailly, Analyse numérique et équations différentielles
— Algèbre linéaire, optimisation, EDP : Grégoire Allaire, analyse numérique et optimisation
— EDP, approximation numérique : Brigitte Lucquin, Équations aux dérivées partielles et leurs
approximations

2. Introduction
La modélisation est la mise en équation d’un problème concret. L’analyse numérique est la
conception et l’étude mathématique d’algorithmes de résolution approchée d’équations. Le calcul
scientifique est l’art de la mise en oeuvre de ces algorithmes en pratique.
2.1. Modélisation. Quelles équations interviennent dans les questions suivantes ?
(1) Quelle sera la température moyenne de la terre dans 50 ans ?
2
 (2) À quelle vitesse un avion décolle-t-il ?
(3) Comment concevoir un avion indétectable par radar ?
(4) Quel est l’angle de la molécule d’eau H2 O ?
(5) Quel est le rendement d’une éolienne ?
(6) ... d’un panneau solaire ?
(7) Comment envoyer un satellite en orbite ?
(8) Comment dimensionner un pont ?
(9) Quelle est la consommation thermique d’un bâtiment ?
(10) Comment investir un porte-feuille d’actifs ?
(11) Quelle est la consommation énergétique d’un circuit électrique ?
(12) Au-delà de quel seuil de contagiosité une maladie devient-t-elle une épidémie ?
2.2. Analyse numérique. L’analyse numérique est l’art de concevoir des algorithmes de calcul
d’une approximation d’un objet bien défini, et de quantifier l’erreur commise. La plupart du temps,
ce sont des procédés
√ itératifs, et on n’obtient jamais la vraie solution. Ex : méthode babylonienne
de calcul de 2 (il y a ≈ 4 millénaires) : xn+1 = 12 (xn + x2n ). Même en arithmétique exacte et en
sachant calculer les racines, on ne peut pas calculer en un nombre fini les racines de polynômes de
degré ≥ 5 (donc par exemple les valeurs propres de matrices de taille ≥ 5).
Exemple de problèmes dont on peut vouloir trouver une solution :
(1) Système d’équations linéaire : Ax = b
(2) Système d’équations non-linéaires : f (x) = 0
(3) EDO : ẋ = f (x)
(4) EDP : ∂t u(x, t) = f (t, x, u, ∂x u, . . . )
(5) EDS : dXt = µ(Xt , t)dt + σ(Xt , t)dBt
(6) Optimisation continue : minimiser f (x), x ∈ RN sous contrainte g(x) = 0, h(x) ≥ 0
(7) Optimisation combinatoire : minimiser f (x), x à valeurs discrètes
(8) Interpolation, extrapolation, intégration de fonctions
Comment trouver des algorithmes ? Souvent, l’analyse donne une première réponse : une preuve
constructive = un algorithme. Vice versa, les algorithmes peuvent fournir des preuves.
Les preuves d’existence des théorèmes suivants sont-elles constructives ? Si oui, quelle est la
méthode numérique associée ?
(1) Existence de l’intégrale de Riemann d’une fonction continue
(2) Existence de l’intégrale de Lebesgue d’une fonction mesurable
(3) Bolzano-Weierstrass
(4) Théorème des valeurs intermédiaires
(5) Une fonction continue sur un compact admet un minimum
(6) Tout polynôme a une racine (complexe)
(7) Cauchy-Lipschitz
(8) x + Ax = b a une unique solution si ∥A∥op < 1
(9) Point fixe de Banach
(10) Point fixe de Brouwer
3
(11) Théorème des fonctions implicites
(12) Lax-Milgram
(13) Théorème d’approximation de Weierstrass
L’analyse numérique ne peut calculer que des quantités bien conditionnées, c’est-à-dire bien
définies et pas trop sensibles aux données. Exemple de problèmes mal conditionnés : Ax = b avec A
presque non inversible, déterminant, valeurs propres de matrices non symétriques, forme de Jordan,
comportement en temps long d’EDO chaotiques, équation de la chaleur en temps négatif...
Chercher à résoudre un problème mal conditionné, c’est poser la mauvaise question. Souvent,
la question finale est bien conditionnée, mais des quantités mal conditionnées apparaissent comme
intermédiaires de calcul : il faut reformuler le problème (eg pivot de Gauss plutôt que formules de
Cramer, forme de Schur plutôt que forme de Jordan...)

2.3. Calcul scientifique. Les points importants sont l’efficacité (en terme de temps d’ordinateur
mais aussi de temps de développement !), la maintenabilité, la robustesse... Les technologies utilisées
seront : langage de programmation Python, avec bibliothèques numpy (calculs sur tableaux), scipy
(méthodes numériques), matplotlib (visualisation) et les notebooks.
Quelques conseils de base de programmation :
— Tester chaque partie du code indépendamment, ne pas se précipiter sur l’écriture d’un code
complet
— Quand des bugs surviennent (inévitablement), ne pas relire vingt fois le code, mais tester
méthodiquement si toutes les propriétés qu’on pense être vraies sont effectivement vraies
— Éviter le copier/coller
— Séparez le code dans des fonctions, et ne pas hésiter à faire des fonctions qui prennent en
paramètre des fonctions
— Faire attention aux noms des variables. Utiliser des identificateurs descriptifs, par exemple for
it in range(Nt), for ix in range(Nx): plutôt que for i in range(N), for j in range(M)

3. Arithmétique machine
On rappelle qu’un ordinateur stocke les nombres réels sous la forme (modulo détails techniques...)
a2b , où a ∈ [1, 2[ est la mantisse, stockée via ses décimales en binaire, et b est l’exposant, stocké
comme un entier signé. Il en résulte que l’erreur entre un nombre x et sa représentation machine
fl(x) vérifie
|x − fl(x)| ≤ ε|x|,

pour x dans une plage très large de valeurs. ε est appelé l’epsilon machine : il est de l’ordre de 2−N ,
avec N le nombre de bits utilisés pour stocker a.
L’arithmétique sur ordinateur fonctionne comme si le résultat exact était calculé, puis arrondi au
nombre représentable le plus proche : calculer x+y sur un ordinateur calcule fl(fl(x)+fl(y)). Il résulte
des arrondis intermédiaires que l’arithmétique machine est commutative, mais non associative :
x + (y + z) fait un arrondi intermédiaire pour y + z, alors que (x + y) + z fait un arrondi pour x + y,
ce qui ne revient pas au même.
Le point important, et parfois contre-intuitif, est que les ordinateurs stockent les nombres avec
une précision relative uniforme : il y a beaucoup plus de nombres entre 0.001 et 0.002 qu’entre
100 et 100.001. En conséquence, la notion importante est celle de précision relative, et non absolue.
Diviser par un petit nombre par exemple n’est, en tant que tel, pas très gênant, même si on commet
ce faisant une grande erreur absolue : l’erreur relative reste de l’ordre de ε. En revanche, soustraire
deux nombres proches est dangereux, car l’erreur relative sur le résultat sera importante.
4
 √ 2
(1) Calculer 2 − 2 en python. Comment expliquer ce résultat ? Chercher sur internet sur com-
bien de bits est représentée la mantisse en python (qui utilise des nombres flottants double
précision par défaut), et vérifier la cohérence.
(2) Par essai et erreur, donner l’ordre de grandeur du plus grand nombre en virgule flottante
(attention, pas entier) représentable en python ? Cela est-il cohérent avec le nombre de bits
sur lequel est représenté l’exposant ?
(3) Calculer l’erreur commise par l’approximation
f (x + h) − f (x)
f ′ (x) ≈
h
par exemple sur f (x) = sin x en x = 1, et la tracer en fonction de h. Quelle valeur de h donne
le meilleur résultat ?
(4) Que peut-t-il
√ se passer si on calcule les racines de ax2 + bx + c = 0 par la formule x =
1 2 2
2a (b − b − 4ac) quand 4|ac| ≪ b ? Comment y remédier ?

4. Algèbre linéaire
L’algèbre linéaire joue un rôle prépondérant en calcul scientifique : la résolution d’équations
linéaires est un des rares ı̂lots de solvabilité dans une mer de problèmes insolubles sous forme
close. Beaucoup de problèmes (optimisation, résolution d’équations non-linéaires, etc) utilisent
une résolution de systèmes linéaires comme étape intermédiaire. Plus généralement, le formalisme
de l’algèbre linéaire fournit des outils puissants pour traiter tout type de phénomènes linéaires
(par exemple, équations différentielles linéaires), qui forment la base de notre compréhension des
systèmes même non-linéaires.
4.1. Normes matricielles. On rappelle que l’espace des matrices RM ×N forme un espace vecto-
riel. Plusieurs normes naturelles peuvent être placées sur cet espace : les plus utiles sont la norme
d’opérateur
p
∥Ax∥ ⟨x, AT Ax⟩
∥A∥op = max = max p .
x̸=0 ∥x∥ x̸=0 ⟨x, x⟩
qui voit une matrice A comme un opérateur de RM vers RN , et la norme de Frobenius
s X s X q
∥A∥F = 2
|Aij | = ∥Aen ∥ = Tr(AT A)
2

i=1,...,M,j=1,...,N i=1,...,N

qui voit une matrice de RM ×N comme un tableau de nombres, isomorphe à RM N . On vérifie que
∥A∥F est égale à la racine de la somme des valeurs propres de AT A (toutes réelles et positives, et
donc les racines sont appelées les valeurs singulières de A), alors que ∥A∥op est égale à la racine
de la plus grande. Quand on écrit ∥A∥ sans précisions, c’est souvent la norme d’opérateur qui est
sous-entendue. On vérifie facilement que ∥AB∥op ≤ ∥A∥op ∥B∥op
4.2. Rayon spectral. Une autre notion utile est le rayon spectral
ρ(A) = max |λ|,
λ∈σ(A)

où σ(A) est le spectre de A (ensemble de ses valeurs propres). Cette notion est particulièrement
utile pour étudier le système dynamique linéaire en temps discret
xn+1 = Axn
avec x0 fixé, qui est d’un grand intérêt dans l’étude de nombreux schémas numériques. On a
5
Théorème. limn→∞ An = 0 si et seulement si ρ(A) < 1.
Ce théorème est facile à démontrer si A est diagonalisable : on a avec les notations usuelles
An = P Dn P −1 , et (An → 0) ⇔ (Dn → 0) ⇔ (λn = 0 pour tout λ ∈ σ(A)). Quand A n’est pas
diagonalisable, on peut procéder par perturbation, ou utiliser la forme normale de Jordan. On a
même plus précisément le théorème de Gelfand
ρ(A) = lim ∥An ∥1/n ,
n→∞
cette limite étant bien définie et ne dépendant pas du choix de la norme.
Attention à ne pas confondre rayon spectral et norme matricielle, ainsi que valeurs propres et
valeurs singulières : les valeurs propres (indépendantes du choix de norme) déterminent le compor-
tement asymptotique des applications répétées de A, alors que les valeurs singulières (dépendantes
de la norme) déterminent l’expansion de la norme par une seule application de A. Un bon exemple
à garder en tête est la matrice
 
0 1
A= ,
0 0
de valeurs propres nulles (donc de rayon spectral nul) mais de valeurs singulières non-nulles. Cet
exemple montre également que ρ(A) n’est pas une norme matricielle.
Un autre bon exemple est
 
0.9 1000
(1) A=
0 0.9
La dynamique de An (exercice) montre une décroissance à l’infini comme 0.9n (comme prédit par
la formule de Gelfand), mais uniquement après une phase transitoire de grande amplitude.
Un cas particulier où tout est plus simple est celui des matrices symétriques (et plus généralement
normales, c’est-à-dire des matrices commutant avec leur transposée, qui se diagonalisent en base
orthonormée). Dans ce cas, les valeurs singulières sont les valeurs absolues des valeurs propres, et
on vérifie aisément en diagonalisant A que ρ(A) = ∥A∥op .
4.3. Conditionnement et stabilité. Imaginons que l’on souhaite résoudre Ax = b, mais que,
pour une raison ou une autre, notre membre de droite 1 est entaché d’une incertitude (ne serait-ce
que dûe à sa représentation numérique). Au lieu de résoudre Ax = b, on résoud en fait A(x + ∆x) =
(b + ∆b). Quelle est la taille de ∆x en fonction de celle de ∆b ? On a aisément
∆x = A−1 ∆b
et donc
∥∆x∥ ≤ ∥A−1 ∥∥∆b∥.
C’est une première réponse, qui donne une première indication sur la sensibilité aux erreurs : plus
A−1 est grande, plus x variera. Cette réponse n’est cependant pas très satisfaisante car elle n’est
pas invariante par échelle : le système 2Ax = 2b n’est pas intrinsèquement plus sensible aux erreurs
que Ax = b. Pour remédier à cela, on préfère travailler en erreur relative et écrire
∥∆x∥ ∥∆b∥
≤ ∥A∥∥A−1 ∥ .
∥x∥ ∥b∥
Le nombre ∥A∥∥A−1 ∥, appelé conditionnement κ(A) de la matrice A (égal au ratio de la plus
grande sur la plus petite valeur singulière) est adimensionnel, et quantifie la perte de précision
relative liée à la résolution du système linéaire. Bien sûr, ce n’est qu’une majoration dans le pire
cas, mais sa simplicité en fait quand même un concept utile.
1. On peut faire la même analyse pour une erreur sur A.
6
 Plus généralement, le conditionnement d’une opération (non nécessairement linéaire) quantifie
l’amplification d’erreur relative intrinsèque à l’opération. Par exemple, le calcul du déterminant
d’une
 matrice
 est mal conditionné en général, parce qu’une petite perturbation relative de la matrice
ε 0
pour ε petit produit une grande perturbation relative de son déterminant. C’est une
0 1/ε
notion à distinguer de celle de stabilité d’un algorithme qui réalise l’opération en question, et qui
réfère à sa capacité à obtenir des résultats aussi bons possibles (étant donné son conditionnement) en
présence d’erreur d’arrondi. L’élimination Gaussienne sans pivot est un algorithme potentiellement
instable, même si le problème lui-même est bien conditionné.
TD/TP
 
ε 0
(1) Calculer le conditionnement de et son déterminant. Même question sur 0.5IN ,
0 1/ε
l’identité en dimension N . Pourquoi le déterminant n’est-il pas une bonne mesure de l’inver-
sibilité d’un système ?
(2) Appliquer la méthode du pivot de Gauss à la main sur le système
    
ε 1 x a
=
1 1 y b
et en déduire que l’algorithme du pivot de Gauss sans échange est potentiellement instable.
(3) Déterminer la matrice correspondant à la discrétisation de l’équation −u′′ = f par différences
finies sur [0, 1] avec condition aux limites de Dirichlet. Montrer qu’elle est symétrique définie
positive. Calculer analytiquement ses valeurs et vecteurs propres ainsi que son conditionne-
ment.

5. Équations linéaires
On s’intéresse au système de N équations à N inconnues Ax = b. Si A est inversible, ce système
a une unique solution. Sinon, il possède soit une infinité de solutions (paramétrisées par le noyau de
A), soit aucune (alternative de Fredholm). La distinction se fait selon que b ∈ ran(A) = ker(AT )⊥
ou non (plus prosaı̈quement : si une combinaison linéaire de lignes donne 0 à gauche de l’équation,
on doit avoir 0 à droite). Dans la suite, on se place dans le cas où A est inversible et Ax = b admet
une unique solution x∗ .

5.1. Systèmes denses. Il existe des formules explicites pour l’inverse d’une matrice, en fonction
de déterminants de sous-matrices : les formules de Cramer. Ces formules sont algébriquement
plaisantes (par exemple, montrent facilement que l’inverse d’une matrice à coefficients rationnels
est à coefficients rationnels), mais analytiquement catastrophiques : les déterminants coûtent cher
à calculer, et sont numériquement instables.
La méthode standard pour résoudre Ax = b (numériquement comme à la main) est le pivot de
Gauss : on soustrait un multiple la première ligne aux autres de sorte à annuler le coefficient de la
première variable. On itère ainsi jusqu’à aboutir à une équation triviale pour la dernière variable,
et on remonte ensuite jusqu’à la première par substitutions successives. Techniquement, il faut
s’assurer que le coefficient de xn à la n-ième étape (le pivot) ne s’annule pas (ou n’est pas trop
petit) : cela peut être fait en échangeant des lignes ou des colonnes.
Cette méthode est souvent exprimée sous la forme d’une factorisation de matrices. Si on ignore
l’échange de lignes ou colonnes, l’algorithme amène A à avoir une forme triangulaire supérieure
par manipulations “triangulaires” successives de lignes (c’est-à-dire qu’à l’étape n, on ne modifie
que les lignes restantes). On peut donc écrire EN EN −1 . . . E1 A = U , où les En sont triangulaires
inférieures et U est triangulaire supérieure. Les matrices triangulaires inférieures étant stables par
7
multiplication et inversion, on peut aussi écrire A = LU avec L triangulaire inférieure. Cette forme
est pratique parce qu’elle permet de calculer facilement l’inverse
À l’étape n de l’algorithme, il faut modifier les N −n+1 éléments des N −n+1 lignes restantes : la
complexité totale est donc O(N 3 ). En pratique, on peut résoudre en quelques dizaines de secondes
sur un PC de bureau des systèmes de taille 10000. Calcul rapide d’ordre de grandeur : 1Ghz = 109
cycles par seconde ; un processeur à 4 Ghz de 8 coeurs qui effectue 4 opérations par cycle fait
donc 1012 opérations en à peu près 10s. Un modèle de climat typique utilise plusieurs millions
d’inconnues, et un pas de temps de l’ordre de l’heure : résoudre un système linéaire par cette
méthode à chaque pas d’une simulation de climat de 50 ans prendrait donc de l’ordre de 10000 ans
sur un ordinateur personnel : même sur des supercalculateurs, c’est impossible.

5.2. La méthode de Richardson. Très souvent, les matrices provenant des problèmes réels (par
exemple, équations aux dérivées partielles, modèles d’optimisation, etc) sont creuses : elles ne
contiennent que O(1) éléments non-nuls par ligne. L’élimination de Gauss ne permet typiquement
pas d’utiliser cette structure de façon optimale : lors des éliminations successives, la matrice se
remplit d’éléments non-nuls.
Une idée plus productive pour résoudre les très grands systèmes linéaires est d’utiliser une
méthode itérative, dont le prototype est la méthode de Richardson
(2) xn+1 = xn − α(Axn − b)
pour un pas α. La quantité Ax − b est le résidu de l’équation Ax − b, et quantifie à quel point on
est proche de la solution : en effet, on a la relation fondamentale erreur-résidu Ax − b = Ax − Ax∗ =
A(x − x∗ ), qui montre que ∥x − x∗ ∥ ≤ ∥A−1 ∥∥Ax − b∥.
L’intérêt de (2) est de n’avoir besoin que d’effectuer un produit matrice-vecteur à chaque itération.
Si A n’a que O(1) entrées non-nulles par ligne, ce produit peut être fait en O(N ) opérations. Clai-
rement, si l’itération converge, elle converge vers la solution du problème x∗ . Cette itération est
très naturelle : étant donné un x candidat, on ne peut pas calculer en général grand chose d’autre
que le résidu, et faire des combinaisons linéaires.
Pour étudier la convergence de la relation de récurrence inhomogène (2), la première étape est
de se ramener au cas homogène en posant en = xn − x∗ :
en+1 = (1 − αA)en
d’où en = (1 − αA)n e0 . L’itération convergera si les valeurs propres de 1 − αA sont plus petites
que 1 en module. C’est possible pour α suffisamment petit si toutes les valeurs propres de A sont
de partie réelle strictement positive (ou, en changeant le signe de α, négative). On peut préciser
cette convergence par exemple dans le cas où A est symétrique définie positive, de valeurs propes
0 < λ1 ≤ λ2 ≤ · · · ≤ λN . Dans ce cas, en projetant sur le i-ième mode propre vi , on obtient
⟨vi , en ⟩ = (1 − αλi )n ⟨vi , e0 ⟩.
On voit que pour optimiser la convergence, il faut choisir α de façon à ce que maxi=1,...,N |1 − αλi |
soit minimal. On voit aisément que le choix optimal est
2
αopt = ,
λ1 + λN
et que dans ce cas le taux de convergence est
2λ1 2
(3) r(αopt ) = 1 − =1−
λ1 + λN 1 + κ(A)
Plus le conditionnement κ(A) est grand, plus r(αopt ) est proche de 1, ce qui indique une convergence
lente.
8
5.3. Préconditionnement. On a vu que plus un problème est mal conditionné, plus il est difficile
à résoudre par la méthode de Richardson. C’est en fait parce que le résidu r = A(x − x∗ ) (ce
qu’on peut calculer) est une image déformée de l’erreur e = x − x∗ (ce qu’on voudrait calculer).
Comment peut-on travailler sur le résidu r pour l’amener à être plus proche de e ? Il arrive parfois
qu’on dispose d’une approximation P de la matrice A qui soit simple, au sens où on peut inverser
l’approximation aisément. On peut alors utiliser la méthode de Richardson avec re = P −1 r dans
l’espoir que re = P −1 Ae soit plus fidèle à e que r :
xn+1 = xn − αP −1 (Ax − b)
La même analyse que précédemment montre que si P −1 A a des valeurs propres proches de 1 (ou plus
généralement d’une constante), la méthode de Richardson préconditionnée convergera rapidement.
P est appelé un préconditionneur.
5.4. Méthodes de Krylov. La méthode de Richardson ci-dessus est intéressante pédagogiquement
car simple à analyser, mais elle est très sous-optimale, parce que c’est une méthode stationnaire
(la même relation, avec les mêmes coefficients, est appliquée à chaque itération) : elle n’est pra-
tiquement jamais utilisée en pratique. Des méthodes plus sophistiquées mélangent l’information
de plusieurs itérées précédentes pour obtenir une meilleure approximation : citons la méthode du
gradient conjugué (applicable aux matrices symétriques définies positives) et la méthode GMRES
(applicable aux matrices quelconques). Leur analyse fait appel à la théorie de l’approximation po-
lynomiale. Ces méthodes, combinées avec un bon préconditionneur, sont celles utilisées en pratique
pour la résolution de grands problèmes linéaires.
5.5. En pratique. En pratique, si la matrice est de petite taille (N ≤ 1000), on utilisera une
méthode directe de type pivot de gauss (numpy.linalg.solve). Si la matrice est grande et creuse
(structures de données de matrices creuses, scipy.sparse), on utilisera une méthode de pivot
de gauss adaptée (scipy.sparse.linalg.spsolve), ou une méthode de Krylov (par exemple
scipy.sparse.linalg.gmres).
TD/TP.
(1) On souhaite trouver une courbe qui reproduit un jeu de données (xn , yn )n=1,...,N . Pour cela,
on va se donner un ensemble de fonctions de base ϕp (x) (par exemple des polynomes) et
chercher une courbe ϕ(x) = Pp=1 αp ϕp (x) qui minimise la fonction objectif
P
 2
XN XP
E= yn − αp ϕp (xn )
n=1 p=1

Par exemple, si ϕ1 (x) = 1, ϕ2 (x) = x, on cherche à faire passer une droite.

(a) Reformuler le problème sous forme vectorielle en l’inconnue α ∈ Rp .
(b) Écrire les conditions d’optimalité
(c) Le coder dans le cas de la régression polynomiale (ϕp (x) = xp−1 ) et le tester sur des
données régulières et bruitées. Que se passe-t-il quand P est grand ?
(d) Refaire avec la fonction objectif
 2
XN P
X P
X
E=  yn − αp ϕp (xn )  +η αp2 .
n=1 p=1 p=1

(2) Montrer la formule de Neumann : si ∥A∥op < 1, alors 1 − A est inversible et (1 − A)−1 =
P∞ k
k=0 A . À quelle méthode numérique peut-on identifier cette formule ?
9
 (3) Écrire explicitement la méthode de Richardson préconditionnée dans le cas où α = 1 et
P = D, la diagonale de A. Montrer que la méthode se ramène à la méthode de Jacobi
Dxn+1 + (A − D)xn = b
(4) Montrer que si A est à diagonale strictement dominante, au sens où pour tout i = 1, . . . , N ,
X
|aii | > |aij |,
j̸=i

la méthode de Jacobi converge. On pourra montrer et utiliser le théorème de Gershgorin

P : les
valeurs propres d’une matrice A sont inclues dans l’union de disques ∪i=1,...,N B(aii , j̸=i |aij |)
(5) Combien d’itérations environ sont nécessaires pour résoudre un système linéaire de matrice
 
1 0
A=
0 100
à une précision proche de la précision machine avec la méthode de Richardson à pas optimal ?
Le vérifier numériquement.
(6) On considère la discrétisation de l’équation de Poisson −u′′ = f par différences finies sur [0, 1]
avec conditions aux limites de Dirichlet.
(a) Résoudre l’équation numériquement avec f (x) = sin(πx).
(b) Déterminer théoriquement et numériquement les vecteurs et valeurs propres de la matrice
A apparaissant dans la discrétisation par différences finies.
(c) Comment croı̂t le nombre d’itérations de la méthode de Richardson avec le nombre de
points de discrétisation ? Le vérifier en pratique.
(d) Que se passe-t-il si on augmente la taille du domaine, à pas de discrétisation fixé ? Même
question pour −u′′ + u = f .
(e) Résoudre analytiquement par fonction de Green l’équation −u′′ = f .
(f) Utiliser cette solution analytique pour construire un préconditionneur P de l’équation
−u′′ + ex u = f . Vérifier que la méthode préconditionnée converge plus rapidement, et
calculer numériquement les valeurs propres de P −1 A.

6. Équations non-linéaires
6.1. Méthode de Newton. Supposons qu’on veuille résoudre une équation non-linéaire f (x) = 0,
où f est une fonction régulière de RN vers RN (mais les arguments ci-dessous sont également valables
dans un espace de Banach). Supposons qu’on ait une estimation xn d’un zéro de f ; comment en
obtenir une meilleure ? On peut linéariser f autour de xn ,
f (x) ≈ f (xn ) + f ′ (xn )(x − xn ),
puis chercher un zéro de cette fonction affine. Cela mène directement à l’itération de Newton
xn+1 = xn − f ′ (xn )−1 f (xn ).
Cette itération converge localement vers des zéros non-dégénérés :
Théorème. Soit f : RN → RN une fonction de classe C 3 , et x∗ un zéro de f non-dégénéré, c’est-
à-dire que f ′ (x∗ ) est inversible. Alors il existe r > 0 tel que, pour tout x0 tel que ∥x0 − x∗ ∥ ≤ r,
l’itération de Newton démarrée de x0 converge vers x∗ . Cette convergence est quadratique, c’est-à-
dire qu’il existe C > 0 tel que
∥xn+1 − x∗ ∥ ≤ C∥xn − x∗ ∥2
10
Démonstration. Soit
G(x) = x − f ′ (x)−1 f (x)
l’application représentant une itération. On a G(x∗ ) = x∗ , et on calcule facilement que G est C 2
dans un voisinage de x∗ , avec G′ (x∗ ) = 0. Il existe donc un voisinage de x∗ et C > 0 tel que, pour
tout x dans ce voisinage,
(4) G(x) − x∗ = G(x) − G(x∗ ) ≤ ∥G′ (x∗ )(x − x∗ )∥ + C∥x − x∗ ∥2 = C∥x − x∗ ∥2 .
Cela montre que, en restreignant le voisinage à avoir un rayon r < 1/C, ce voisinage est stable par
G, et ∥xn+1 − x∗ ∥ ≤ rC∥xn − x∗ ∥ : la méthode de Newton converge. La convergence quadratique
suit par (4). □
En pratique, la méthode de Newton est souvent très coûteuse ; on l’approxime par des méthodes
approchant f ′ (x)−1 de façon plus ou moins explicite, comme les méthodes de Broyden. Le choix le
plus simple xn+1 = xn + f (xn ) est connu sous le nom de méthode de Picard, et converge localement
si 1 + f ′ (x) a des valeurs propres plus petites que 1 en module.
TD/TP.
(1) Montrer une variante du théorème de convergence de la méthode de Newton : si f est C 3
dans une boule de rayon δ autour de x0 et f ′ (x0 ) est inversible, alors il existe ε > 0 tel que,
si ∥f (x0 )∥ ≤ ε, la méthode de Newton converge vers un zéro de f . L’utiliser pour montrer le
théorème des fonctions inverses.
(2) Comment modifier la preuve ci-dessus si f n’est que de classe C 2 ? Que peut-on dire si f n’est
que de classe C 1 ?
(3) Tracer la fonction W de Lambert, définie par w(x)ew(x) = x
(4) Montrer que, pour une fonction f : C → C analytique, la méthode de Newton zn+1 =
zn − ff′(z n) 2
(zn ) est identique à celle qu’on obtiendrait si on considérait C comme R . Tracer
numériquement les bassins d’attraction (l’ensemble des conditions initiales pour lesquelles la
méthode de Newton converge vers chacune des racines) de l’équation z 3 = 1.
(5) Comment pourriez-vous utiliser le résultat de la question 1 pour démontrer mathématiquement
l’existence d’une solution à une équation non-linéaire ?

7. Optimisation
7.1. Optimisation sans contraintes. On considère le problème
min f (x)
x∈RN
avec f une fonction régulière. On rappelle qu’on a existence d’un minimum global si f est coercive
(est infinie à l’infini), et unicité si f est strictement convexe (il existe c > 0 tel que ∇2 f (x) ≥ c
pour tout x ∈ RN ). À un minimum x∗ , on a ∇f (x∗ ) = 0 et ∇2 f (x∗ ) ≥ 0.
Le problème de trouver le minimiseur global d’une fonction non-convexe est extrêmement dur,
et on doit souvent se contenter en pratique d’un minimiseur local, et d’une méthode itérative. Au
point x, dans quelle direction h aller pour diminuer f ? On a
f (x + h) = f (x) + ∇f (x) · h + O(∥h∥2 )
donc la direction de plus forte descente est h = −∇f (x). Cela suggère la méthode du gradient :
(5) xn+1 = xn − α∇f (xn )
avec un pas α à choisir.
11
Théorème. Si f est C 2 et vérifie c ≤ ∇2 f (x) ≤ C pour tout x ∈ RN pour deux constantes c, C > 0.
Alors il existe α0 tel que pour tout pas α ≤ α0 et point de départ x0 ∈ RN , l’itération (5) converge
vers l’unique minimum global de f .
Démonstration. On rappelle que c ≤ ∇2 f (x) implique que f est strictement convexe, donc coercive.
Il existe un unique minimiseur x∗ , unique solution de ∇f (x) = 0.
On a
f (xn+1 ) = f (xn ) − α∥∇f (xn )∥2 + O(Cα2 ∥∇f (xn )∥2 ).
En choisissant α0 suffisamment petit, on a
f (xn+1 ) ≤ f (xn ) − α2 ∥∇f (xn )∥2 ,
d’où il suit que f (xn ) a une limite, et que ∇f (xn ) → 0. Ce n’est pas encore suffisant pour obtenir
la convergence. Tous les points d’accumulation de xn vérifient ∇f (x) = 0, et sont donc égaux à
x∗ . Supposons maintenant que xn ne converge pas vers x∗ ; cela implique qu’il existe ε > 0 et
une sous-suite qui reste à une distance ε de x∗ . Par compacité des sous-ensembles de niveau de f
(conséquence de la coercivité), on peut en extraire encore une sous-suite qui converge ; comme la
limite ne peut être que x∗ , c’est une contradiction. □
Cette preuve utilise l’hypothèse très forte que f est strictement convexe. Dans le cas contraire,
mais quand f reste coercive, on a génériquement convergence, par exemple quand les points critiques
de f sont non dégénérés.
Comment choisir α en pratique ? Le choix d’un pas fixe est sous-optimal. Il existe une grande
variété d’algorithmes de recherche linéaire (c’est-à-dire de choix de α). Par exemple, un algorithme
consiste à l’étape n à choisir un pas initial, puis à le réduire tant que la condition d’Armijo (de
“descente suffisante”)
f (xn − α∇f (xn )) ≤ f (xn ) − αβ∥∇f (xn )∥2
n’est pas satisfaite, pour une valeur de β ∈]0, 1[ fixée.
Quelle est la vitesse de convergence de la méthode de gradient ? Supposons qu’on démarre (ou
qu’on arrive après plusieurs pas de la méthode) au voisinage d’un minimum local non-dégénéré.
Alors on peut linéariser l’équation en calculant la jacobienne de l’application d’itération x 7→
f (x) − α∇f (x), qui vaut 1 − α∇2 f (x). On retrouve alors exactement l’analyse de la méthode de
Richardson, avec A = ∇2 f (x). La conclusion est, comme pour la méthode de Richardson, que le
taux de convergence se dégrade avec l’augmentation du conditionnement de A.  
a 0
Que veut dire cette convergence en pratique ? On voit sur l’exemple de A = (faire un
0 b
dessin, en représentant les lignes de niveau de f ) que la méthode converge en zigzag lorsque a et b
sont différents. Pour corriger ce zigzag, trois solutions :
— Le préconditionnement xn+1 = xn − αP −1 f (xn ) ;
— L’utilisation de la méthode de Newton xn+1 = xn − ∇2 f (xn )−1 ∇f (xn ), très efficace locale-
ment, mais coûteuse et convergeant également vers les points critiques non minima ;
— L’utilisation de méthodes de quasi-Newton visant à approximer ∇2 f (x)−1 sans le calculer,
comme le gradient conjugué ou l’algorithme BFGS.

7.2. Optimisation sous contrainte. On rappelle dans un cadre simple la méthode des multi-
plicateurs de Lagrange. Supposons par simplicité qu’on veuille optimiser f (x) sous la contrainte
g(x) = 0. Supposons également que g soit C 2 et ∇g(x) ̸= 0 dès que g(x) = 0, ce qui signifie que la
contrainte g(x) définit localement une sous-variété de dimension N − 1. Le plan tangent en x est
défini par l’équation ∇g(x)T h = 0.
12
 Si on est à un minimum x∗ du problème, on doit avoir f (x) ≥ f (x∗ ) pour tout y proche de
x tel que g(x) = 0. Soit h dans le plan tangent en x∗ . En appliquant le théorème des fonctions
implicites à l’équation g(x∗ + th + u∇g(x∗ )) = 0 en u à paramètre t, on obtient une courbe u(t)
qui est C 1 et qui vérifie u(0) = 0, u′ (0) = 0. Cela donne un chemin C 1 sur la variété x(t) =
x∗ + th + u(t)h⊥ tel que x(0) = x∗ et x′ (0) = h. Il suit que ∇f (x)T h = 0 pour tous ces h. On
a donc que ∇f (x∗ ) est orthogonal à l’orthogonal de ∇g(x∗ ), c’est-à-dire qu’il existe λ ∈ R tel
que ∇f (x∗ ) = λ∇g(x∗ ). Ce λ est appelé multiplicateur de Lagrange. Géométriquement, il traduit
simplement le fait que le gradient est orthogonal au plan tangent (donc appartient à l’espace normal,
c’est-à-dire Im(∇g(x∗ ))).
Pour résoudre un problème d’optimisation sous contrainte, on peut s’attaquer aux conditions
d’optimalité ∇f (x) = λ∇g(x), g(x) = 0 vu comme un système non-linéaire de N + 1 équations en
N +1 variables, et utiliser par exemple la méthode de Newton, mais on n’obtient que la convergence
locale, et vers tous les points critiques, pas nécessairement minima. On peut également utiliser
l’algorithme du gradient projeté : on se dirige dans la direction du gradient projeté sur le plan
tangent, puis on se ramène sur la variété des contraintes
(1) Donner la solution du problème min x2 + y 2 sous contrainte 2x + y = 1
(2) Calculer le gradient et la Hessienne de x 7→ 21 xT Ax − xT b avec A symétrique définie positive
(3) Montrer le théorème d’approximation par des matrices de rang 1 (cas particulier du théorème
d’Eckart-Young) : pour toute matrice A ∈ RN ×M ,
min ∥A − uv T ∥F = σ2 (A)
u∈RN ,v∈RM

où σ2 (A) est la deuxième plus grande valeur singulière de A, avec le minimum atteint lorsque
u et v sont les vecteurs singuliers à gauche et à droite associés à la plus grande valeur singulière
(c’est-à-dire que AT Au = σ1 (A)2 u, AAT v = σ1 (A)2 v).
(4) Appliquer la méthode du gradient à la fonction 12 xT Ax − xT b avec A symétrique définie
positive, et montrer qu’elle s’identifie à l’itération de Richardson
 
10 0
(5) Appliquer la méthode du gradient à la fonction ci-dessus pour A = , b = 0, et tracer
0 1
les itérés
(6) Donner l’algorithme du gradient projeté pour une fonction f (x) soumise à la contrainte
∥x∥2 = 1
(7) Implémenter la méthode du gradient projeté pour le problème de séparation de source (voir
code fourni) : on se donne d sources de longueur N stockée dans une matrice X ∈ RN ×d
(centrée et normalisée), on cherche un vecteur de poids w ∈ Rd avec ∥w∥ = 1 qui minimise
f (x) = E(G(Xw))
avec G(x) = log(cosh(x)) (qui mesure la déviation par rapport à la gaussianité).

8. EDO
On considère l’équation
ẋ = f (x)
avec f : RN → RN pour simplicité. La théorie s’étend au cas non-autonome (f (x, t)), au cas où x
vit dans un espace de Banach, et aux équations d’ordre plus élevé (par exemple via le changement
de variable yn = ẋn ).
On rappelle le théorème de Cauchy-Lipschitz : si f est localement Lipschitz (par exemple, C 1 )
en x0 , on a existence et unicité pour un temps fini. Si f est globalement Lipschitz (ou qu’elle est
13
localement Lipschitz et qu’on arrive à montrer que la solution reste bornée), la solution existe pour
tout temps.
8.1. Méthode d’Euler. La méthode d’Euler est
xn+1 = xn + ∆tf (xn )
Théorème. Soit f : RN → RN uniformément C 1 (i.e. f est globalement C 1 et sa différentielle est
uniformément bornée), x0 ∈ RN , (x(t))t∈R l’unique solution de ẋ = f (x), x(0) = x0 , et (xn )n∈N la
suite définie par
xn+1 = xn + ∆tf (xn ).
Alors, pour tout T ∈ R, il existe CT > 0 tel que, pour tout N ∈ N, si ∆t = T /N , alors
sup |xn − x(n∆t)| ≤ CT ∆t
n=0,...,N

Démonstration. On définit l’erreur de troncation locale (erreur de consistance, par mauvaise tra-
duction de l’anglais “consistent”, qui signifie “cohérent”)
tn = x(tn+1 ) − (x(tn ) + ∆tf (x(tn ))),
à différencier de l’erreur
en = xn − x(tn )
(faire le dessin fondamental : x(tn ), x(tn+1 ), xn , xn+1 , tn , en ). Avec un développement limité à l’ordre
2 de x(t), on montre aisément que |tn | ≤ C∆t2 avec C > 0 ne dépendant pas de ∆t.
On a la relation de récurrence
en+1 = xn+1 − x(tn+1 )
= (xn + ∆tf (xn )) − (x(tn ) + ∆tf (x(tn ))) − tn
|en+1 | ≤ (1 + ∆tL)|en | + |tn |
où L est la constante de Lipschitz de f , qui est une estimée de stabilité.
On résoud maintenant l’inégalité récurrente (“lemme de Gronwall discret”)
an+1 ≤ αan + β
en montrant que an ne peut pas croı̂tre plus vite que dans le cas d’égalité. Soit bn+1 = αbn + β. On
pose la différence dn = an − bn , on a dn+1 ≤ αdn donc dn ≤ 0.
Résolvons maintenant la récurrence inhomogène bn+1 = αbn + β. On utilise la méthode de
variation de la constante : bn = cn αn . On remplace :
cn+1 = cn + βα−(n+1)
d’où
n
X
cn = b0 + βα−k
k=1
n
X
bn = αn b0 + β αn−k
k=1
n−1
′
X
= αn b0 + β αk
k=0
αn −
1
= αn b0 + β
α−1
14
Avec α = (1 + ∆tL), β = C∆t2 , on obtient
|en | ≤ C∆t((1 + ∆tL)n − 1) ≤ C∆teLtn
□
8.2. Commentaires et extensions.
— D’un point de vue mathématique, la méthode d’Euler peut être utilisée comme preuve alterna-
tive d’existence de solution : schématiquement, on interpole linéairement les points xn pour
construire une solution continue approchée x eN (t), et on utilise une méthode de compacité
(théorème d’Arzela Ascoli) pour extraire une sous-suite convergente.
— Consistance (on résoud l’équation de façon approchée à chaque étape) + stabilité (ampli-
fication d’erreur raisonnable) = convergence
— La notion de stabilité est la même que celle utilisée pour prouver la stabilité d’une EDO par
rapport à des erreurs sur x0 et/ou f .
— Le schéma de preuve marche aussi pour d’autres méthodes, comme celle ci-dessus : on a
besoin de borner l’erreur de consistance, et que la fonction d’incrément xn → xn+1 soit
Lipschitzienne de constante 1 + ∆tL′ . Par exemple, pour une méthode de type Runge-Kutta
à l’ordre 2
∆t
xn+1 = xn + (f (xn ) + f (xn + ∆tf (xn ))),
2
en refaisant le même raisonnement, on obtiendrait une erreur totale d’ordre ∆t2 .
— L’hypothèse d’être globalement Lipschitz est très forte et permet d’obtenir facilement l’exis-
tence d’une solution globale en temps. On n’utilise en fait dans la preuve que le caractère
Lipschitzien autour de la solution.
— En pratique, eLt peut être énorme. La méthode d’Euler est imprécise en temps long, mais la
dynamique sous-jacente potentiellement aussi... systèmes chaotiques, effet papillon.
— La stabilité de l’EDO n’implique pas la stabilité du schéma numérique. Par exemple, le schéma
multi-pas
xn+2 + 4xn+1 − 5xn = ∆tf (xn )
est consistant, mais instable : si on prend f = 0, on voit que la solution est donnée par
α + β(−5)n . Si la condition initiale (x0 , x−1 ) est entachée d’une erreur d’ordre h, β sera
d’ordre h, et l’erreur h(−5)n ne tend pas vers 0 pour tout n = 1, . . . , T /h.
— En pratique :
(1) Ne pas utiliser la méthode d’Euler, utiliser une méthode d’ordre plus élevé
(2) Utiliser une méthode à pas variable
(3) Ne pas coder à la main (risque d’erreur)
— Notion de raideur (stiffness) : ẋ = −λx avec λ > 0, solution explicite avec la méthode
d’Euler : besoin de ∆t ≤ λ2 pour ne pas avoir d’explosion. Extension à ẋ = Ax. Problèmes
pratiques avec les oscillations/amortissements rapides (eg EDP). Méthode d’Euler implicite
xn+1 = xn + ∆tf (xn+1 ). Difficultés d’implémentation.
TD/TP.
(1) Refaire l’analyse de la méthode d’Euler à la main quand f (x) = Ax
(2) Refaire l’analyse pour la méthode de Heun
x
en+1 = xn + ∆tf (xn )
f (xn ) + f (e
xn+1 )
xn+1 = xn + ∆t
2
15
 (3) Coder la méthode d’Euler et de Heun sur l’oscillateur harmonique ẍ + ω 2 x = 0, valider son
ordre.
(4) Interpréter et simuler les équations proie-prédateur de Lotka-Volterra
ẋ = x(α − βy)
ẏ = y(δx − γ)
Calculer le point fixe non trivial, estimer la période des solutions par linéarisation, et l’utiliser
pour vérifier votre implémentation.
(5) Faire la même chose avec l’équation d’un pendule ẍ + ω 2 sin x = 0. Coder la méthode d’Euler.
Comparer son efficacité (erreur par rapport au nombre d’évaluation du membre de droite)
par rapport à scipy.integrate.solve ivp. Tracer la courbe de la période en fonction de
l’amplitude initiale (comment ?).

9. Introduction aux EDP

9.1. Modélisation. Dérivation physique des équations de transport, de la chaleur, des ondes.
Équation de la chaleur et mouvement brownien. Gravitation/électrostatique et équation de Poisson.

9.2. Classification. Au contraire des EDO, les EDP ne possèdent pas de théorie systématique.
On se limite aux EDP linéaires à coefficients constants.
Ordre 1 :
a∂x u + b∂y u = 0
Adimensionalisation. Modélisation : transport. Solution explicite. Condition initiale/aux limites.
Ordre 2 :
a∂xx u + 2b∂xy u + c∂yy u = 0
 
T a b
Formellement, ∇ A∇u = 0 avec A = . On peut diagonaliser A et obtenir ses valeurs propres
b c
λ (≥ 0 sans perte de généralité) et µ réels, ce qui permet de classifier, par analogie avec les courbes
quadriques, en
— λ > 0, µ > 0 : elliptique (type Laplace/Poisson)
∂xx u + ∂yy u = 0
— λ > 0, µ < 0 : hyperbolique (type ondes)
∂yy u = ∂xx u
— λ = 0 : parabolique (type chaleur)
∂y u = ∂xx u
Dans les cas des équations de transport, d’ondes et de la chaleur, on voit y comme une variable
temporelle. C’est le type des équations qui dicte les conditions initiales/aux limites qu’il faut mettre
pour obtenir une unique solution ; il n’est par exemple pas naturel de considérer y comme une
variable temporelle d’une équation elliptique, car il faudrait dans ce cas donner une condition
initiale et une condition finale.
On peut combiner ces types d’équations. Par exemple, réaction-advection-diffusion :
∂t u = α∂xx u + β∂x u + f (u)
16
9.3. Solution par Fourier. Une méthode très puissante pour la résolution d’EDP linéaires à coef-
ficients constants est l’utilisation des séries (dans le cas d’un intervalle borné) ou de la transformée
(dans le cas de l’espace entier) de Fourier. Dans le cas d’un intervalle borné, il faut que les conditions
au bord s’y prêtent, et on considérera ici des conditions aux limites périodiques. D’autres transfor-
mations dérivées de Fourier peuvent être pertinentes pour d’autres types de conditions aux limites.
Pour des domaines multidimensionnels plus compliqués, les exponentielles complexes doivent être
remplacés par les modes propres, ce qui sort du cadre de ces notes.
Considérons par exemple l’équation de la chaleur
∂t u = ∂xx u, u(0, t) = u(2π, t), ∂x u(0, t) = ∂x u(2π, t), u(x, t = 0) = u0 (x).
On suppose par simplicité que la condition initiale u0 est C ∞ et 2π-périodique. On a alors
1 X
u0 (x) = √ cn (0)einx ,
2π n∈Z
et cn (0) décroı̂t plus vite que tout polynome inverse en |n|. Il est naturel de chercher u(x, t) sous
la forme
1 X
u(x, t) = √ cn (t)einx .
2π n∈Z

À condition que cn décroisse suffisamment vite (ce que nous vérifions a posteriori ), on peut dériver
sous le signe somme et obtenir
ċn (t) = −n2 cn (0)
2
donc cn (t) = e−n t cn (0), et
1 X 2
u(x, t) = √ cn (0)e−n t einx .
2π n∈Z
ce qui fournit une solution de l’EDP. On voit que les coefficients de Fourier sont atténués, d’autant
plus qu’ils sont de fréquence élevée. En temps long, u tend vers la fonction constante (la température
s’homogénéise).
Cette résolution explicite est permise par le fait que chaque mode propre einx soit une fonction
propre du second membre et donc évolue de façon indépendante. C’est une propriété générale des
équations linéaires invariantes par translation.
On pourrait aussi utiliser le théorème de la convolution pour reformuler u(x, t) comme la convo-
lution de u0 et d’une certaine
√ fonction, la fonction de Green (ici, une gaussienne périodisée de
largeur proportionnelle à t).
TD/TP.
(1) Faire la même étude avec l’équation de transport.
(2) De même avec l’équation de Schrödinger i∂t ψ = ∂x2 ψ. Parmi les équations ci-dessus, de quelle
équation est-elle la plus proche ?
(3) Résoudre numériquement l’équation de la chaleur sur [0, 2π] avec conditions périodiques et
condition initiale ecos x
(a) Montrer que la formule de quadrature
Z 2π N −1
2π X
f 2π Nk


f (x)dx ≈
0 N
k=0

est exacte quand f (x) = einx avec |n| < N .

17
 (b) Estimer la décroissance des coefficients de Fourier de ecos x .
(c) Déduire de ces deux points un ordre de grandeur du nombre de points nécessaires pour
obtenir une approximation numérique précise à 10−6 près du n-ième coefficient de Fourier
de ecos x . Comparer cet ordre de grandeur à celui donné par l’estimation naı̈ve en O(1/N )
de l’erreur de quadrature de la méthode de Riemann.
(d) Résoudre numériquement l’équation de la chaleur avec u0 (x) = ecos x . Visualiser l’évolution
de la solution et de ses coefficients de Fourier au cours du temps. Vous pouvez utiliser les
fonctions matplotlib pcolormesh (pour visualiser u(x, t) en couleur), plot surface (pour
visualiser u(x, t) en 3D), et FuncAnimation (pour animer u(x) au cours du temps).
(e) Que se passe-t-il si vous essayez de résoudre l’équation de la chaleur en temps négatif ?
(f) Que se passe-t-il si vous ajoutez un terme en ∂x u à l’EDP ?

10. EDP elliptiques

L’EDP elliptique la plus simple est
(6) − u′′ (x) = f (x)
(7) u(0) = u(1) = 0
avec f continue (par exemple). C’est en fait une équation différentielle (une seule variable), mais
de type “boundary value problem” (BVP) où les valeurs au bord sont données, plutôt que “initial
value problem” (IVP, problème de Cauchy), où la condition initiale est donnée. Cela la rend plus
proche d’une EDP que d’une EDO.
10.1. Caractère bien posé. Dans ce cas très précis, l’équation est facile à résoudre explicitement
(au sens où elle se ramène à un calcul de primitive) : on trouve une primitive seconde u e(x) de f (x),
et la solution générale est donc u(x) = u e(x) + a + bx. Les conditions u(0) = u(1) = 0 donnent un
système de deux équations linéaires indépendantes pour a et b, qui peuvent donc être résolues. On
obtient immédiatement par cette méthode que u est C 2 , et même, si f est C k , C k+2 (en dérivant
l’équation k fois).
Noter que cette méthode est très spécifique à cette équation, et ne peut pas par exemple traiter
d’équation de type −u′′ +u = f ou −∆u = f en plusieurs dimensions. Dans un cadre plus général, on
peut utiliser des méthodes variationnelles. L’unicité est facile à obtenir par un argument d’énergie :
si u et v sont deux solutions, alors w = u − v satisfait
(8) −w′′ = 0, w(0) = w(1) = 0.
On multiplie par w et on intègre, ce qui donne w′2 = 0, donc w′ = 0, ce qui combiné avec les
R
conditions
R ′2 auR bord donne w = 0. L’existence est plus subtile et peut être obtenue par minimisation
1 1
de 2 u − uf dans l’espace H0 (0, 1) (ou le théorème de Lax-Milgram, une généralisation de
cette approche), ce qui sort du cadre de ces notes.
10.2. Discrétisation. Cette équation peut être discrétisée de nombreuses manières ; on considère
ici les différences finies, qui consistent à poser un maillage xn = n/N de pas h = 1/N et à chercher
UN = (un )n=1,...,N tel que
−un+1 − un−1 + 2un
(9) = f (xn ),
h2
avec la convention u0 = uN = 1. Cette méthode se reformule également en
(10) AN U N = F N
avec (FN )n = f (xn ) et AN une matrice tridiagonale avec −1/h2 sur la sous et sur diagonale et 2/h2
sur la diagonale.
18
10.3. Analyse d’erreur. L’objectif est de montrer une borne sur l’erreur
(11) EN = UN − PN (u),
où PN est l’opérateur d’évaluation sur la grille : (PN (u))n = u(xn ).
La seule chose qu’on connaisse sur UN est AN UN , donc on écrit
(12) AN EN = AN UN − AN PN (u)
(13) = FN − AN PN (u)
C’est l’équation centrale : elle énonce que l’erreur EN vérifie une équation du même type que UN ,
mais avec en membre de droite le résidu de la solution exacte dans le schéma numérique (l’analogue
de ce que l’on avait appelé l’erreur de troncation locale pour les équations différentielles ordinaires).
Nous avons déjà croisé ce type de relation “erreur-résidu”, qui jouent un rôle fondamental en analyse
numérique.
On va procéder par la méthode de consistance et stabilité : on montre que le résidu est petit
(consistance) et que cela implique que EN est petit (stabilité).
Consistance. Le résidu est relativement aisé à borner : si on suppose par exemple que f est C 2 ,
alors u est C 2 , et on a
−u(xn+1 ) − u(xn−1 ) + 2u(xn )
(14) f (xn ) = −u′′ (xn ) = + O(h2 ),
h2
où le O est uniforme en n : explicitement, il existe une constante de consistance C1 indépendant
de N tel que ∥FN − AN PN (u)∥∞ ≤ C1 h2 .
Stabilité. On a vu par diagonalisation explicite de AN sur les modes sin( πkn N ) que sa plus petite
1 π
valeur propre h2 (2 − 2 cos( N ))est bornée inférieurement, indépendamment de N . Comme AN est
symétrique, on a pour tout V ∈ RN −1
N −1
X 1
(15) A−1
N V = ⟨Pk , V ⟩Pk
λk
k=1
avec Pk une base orthonormée de vecteurs propres de AN associés aux valeurs propres λk . Cela
implique par Pythagore que
v
uN
−1
uX 1
(16) ∥AN V ∥ = t 1
⟨p , V ⟩2 ≤ min
λ2k k
∥V ∥2 ≤ C2 ∥V ∥2
k |λk |
k=1
avec C2 une constante de stabilité ne dépendant pas de N .
Convergence. On a enfin
√ √
(17) ∥EN ∥ ≤ ∥AN ∥op ∥FN − AN U ∥ ≤ C2 N ∥FN − AN U ∥∞ ≤ C1 C2 N h2
soit, explicitement,
v
u
u1 XN
(18) t |un − u(xn )|2 ≤ C1 C2 h2 .
N
n=1

C’est une estimation d’erreur en norme L2 discrète.

Commentaires :
— L’estimation d’erreur est asymptotiquement optimale (on ne s’attend pas à un meilleur ordre
en général)
√ ne donne pas une estimation asymptotiquement optimale de |un − u(xn )| ponctuellement
— Elle
( N h2 au lieu de h2 ) : pour cela, il faut travailler en norme ∞, donc borner la norme de A−1
N
dans cette norme : cela peut se faire avec un principe du maximum discret
19
 — La notion de stabilité du schéma est que, si AN UN = FN et AN U eN = FeN avec FN et FeN
proches, alors UN et UeN sont proches. Cette stabilité existe également au niveau continu.
— Faire l’analyse numérique nécessite de 1) établir que l’équation initiale est bien posée, et
en particulier que sa solution dépend de façon continue des données (stabilité continue) 2)
montrer que la discrétisation respecte suffisamment la structure du problème pour qu’on
puisse prouver la stabilité au niveau discret 3) montrer la consistance.
TD/TP
(1) Refaire la même analyse pour des conditions au bord périodiques u(0) = u(1), u′ (0) = u′ (1).
(2) On considère un schéma général
αu(xn ) + β(u(xn+1 ) + u(xn−1 )) + γ(u(xn+2 ) + u(xn−2 ))
u′′ (xn ) ≈ .
h2
— Montrer que le schéma est consistent si α + β + γ = 0, α + β + 4γ = 2
— Montrer que, si on utilise ce schéma pour discrétiser l’équation −u′′ = f sur [0, 1] avec
conditions au bord périodiques, alors les valeurs propres de AN sont données par
α + 2β cos(k) + 2γ cos(2k)
λ(k) =
h2
2πZ
avec k ∈ N
— Montrer qu’il existe des valeurs de α, β, γ pour lesquelles λ(k) n’est pas nécessairement
positive sur R. En déduire que le schéma de différences finies peut être instable.
(3) Refaire la même analyse pour −u′′ + u = f avec f de classe C 2 , en admettant que cette
équation admette une unique solution de classe C 4 . Même question pour −u′′ + g(x)u =
f , avec g(x) positive, de classe C ∞ et régulière. Même question pour −u′′ + u′ = f , en
discrétisant u′ par des différences finies centrées u′ (x) ≈ (u(x+h)−u(x−h))/(2h). Pour cette
dernière question, on pourra montrer et utiliser le théorème de Lax-Milgram discret : si S est
symétrique définie positive avec S ≥ α au sens des matrices symétriques et A antisymétrique,
alors S + A est inversible et ∥(S + A)−1 ∥op ≤ 1/α
(4) Montrer l’unicité de solutions de l’équation −u′′ + u = f avec conditions de bord de Neumann
u′ (0) = u′ (1) = 0. Proposer une méthode de résolution numérique et montrer la convergence.
Que dire de l’équation −u′′ = f avec ces conditions de bord ?
(5) Résoudre analytiquement et numériquement l’équation −u′′ = 1 avec conditions au bord de
Dirichlet sur [0, 1]. Calculer et tracer quelques modes propres du Laplacien. Mêmes questions
pour −∆u = 1 sur [0, 1]2 . Pour simplifier la manipulation d’indices, on pourra utiliser des
tableaux quadri-dimensionnels et la fonction reshape.
(6) Empiriquement, comment se comporte la solution de −∆u + αχD (x)u = f quand α devient
grand, si χD est l’indicatrice d’un sous-ensemble du domaine de résolution de l’équation ? En
utilisant cette méthode, calculer des approximations des fonctions propres du Laplacien sur
un domaine en forme de L.
(7) Résoudre numériquement l’équation
(19) −εu′′ + u′ = 1, u(0) = 0, u(1) = 0.
Que constatez-vous ? Comment faut-il choisir ∆x par rapport à ε ?

11. EDP d’évolution

On considère une équation aux dérivées partielles d’évolution, dont le prototype est l’équation
de la chaleur
∂u ∂2u
(20) =
∂t ∂x2
20
toujours sur [0, 1], avec des conditions au bord de Dirichlet et une condition initiale u0 régulière
(par exemple C ∞ ). Cette équation est bien posée et la solution est régulière : cela peut se voir
ici avec une décomposition en série de Fourier (sur les modes sin(nπx)). Dans le cadre général, il
existe beaucoup de méthodes pour montrer le caractère bien posé (calcul fonctionnel des opérateurs
auto-adjoints, théorie des semi-groupes). On a notamment
Z Z
d
(21) ∥u∥L2 (Ω) = 2 u∆u = −2 |∇u|2 ≤ 0.
2
dt Ω Ω
Il s’ensuit que la norme L2 de u ne peut que décroı̂tre au cours du temps, ce qui traduit le caractère
dissipatif de l’équation (et donc, l’équation étant linéaire, une certaine forme de stabilité : si u0 et
v0 sont proches, alors u(t) et v(t) le sont).
Cette équation nécessite une discrétisation en espace et en temps. On commence par étudier
l’équation semidiscrète
dUN
(22) = AN UN , (UN (0))n = u0 (xn )
dt
avec les mêmes notations que précédemment.
11.1. Discrétisation en espace. Comme précédemment, on écrit une équation pour l’erreur
EN (t) = UN (t) − PN (u(t)) :
dEN
(23) = AN UN − PN (u′′ )
dt
(24) = AN EN + AN PN (u) − PN (u′′ )
| {z }
rN (t)

C’est encore une relation erreur-résidu : cette fois, le résidu apparaı̂t comme terme source de
l’équation. On utilise le principe de Duhamel (variation de la constante) pour obtenir, comme
EN (0) = 0,
Z t
′
(25) EN (t) = e−AN (t−t ) rN (t′ )dt′
0
Comme AN est positive définie, ∥e−AN t ∥op ≤ 1 (le montrer par une diagonalisation de AN en base
orthonormée) et on a donc
(26) EN (t) ≤ t sup ∥rN (t′ )∥
t′ ∈[0,t]

(une analyse plus fine montrerait qu’on peut remplacer le facteur t par une constante). Il ne reste
plus qu’à borner ∥rN (t)∥, ce qui se fait comme précédemment.
11.2. Discrétisation en temps. Il reste maintenant à discrétiser en temps l’équation
dUN
(27) = AN UN , (UN (0))n = u0 (xn ).
dt
Comme on l’a déjà vu, la matrice AN a des valeurs propres d’ordre 1/(∆x)2 ; plus précisément,
sa valeur propre maximale est inférieure à 4/(∆x)2 . C’est donc une équation raide (possédant des
échelles de temps très différentes). Comme on l’a vu au chapitre EDO, si on discrétise cette équation
par une méthode d’Euler, il est nécessaire de prendre un pas de temps
2 1
∆t < < (∆x)2
λmax 2
pour ne pas que la solution explose. C’est la condition CFL (Courant-Friedrichs-Lewy) pour
l’équation de la chaleur. Cette condition signifie en particulier que la convergence du schéma total
(discrétisé en espace et en temps) ne converge pas inconditionnellement si ∆x, ∆t → 0 : il faut que
21
la convergence se fasse en respectant la condition CFL. Comme vu au chapitre EDO, une alternative
est d’utiliser un schéma implicite, qui garantit la convergence si ∆x, ∆t → 0 dans n’importe quel
ordre.
TD/TP
— Refaire la même analyse pour l’équation
∂u ∂4u
=− 4
∂t ∂x
Que devient la condition CFL ?
— Que se passe-t-il si dans l’analyse d’erreur on utilise la majoration
d∥EN ∥
(28) ≤ ∥AN ∥∥EN ∥ + ∥rN (t)∥
dt
et on essaie d’appliquer un lemme de Gronwall ?
— Refaire la même analyse sur l’équation
∂u
(29) = u′′ + u.
∂t
L’analyse marche-t-elle sur l’équation suivante ?
∂u
(30) = u′′ + u′
∂t
— Charger une image grâce à la commande image.imread du module matplotlib, et calculer
la solution de l’équation de la chaleur avec comme condition initiale l’image vue comme une
fonction de [0, 1]2 dans R. Que constatez-vous ?

22