0% ont trouvé ce document utile (0 vote)

46 vues124 pages

Methodes Numeriques-CIFAD

Ce document présente un cours de méthodes numériques pour les étudiants en mathématiques, axé sur l'utilisation pratique des algorithmes en tenant compte des contraintes de l'arithmétique à précision finie. Il couvre des sujets tels que la formule de Taylor, la représentation des nombres, la résolution de systèmes linéaires, et les équations différentielles ordinaires, tout en intégrant des exercices pour l'auto-évaluation. Le cours vise à harmoniser l'arithmétique finie avec les algorithmes pour optimiser l'utilisation des calculateurs numériques.

Transféré par

marcellinbodjrenou82

Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.

Formats disponibles

Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd

0% ont trouvé ce document utile (0 vote)

46 vues124 pages

Methodes Numeriques-CIFAD

Transféré par

marcellinbodjrenou82

Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.

Formats disponibles

Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd

UNIVERSITE DES SCIENCES ET TECHNOLOGIES

DE CÔTE-D’IVOIRE

LICENCE
DE
MATHÉMATIQUES

TROISIÈME ANNÉE

Méthodes numériques

• Formation À Distance
• Formation Présentielle
• Formation Continue
Table des matières

Avant-Propos 3

1 Un outil de base : la formule de Taylor 4

1.1 Théorème de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Les formes usuelles du théorème de Taylor . . . . . . . . . . . . . 5
1.3 Applications à l’approximation numérique des fonctions . . . . . 8

2 Représentation des nombres et notion d’erreurs 11

2.1 Représentation d’un réel . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.1 Numération à base β . . . . . . . . . . . . . . . . . . . . . 12
2.1.2 Changement de base de numération . . . . . . . . . . . . 13
2.1.3 Arithmétique à précision finie . . . . . . . . . . . . . . . 17
2.1.4 La norme IEEE . . . . . . . . . . . . . . . . . . . . . . . . 18
2.1.5 Opérations élémentaires . . . . . . . . . . . . . . . . . . . 20
2.1.6 La perte de chiffres significatifs . . . . . . . . . . . . . . . 21
2.1.7 Comment corriger la perte de précision dans une soustrac-
tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.2 Notion d’erreur . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.2.1 Erreurs de modélisation . . . . . . . . . . . . . . . . . . . 23
2.2.2 Erreurs de Troncature . . . . . . . . . . . . . . . . . . . . 24

3 Résolution des systèmes linéaires 26

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.2 Méthodes directes de résolution de systèmes linéaires Au = b . . 27
3.2.1 Systèmes triangulaires . . . . . . . . . . . . . . . . . . . . 27
3.2.2 Elimination naı̈ve de Gauss . . . . . . . . . . . . . . . . . 29
3.2.3 Triangularisation . . . . . . . . . . . . . . . . . . . . . . . 29
3.2.4 Factorisation LU . . . . . . . . . . . . . . . . . . . . . . . 32
3.2.5 Elimination de Gauss avec stratégie de pivot . . . . . . . 36
3.2.6 Système linéaires spéciaux : factorisations de Cholesky,
LDM t ... . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.2.7 Analyse d’erreur d’arrondi . . . . . . . . . . . . . . . . . 48
3.3 Méthodes itératives de résolution de systèmes linéaires Au = b . 51
3.3.1 Contexte général . . . . . . . . . . . . . . . . . . . . . . . 51

1
3.3.2 Méthodes itératives standard : Jacobi, Gauss Seidel, méthodes
de relaxation . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.3.3 Convergence des méthodes . . . . . . . . . . . . . . . . . 60
3.3.4 Effet de l’arithmétique finie . . . . . . . . . . . . . . . . . 60
3.4 Accélération de convergence . . . . . . . . . . . . . . . . . . . . 63

4 Résolution numérique des équations différentielles ordinaires 64

4.1 Equation différentielle et classification . . . . . . . . . . . . . . . 64
4.1.1 Equation différentielle et solution . . . . . . . . . . . . . . 65
4.1.2 Problème de Cauchy . . . . . . . . . . . . . . . . . . . . . 66
4.1.3 Problème de Cauchy d’ordre 1 . . . . . . . . . . . . . . . 68
4.2 Résolution d’un problème de Cauchy . . . . . . . . . . . . . . . 70
4.2.1 Résultats d’existence et unicité . . . . . . . . . . . . . . . 70
4.2.2 Résolution numérique des EDO et l’intégration . . . . . . 71
4.3 Méthodes de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.3.1 Méthode d’Euler . . . . . . . . . . . . . . . . . . . . . . . 73
4.3.2 Majoration d’erreur de la méthode d’Euler . . . . . . . . 76
4.3.3 Méthode de Taylor d’ordre n . . . . . . . . . . . . . . . . 76
4.4 Méthodes de Runge Kutta . . . . . . . . . . . . . . . . . . . . . . 83
4.4.1 Rappel de la formule de Taylor à deux variables . . . . . 84
4.4.2 Méthodes de Runge Kutta d’ordre 2 . . . . . . . . . . . . 85
4.4.3 Méthodes de Runge Kutta d’ordre 4 . . . . . . . . . . . . 86
4.5 Méthodes multipas . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.5.1 Méthodes d’Adams-Bashford . . . . . . . . . . . . . . . . 95
4.5.2 Méthodes d’Adams-Moulton . . . . . . . . . . . . . . . . 98
4.6 Méthodes prédicteur-correcteur . . . . . . . . . . . . . . . . . . . 101
4.7 Notion de stabilité de solution des ODE . . . . . . . . . . . . . . 101
4.8 Les systèmes d’équations différentielles . . . . . . . . . . . . . . . 103
4.9 Travaux Pratiques . . . . . . . . . . . . . . . . . . . . . . . . . . 104

5 Rappel et complément d’algèbre linéaire 105

5.1 Principales définitions d’algèbre linéaire . . . . . . . . . . . . . . 105
5.1.1 Vecteurs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
5.1.2 Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
5.1.3 Operations sur les matrices . . . . . . . . . . . . . . . . . 109
5.1.4 Matrices carrées . . . . . . . . . . . . . . . . . . . . . . . 110
5.1.5 Réduction de matrices . . . . . . . . . . . . . . . . . . . . 114
5.2 Normes vectorielles et matricielles . . . . . . . . . . . . . . . . . 114
5.2.1 Normes vectorielles . . . . . . . . . . . . . . . . . . . . . . 114
5.2.2 Normes matricielles . . . . . . . . . . . . . . . . . . . . . 115
5.2.3 Suites vectorielles, suites matricielles . . . . . . . . . . . . 117

6 Conditionnement 119
6.1 Conditionnement d’une matrice . . . . . . . . . . . . . . . . . . 119
6.2 Conditionnement d’un système linéaire . . . . . . . . . . . . . . 121

2
Avant-Propos

Les calculateurs numériques ont envahi notre quotidien comme le témoignent

nos artefacts téléphoniques, nos automobiles et j’en passe. Donc d’une façon ou
d’une autre nous sommes confrontés à des expériences numériques et le marché
foisonne de logiciels offrant des services plus ou moins adaptés. Dans un tel
environnement quelle pertinence peut avoir un cours de méthodes numériques,
disons plutôt un cours de calcul numérique pour mieux cadrer le sujet ?
Cette introduction au calcul numérique n’est pas un cours d’analyse numérique
qui met généralement l’accent sur les développements théoriques au détriments
des démarches pratiques d’utilisation des méthodes étudiées. Le présent cours est
surtout centré sur la prise en compte des contraintes imposées par l’arithmétique
à précision finie dans la pratique du calcul. Un calculateur numérique est un outil
puissant mais imparfait de nature et notre objectif ici est d’en tirer un meilleur
parti. Ce pari sera gagné si à l’issue du cours, la nécessaire harmonie du tandem
arithmétique finie/algorithme est bien perçue.
Une calculatrice électronique est indispensable dans les manipulations lors
des activités pédagogiques, évaluations incluses. Il est à noter que ce n’est pas
la calculatrice qui fera le calcul, il ne sera qu’une aide au calcul et vous serez le
simulateur de la machine numérique.
La maı̂trise d’un environnement de programmation mathématique comme le
MatLab ou le SciLab facilitera la prise en main des algorithmes étudiés tout au
long de ce cours. SciLab est un produit en libre diffusion que nous recommandons
particulièrement car il n’exige qu’un investissement dans l’apprentissage. Il est
disponible sur le site web www.scilab.org.
Cette note de cours comporte quatre chapitres qui couvrent des thèmes dont
le choix n’est motivé que par la poursuite de l’harmonie du tandem arithmétique
finie/algorithme. On trouvera en annexe deux chapitres complémentaires par
souci d’auto-suffisance du document. Des exercices sont disséminés dans tout le
texte pour permettre l’auto-évaluation de l’apprentissage.
Vos critiques et suggetions seront utiles pour améliorer ce document.

Dr. Auélien GOUDJO

[email protected], [email protected]

3
Chapitre 1

Un outil de base : la
formule de Taylor

La formule de Taylor constitue un outil fondamental d’analyse tout au long

de ce cours. Elle présente un interèt théorique et pratique et on la retrouve sous
des formes très variées.

1.1 Théorème de Taylor

Sa formulation universelle est celle dite avec reste intégral. Cette forme avec
reste intégral se prête bien à une généralisation aux cas des fonctions vectorielles
et des fonctions à plusieurs variables. Elle peut s’énoncer :
I ∗ . Soit
I . Soit m ∈ N
Théorème 1.1.1 Soit I un intervalle ouvert non vide de R
m
f ∈ C (I) et x ∈ I, alors pour tout y ∈ I on a
m−1 y
(y − x)k (k) (y − t)m−1 (m)
X Z
f (y) = f (x) + f (t)dt
k! x (m − 1)!
k=0

Remarque :
En posant h = y − x, on a
m−1 x+h
hk (k) (x + h − t)m−1 (m)
X Z
f (x + h) = f (x) + f (t)dt
k! x (m − 1)!
k=0

m−1 1
hk (k) hm
X Z
f (x + h) = f (x) + (1 − λ)m−1 f (m) (x + λh)dλ
k! (m − 1)! 0
k=0

Remarque :
In
f ∈ C m (I) s’entend une fonction m fois continûment dérivable à valeur dans R

4
Preuve
On procède par réccurence sur m.
Pour m = 1, on a Z y
f (y) = f (x) + f 0 (t)dt
x
puisque f est une primitive de f 0 .
Pour m > 1, considérons
Z y Z y
(y − t)m−1 (m) (y − t)m−1 (m−1)
Rm = f (t)dt = df (t)
x (m − 1)! x (m − 1)!
En posant ( m−1
u(t) = (y−t)
(m−1)!
dv(t) = f (m) (t)dt
on a ( m−2
du(t) = − (y−t)
(m−2)! dt
(m−1)
v(t) =f (t)
et
h iy
Ry (y−t)m−1 (m) (y−t)m−1 (m−1) Ry m−2

x (m−1)! f (t)dt = (m−1)! f (t) + x (y−t)

(m−2)! f
(m−1)
(t)dt
m−1
t=x m−2
(y−x) (m−1)
R y (y−t) (m−1)
= − (m−1)! f (x) + x (m−2)! f (t)dt

On a donc la relation
(
R1 = f (y) − f (x)
m−1
Rm = − (y−x)
(m−1)! f
(m−1)
(x) + Rm−1 ∀m > 1

Donc par récurrence on a

m−1 m−1
X (y − x)k (k) X (y − x)k
Rm = − f (x) + R1 = − f (k) (x) + f (y) − f (x)
(k)! (k)!
k=1 k=1

d’où
m−1 y
(y − x)k (k) (y − t)m−1 (m)
X Z
f (y) = f (x) + f (t)dt
(k)! x (m − 1)!
k=0

1.2 Les formes usuelles du théorème de Taylor

La formule de Taylor avec reste de Young est très utilisée avec deux présentations.

Corollaire 1.2.0.1 Soit f ∈ C m (I), alors pour tout x 6= y ∈ I, on a

m−1
X (y − x)k (k)
f (y) = f (x) + O((y − x)m )
k!
k=0

5
Preuve

En effet
m−1
X (y − x)k (k)
f (y) = f (x) + Rm
k!
k=0
où y
(y − t)m−1 (m)
Z
Rm = f (t)dt
x (m − 1)!
Comme f (m) est bornée parce que continue, alors il existe c > 0 tel que f (m) (t) ≤
cm! pour tout t ∈ (x, y). Ainsi
Z y
(y − t)m−1 m
|Rm | ≤ cm! dt = c |y − x|
x (m − 1)!
On déduit alors que
Rm
≤c
(y − x)m
d’où le resultat
Corollaire 1.2.0.2 Soit f ∈ C m (I), alors pour tout x 6= y ∈ I, on a
m
X (y − x)k
f (y) = f (k) (x) + o((y − x)m )
k!
k=0

Preuve
En effet d’après le théorème de Taylor, on a
m−1
X (y − x)k (k)
f (y) = f (x) + Rm
k!
k=0

avec y
(y − t)m−1 (m)
Z
Rm = f (t)dt
x (m − 1)!
(m)
Puisque f est continue en x on a
m−1
Ry Ry (y−t)m−1
Rm = x (y−t)
(m−1)! f
(m)
(x)dt + x (m−1)! {f
(m)
(t) − f (m) (x)}dt

(y−x)m (m) Ry (y−t)m−1 (m)

= m! f (x) + x (m−1)! {f (t) − f (m) (x)}dt

Or f (m) continue en x alors pour tout > 0 il existe δ > 0 tel que pour tout
t on a
|t − x| < δ ⇒ f (m) (t) − f (m) (x) < m!
Ainsi
y
(y − x)m (m) (y − t)m−1
Z
m
Rm − f (x) < m! dt = |y − x|
m! x (m − 1)!

6
Donc
m
X (y − x)k m
f (y) − f (k) (x) < |y − x|
k!
k=0
D’où la conclusion
La formulation la plus pratique de la formule de Taylor est celle avec reste
de Lagrange qui s’énonce :
Corollaire 1.2.0.3 Soit f ∈ C m (I) à valeur réelle, alors pour tout x 6= y ∈ I,
il existe c ∈ (x, y) tel que
m−1
X (y − x)k (k) (y − x)m (m)
f (y) = f (x) + f (c)
k! m!
k=0

Preuve

On utilisera à cet effet la seconde formule de la moyenne énoncée dans la

proposition 1.2.1 . En effet en posant
(y − t)m−1
φ(x) =) , ψ(t) = f (m) (t)
(m − 1)!
on a Ry (y−t)m−1 (m) Ry Ry
x (m−1)! f (t)dt = x
φ(t)ψ(t)dt = ψ(c) x
φ(t)dt

Ry (y−t)m−1
= f (m) (c) x (m−1)! dt

m
= f (m) (c) (y−x)
m!
Remarque :
La formule de Taylor avec reste de Lagrange n’est pas applicable aux fonctions
vectorielles. Prenons pour exemple f (t) = eıt et m = 1. Si la formule est vraie
alors pour tout x ∈ R
I on aurait l’existence d’un c ∈ (0, x) donnant
f (x) − f (0) = xf 0 (c)
c’est à dire
eıx − 1 = ıxeıc
Ainsi pour x = 2π on a
0 = 2πıeıc
En d’autres termes on a
2π = |2πıeıc | = 0
ce qui est absurde.
Proposition 1.2.1 ∀φ, ψ ∈ C[a, b] avec φ ≥ 0, il existe c ∈ [a, b] tel que
Z b Z b
ψ(t)φ(t)dt = ψ(c) φ(t)dt
a a

7
Preuve

Comme ψ est continue sur [a, b] alors ψ([a, b]) = [m, M ] Donc

mφ(t) ≤ φ(t)ψ(t) ≤ M φ(t), ∀t ∈ [a, b]

Ainsi Z b Z b Z b
m φ(t)dt ≤ φ(t)ψ(t)dt ≤ M φ(t)dt
a a a
Rb
En supposant φ non identiquement nulle, on a a
φ(t)dt > 0 et
Rb
a
φ(t)ψ(t)dt
m≤ Rb ≤M
a
φ(t)dt

D’après le théorème des valeurs intermédiaires, il existe c ∈ [a, b] tel que

Rb
a
φ(t)ψ(t)dt
ψ(c) = Rb
a
φ(t)dt

Ce qui achève la démonstration

1.3 Applications à l’approximation numérique

des fonctions
La plupart des bibliothèques numériques disponibles sur nos calculatrices
implémentent les fonctions classiques en utilisant la formule de Taylor car on se
ramène à une évaluation de polynôme.
Exemple √
Donner une √ approximation de 1.2 en utilisant le développement de Taylor à
l’ordre 3 de 1 +√x au voisinage de 0.
Pour f (x) = 1 + x, f ∈ C ∞ (] − 1 , ∞[) et on a
√
f (x) = 1 + x f (0) = 1
0
1 − 21 0
f (x) = 2 (1 + x) f (0) = 12
3
−
f (2) (x) = − 41 (1 + x) 2 f (2) (0) = − 14
5
−
f (3) (x) = 38 (1 + x) 2 f (3) (0) = 38
7
15 −
f (4) (x) = − 16 (1 + x) 2

Ainsi d’après le corollaire 1.2.0.3 il existe c ∈]1 , 1.2[ tel que

√ 2 3 4
1.2 = f (0.2) = f (0) + 0.2f 0 (0) + 0.2
2 f
(2)
(0) + 0.2
6 f
(3)
(0) + 0.224 f
(4)
(c)
− 72
= 1 + 0.1 − 0.005 + 0.0005 − 0.0000625(1 + c)
−7
= 1.0955 − 0.0000625(1 + c) 2

8
√
d’où 1.2 ≈ 1.0955 dont la valeur exacte à 20 décimales est 1.09544511501033215
Exemple Z 0.1
Donner une approximation de cos xdx à l’aide d’un développement de Tay-
0
lor à l’ordre 3 au voisinage de 0.
Pour x ∈]0 , 0.1[, d’après le corollaire 1.2.0.3 il existe c(x) ∈]0 , x[ tel que
1 1
cos x = 1 − x2 + x4 cos(c(x))
2 24
R 0.1 R 0.1 1 2
1
R 0.1 4
0
cos xdx = 0 1 − 2 x dx + 24 0
x cos(c(x))dx
3 0.1
h i
0.1
= x − x6 1
x4 cos(c(x))dx
R
+ 24 0
0
−3 0.1
= 599 1
x4 cos(c(x))dx
R
6 10 + 24 0
Z 0.1
d’où cos xdx ≈ 0.09983̄ dont la valeur exacte est sin(0.1) qui, à 20 décimales,
0
est 0.09983341664682815
La précision des résultats soulèvent deux questions que nous allons survoler
ici à travers des exemples illustratifs.
La question directe consiste à établir une majoration d’erreur commise sachant
l’ordre de développement
Exemple √
Reprenant l’approximation de 1.2 par le développement de Taylor à l’ordre 3
√ −7
de f (x) = 1 + x, le terme d’erreur est R = −0.0000625(1 + c) 2 . Donc
− 72
|R| ≤ max 0.0000625(1 + c) = 0.0000625
0≤c≤0.2

On peut comparer à
√
1.2 − 1.0955 = |1.09544511501033215 − 1.0955| = 0.00005488498966777
Z 0.1
Exercice 1 Donner une majoration de l’erreur lors de l’approximation de cos xdx
0
à l’aide d’un développement de Taylor à l’ordre 3 au voisinage de 0.
La deuxième question est plutôt indirecte. Il s’agit de déterminer l’ordre de
développement de Taylor pour un seuil d’erreur admissible.
Exemple
A quel ordre n faut-t-il faire le développement de Taylor de cos x au voisinage
de 0 pour évaluer cos 0.1 avec une erreur inférieure à 10−8 ?
D’après le corollaire 1.2.0.3 on a
n k n
X (0.1) π (0.1) π
cos 0.1 = cos (k ) + Rn avec Rn = cos (c + n )
k! 2 n! 2
k=0

Comme n n
(0.1) π (0.1)
|Rn | ≤ max cos (c + n ) ≤
n! 0≤c≤0.1 2 n!

9
pour avoir |Rn | ≤ 10−8 il suffit que
n
(0.1)
≤ 10−8 ⇔ 108 ≤ n!10n ⇔ n ≥ 6
n!
car 5!105 = 12 106 < 108 < 72 107 = 6!106

Exercice 2 Combien de termes du développement de Taylor de f (x) = ex au

voisinage de 0 pour que l’approximation de e soit exacte sur 10 décimales.

Remarque :
En pratique, pour réduire la propagation des erreurs d’arrondis que nous verrons
dans le chapitre 2 relatif à la représentation des nombres, on utilise l’algorithme
de Horner pour évaluer les polynômes.
Xn
Pour un polynôme P de degré n ∈ N I ∗ donné par P (x) = ai xi , l’algo-
i=0
rithme de Horner se ramène l’évaluation de P (x) au calcul des termes de la
n
suite finie (pk )k=0 donnée par la formule de récurrence :

p0 = an
pk = pk−1 x + an−k avec k = 1, . . . , n

pn est la valeur de P (x).

On parle aussi de l’algorithme de la multiplication emboı̂tée.

Exercice 3 Utiliser √le développement de Taylor

√ à l’ordre
√ 2 au voisinage de 0
de la fonction x 7→ 1 + x pour calculer 1.001 et 0.999. Préciser l’erreur
commise lors de ces evaluations.

Exercice 4 A l’aide de deux polynômes de Taylor de degré 2 calculer une ap-

proximation du réel √
K = 1.05 + 0.963
Donner une majoration de l’erreur commise

10
Chapitre 2

Représentation des nombres

et notion d’erreurs

La conscience que nous devons avoir de la représentation des nombres sur

les outils de calcul que sont les calculatrices électronique et les ordinateurs est
la base du bon usage que nous pouvons en faire. Ce chapitre est destiné à faire
ressortir les limitations naturelles de ces outils.

2.1 Représentation d’un réel

Dans le système de numération décimale traditionnelle qui est une numération
de position, nous utilisons dix symboles ou chiffres qui sont 0, 1, 2, 3, 4, 5, 6, 7,
8, 9.
La suite de chiffres an an−1 . . . a2 a1 a0 représente le nombre entier
0 1 2 n−1
an an−1 . . . a2 a1 a0 =a0 × 10 + a1 × 10 + a2 × 10 + . . . + an−1 × 10 + an × 10n
P n k
= k=0 ak 10

Exemple

47295 = 40000 + 7000 + 200 + 90 + 5

= 4 × 104 + 7 × 103 + 2 × 102 + 9 × 101 + 5 × 100
La suite de chiffres 0.b1 b2 b3 . . . représente le nombre décimal
−1
0.b1 b2 b3 . . . 1 × 10
= bP + b2 × 10−2 + b3 × 10−3 + . . .
∞
= k=1 bk 10−k

Exemple
6 2 7 5
0.6275 = 10 + 100 + 1000 + 10000
= 6 × 10 + 2 × 10 + 7 × 10−3 + 5 × 10−4
−1 −2

11
De façon générale, en numération décimale, un nombre réel est de la forme
n
X ∞
X
(an an−1 . . . a2 a1 a0 .b1 b2 b3 . . .)10 = ak 10k + bk 10−k
k=0 k=1
Pn k
P∞ −k
où k=0 ak 10 est dite partie entière et k=1 bk 10 est dite partie fraction-
naire

2.1.1 Numération à base β

Bien que le système décimal soit le système de numération traditionnelle,
l’essor de l’informatique a permis la vulgarisation des systèmes binaire (base 2),
octal (base 8) et hexadécimal (base 16).
Dans un système à base β, les chiffres utilisés sont 0, 1, . . . , β − 1. Lorsque
β > 10 il est alors nécessaire d’introduire des symboles pour représenter les
chiffres 10, 11, . . . , β − 1. On utilise un point comme séparateur de la partie
entière de la partie fractionnaire. Ainsi dans le système à base β un nombre est
représenté comme suit
n
X ∞
X
(an an−1 . . . a2 a1 a0 .b1 b2 b3 . . .)β = ak β k + bk β −k
k=0 k=1

Exemple

(21467)8 = 2 × 84 + 1 × 83 + 4 × 82 + 6 × 8 + 7
= 7 + 8(6 + 8(4 + 8(1 + 8(2)))) = 9015
Exemple

(0.36207)8 = 3 × 8−1 + 6 × 8−2 + 2 × 8−3 + 0 × 8−4 + 7 × 8−5

= 8−5 (7 + 8(0 + 8(2 + 8(6 + 8(3))))) = 32768
15495

= 0.47226987 . . .
Remarque :
Dans ces deux exemples on voit le rôle que peut jouer l’algorithme de Horner
pour l’évaluation d’un polynôme. Cet algorithme consiste à transformer un
polynôme en une suite emboı̂tée de produit de la forme :
p(x) = a0 + a1 x + a2 x2 + . . . + an−1 xn−1 + an xn
= a0 + x(a1 + x(a2 + . . . + x(an−1 + x(an ))))
Cet algorithme se traduit par le pseudocode :
real array (ai )0:n
integer i, n
real p, x
p ← an
for i = n − 1 to 0 step -1 do
p ← ai + x . p
end for

12
Exercice 5 Donner la valeur décimale du réel x dans chacun des cas :
1. x = (4506)7
2. x = (2.1211)3
3. x = (361.002)8
4. x = (1253.75)9
5. x = (0.3042)5

2.1.2 Changement de base de numération

Le changement de base de numération suppose en général que nous passons
d’une base β à une base α. Du fait que la base 10 est notre base naturelle,
nous procèderons en deux séquences en transitant par la base 10. Nos aurons
un passage de la base β à une base 10, suivi d’un passage de la base 10 à une
base α.
La première séquence est explorée dans la sous section précédente en utilisant
l’algorithme de Horner.
Regardons maintenant le passage de la base 10 à une autre.
Soit x = (an an−1 . . . a2 a1 a0 .b1 b2 b3 . . .)α = I(x) + F(x) avec
Pn
I(x) = (an an−1 . . . a2 a1 a0 )α = k=0 ak αk
et P∞
F(x) = (0.b1 b2 b3 . . .)α = k=1 bk α−k

Soit β une autre base. Nous recherchons N = I(x) et F = F(x) sous la

forme
N = (cm cm−1 . . . c2 c1 c0 )β = c0 + β(c1 + β(c2 + . . . + β(cm−1 + β(cm ))))
F = (0.d1 d2 . . . dl . . .)β = β −1 (d1 + β −1 (d2 + . . . + β −1 (dl−1 + β −1 (dl + . . .))))

Conversion de la partie entière

On observe que les Ci sont obtenus comme reste de division euclidienne
répétée par β.

Dividende Quotient Reste

N q0 = c1 + β(c2 + . . . + β(cm−1 + β(cm ))) c0
q0 q1 = (c2 + . . . + β(cm−1 + β(cm )) c1
q1 q2 = (c3 + . . . + β(cm−1 + β(cm )) c2
.. .. ..
. . .
qm−2 qm−1 = cm cm−1
qm−1 qm = 0 cm

Exemple

13
Convertir l’entier N = 397 en base 5
Dividende Quotient Reste
397 79 2
79 15 4
15 3 0
3 0 3

On conclut que N = (3042)5

Conversion de la partie fractionnaire

On observe que les di sont obtenus par un processus de multiplication par
β. En effet :

F0 = F F1 = F(βF0 ) d1 = I(βF0 )
F1 F2 = F(βF1 ) d2 = I(βF1 )
F2 F3 = F(βF2 ) d2 = I(βF2 )
.. .. ..
. . .
Fl−2 Fl−1 = F(βFl−2 ) dl−1 = I(βFl−2 )
Fl−1 Fl = F(βFl−1 ) dl = I(βFl−1 )
.. .. ..
. . .
Exemple
Convertir F = 0.5318 en base 5.

Fi βFi Fi+1 = F(βFi ) di+1 = I(βFi )

0.5318 2.659 0.659 2
0.659 3.295 0.295 3
0.295 1.475 0.475 1
0.475 2.275 0.275 2
0.275 1.375 0.375 1
0.375 1.875 0.875 1
0.875 4.375 0.375 4
0.375 1.875 0.875 1
0.875 4.375 0.375 4
... ... ... ...
d’où F = (0.231211414...)5 avec la séquence 14 terminale qui se répète indéfiniment.

Exercice 6 Pour le réel x faire les conversions dans chacun des cas :
1. x = (647)10 = (. . .)7 = (. . .)2
2. x = (. . .)10 = (34.012)5 = (. . .)3
3. x = (. . .)10 = (. . .)7 = (3A.B)12

14
Conversion base 10 ↔ 8 ↔ 2
Le passage du système décimal au sytème octal utilise les processus décrits
plus haut. Les conversions du système octal au binaire et vice-versa et du
système hexadécimal au binaire et vice-versa se font de façon beaucoup plus
simple et par substitution suivant des tableau de correspondance.
En effet
P3m+2
(c3m+2 c3m+1 c3m . . . c2 c1 c0 )2 = k=0 ck 2k
Pm
= i=0 c3i+2 23i+2 + c3i+1 23i+1 + c3i 23i
Pm
= i=0 c3i+2 22 + c3i+1 21 + c3i 20 23i
Pm
= i=0 c3i+2 22 + c3i+1 21 + c3i 20 8i
Pm
= i=0 di 8 i

avec di = c3i+2 22 + c3i+1 21 + c3i 20 ∈ {0, 1, 2, . . . , 7}

P4m+3
(c4m+3 c4m+2 c4m+1 c4m . . . c3 c2 c1 c0 )2 = k=0 ck 2 k
Pm
= i=0 c4i+3 24i+3 + c4i+2 24i+2 + c4i+1 24i+1 + c4i 24i
Pm
= i=0 c4i+3 23 + c4i+2 22 + c4i+1 21 + c4i 20 24i
Pm
= i=0 di 16i

avec di = c4i+3 23 + c4i+2 22 + c4i+1 21 + c4i 20 ∈ {0, 1, 2, . . . , 9, A, B, . . . , F }

Les tables de correspondance sont :

Octal/binaire
Octal 0 1 2 3 4 5 6 7
binaire 000 001 010 011 100 101 110 111
Hexadécimal/binaire
hexadécimal 0 1 2 3 4 5 6 7
binaire 0000 0001 0010 0011 0100 0101 0110 0111
hexadécimal 8 9 A B C D E F
binaire 1000 1001 1010 1011 1100 1101 1110 1111
Exemple
Convertir x = (347.12537)8 en base 2 et en base 16.
On a :
x = (347.12537)8
= (11100111.001010101011111)2
= (E7.2ABE)16

15
Exemple
Convertir x = (2A.EF 1)16 en base 2 et en base 8.
On a :
x = (2A.EF 1)16
= (101010.111011110001)2
= (52.7341)8
Remarque :
L’approche intuitive développée dans la présente section est soutendue par le
théorème 2.1.1 qui s’énonce comme suit :
Théorème 2.1.1 Soit α ∈ N I ∗ − {1}. Soit x ∈ R I ∗+ . Il existe un unique n ∈ Z
et une unique suite (bk )k≥1 d’éléments de {0, 1, 2, . . . , α − 1} tel quel b1 6= 0 et
∞
X bk
x = αn
αk
k=1

I ∗ , on a l’encadrement :
Plus précisement si t ∈ N
t t
!
X
n bk X bk 1
xt = α ≤ x < x̄t = αn + t
αk α k α
k=1 k=1
n
On note x = (0.b1 b2 . . . bt . . .)α × α et on parle de notation scientifique nor-
malisée.
Preuve
I ∗+ , comme RI ∗+ , × , ≤ est un groupe archimédien alors il existe

α >1 e t x∈R
unique n ∈ Z tel que
αn−1 ≤ x < αn
x x
On en déduit que 1 ≤ n−1 < α c’est à dire que n−1 ∈ [1 , α[= ∪α−1 j=1 [j , j + 1[.
α α
x
Ainsi il existe un unique b1 ∈ {1, 2, . . . , α − 1} tel que n−1 ∈ [b1 , b1 + 1[
α
et on a
x − b1 αn−1
b1 αn−1 ≤ x < (b1 + 1)αn−1 ⇔ 0≤ n−2
<α
αn−1
x − b1 α α−1
⇔ ∈ [0 , α[= ∪j=0 [j , j + 1[
αn−2
x − b1 αn−1
Il existe donc un unique b2 ∈ {0, 1, . . . , α − 1} tel que ∈ [b2 , b2 + 1[
αn−2
et on a
b1 αn−1 + b2 αn−2 ≤ x < b1 αn−1 + (b2 + 1)αn−2
m
n b 1 b2 n b1 b2 1
α + 2 ≤x<α + 2+ 2
α1 α α1 α α
! m !
2 2
X b i
X bi 1
αn i
≤ x < αn + 2
i=1
α i=1
αi α

16
I ∗ on a l’encadrement
On déduit alors par récurrence que pour tout t ∈ N
t
! t
!
n
X bi n
X bi 1 bi ∈ {0, 1, . . . , α − 1}
xt = α ≤x<α + t = x̄t pour
α i α i α et b1 6= 0
i=1 i=1

On vérifie facilement que la suite (xt )t≥1 est croissante et majorée par x et que
la suite (x̄t )t≥1 est décroissante et minorée par x. Elles convergent donc dans R
I
avec lim xt ≤ x ≤ lim x̄t .
t→∞ t→∞
1
Par ailleurs on a |xt − x̄t | = t → 0. Ainsi lim xt = lim x̄t = x.
α t→∞ t→∞

2.1.3 Arithmétique à précision finie

Du fait de l’arithmétique à précision finie ou arithmétique flottante, une
machine ne peut représenter qu’une partie F de R I . Cette partie F est caractérisé
par quatre entiers : la base β, la précision t, la limite inférieure L et la limite
supérieure U . En représentation scientifique normalisée on a G = F ∪ {0} avec

I , |f | = 0.d1 d2 . . . dt × β e ,
F = {f ∈ R di ∈ {0, 1, . . . β − 1}, d1 6= 0, L ≤ e ≤ U }

Remarque :
Les éléments non nul de G (ensemble des nombres-machine) sont de la forme :

f = ∓0.d1 d2 . . . dt × β e , 0 ≤ di < β, d1 6= 0, L≤e≤U

Remarque :
si 0 6= f ∈ G alors m ≤ |f | ≤ M avec m = β L−1 et M = β U (1 − β −t )
Pour tout x ∈ F il existe r > 0 et e ∈ Z tel que |x| = rβ e avec ≤ r < 1−β −t
Soit Ĝ le sous ensemble de RI défini par

Ĝ = {x ∈ R
I , m ≤ |x| ≤ M }

Le processeur I →G
arithmétique est la donnée de Ĝ et de l’opérateur f l : R
defini par

 c∈G tel que |c − x| minimal : arithmétique arrondie
f l(x) = ou
c∈G tel que |c − x| minimal et |c| ≤ |x| : arithmétique tronquée


Ainsi
f l(x) = x(1 + ) , || ≤ u
avec u la précision machine c’est à dire la plus grande erreur relative que l’on
puisse commettre en représentant un nombre réel sur un ordinateur.
Exemple
Lister les éléments de l’ensemble A des réels strictement positifs de la forme
0.b1 b2 × 2k avec bi ∈ {0, 1} et k ∈ {−1, 0, 1}

17

b1 b2
x ∈ A si x = + 2k avec (b1 , b2 ) ∈ {(0 , 1), (1 , 0), (1 , 1)} et k ∈
2 4
{−1, 0, 1}.
1 1 3
Pour k = 0 on a x ∈ , ,
4 2 4

1 1 3
Pour k = −1 on a x ∈ , ,
8 4 8
1 3
Pour k = 1 on a x ∈ , 1,
2 2
1 1 3 1 3 3
En résumé L = , , , , , 1,
8 4 8 2 4 2
Exemple
Que sont les éléments positifs de l’arithmétique binaire flottante dont la forme
normalisée est 0.b1 b2 × 2k avec k ∈ {−1, 0, 1} ?
1 3 1 3 3
En exploitant les calculs précédents on trouve F+ = , , , , 1,
4 8 2 4 2

Théorème 2.1.2 Dans une représentation en base β, la précision machine

vérifie  1 1−t
 2β en arithmétique arrondie
u≤ ou
 1−t
β en arithmétique tronquée

t étant le nombre de chiffres de la mantisse.

Preuve
Soit x un nombre quelconque. Sa représentation exacte en base β est de la
forme :
x = 0, d1 d2 d3 . . . dt dt+1 dt+2 × β e
Par troncature l’erreur absolue est :
∆x = 0, 000 . . . dt dt+1 dt+2 . . . × β e
= 0, dt+1 dt+2 . . . × β e−t
≤ 0, (β − 1)(β − 1)(β − 1) . . . × β e−t

L’erreur relative satisfait donc :

|∆x| 0,(β−1)(β−1)(β−1)...×β e−t
Er (x) = |x| ≤ 0,d1 d2 d3 ...dt dt+1 dt+2 ×β e

0,(β−1)(β−1)(β−1)...×β e−t
≤ 0,1000...×β e

≤ 1 × β 1−t

2.1.4 La norme IEEE

Selon le standard IEEE (Institut for Electrical and Electronic Engeineers),
tous les construiteurs d’ordinateurs utilisent depuis 1985 les mêmes formats de

18
représentation binaire. Par exemple selon la norme IEEE 754, les coprocesseurs
arithmétiques implémentent un réel sur un mot de 32 bits en simple précision
et un mot de 64 bits en double précision.
– Le premier bits indique le signe et est noté s
– les 8 bits suivants en simple précision (les 11 bits suivants en double
précision), indiquent l’exposant et est noté c et dénommé caractéristique
– les 23 bits restants en simple précision (les 52 bits restants en double
précision), indiquent la partie fractionnaire ou mantisse et est noté f .
La base utilisé par le format IEEE est la base 2. Ainsi en double précision
ce système à virgule flottante a une partie fractionnaire f d’une précision d’au
moins 16 chiffres décimaux pendant que l’exposant c va de 0 à 2047. Afin de
représenter les nombres de petite valeur absolue on procède à une normalisation
des exposants pour aller de −1023 à 1024. Dans ce système un nombre à virgule
flottante est sous la forme normalisée :
s
(−1) ∗ 2c−1023 ∗ (1 + f )
En simple précision on a la forme normalisée est donnée par :
s
(−1) ∗ 2c−127 ∗ (1 + f )
avec 0 ≤ c ≤ 255 et 1 ≤ (1.f )2 = (1 + f )2 ≤ 2 − 2−23
Exemple
Considérons le nombre machine :
0 10000000011 10111001000100000000000000000000000000000000000
On a
s=0
c = 1. 210 + 1. 21 + 1. 20 = 1024 + 2 + 1 = 1027
1 3 4 5 8 1 12
f = 1. 12 + 1. 12 + 1. 21 + 1. 12 + 1. 12 + 1.

2

Ce nombre machine représente alors le nombre décimal

s
(−1) ∗ 2c−1023 ∗ (1 + f ) = 27.56640625
Il est à remarquer que ce nombre machine ne représente pas que le décimal
27.56640625. Il représente également tout un intervalle centré sur le décimal
27.56640625.
Les flottants sont souvent donnés sous forme hexadécimale. Ainsi les nombres
machine x = [3A00E100]16 et y = [C1008000]16 représentent respectivement les
nombres machine 32 bits :
32977
– x=
134217728
65954
pour qui s = 0, c = (01110100)2 = 116, 1 + f = (1.01A2)16 =
164
257
– y=−
32
257
pour qui s = 1, c = (10000010)2 = 130, 1 + f = (1.01)16 = 2
16

19
2.1.5 Opérations élémentaires
Les opérations élémentaires sont l’addition, la soustraction, la multiplication
et la division. En Arithmétique flottante elles s’effectuent de la façon suivante :
x+y → f l(f l(x) + f l(y))
x−y → f l(f l(x) − f l(y))
x×y → f l(f l(x) × f l(y))
x÷y → f l(f l(x) ÷ f l(y))
Exemple
Si l’on prend la précision t = 4 et la base β = 10, alors on a
1
3 ×3 → f l(f l( 13 ) × f l(3))
= f l((0.3333 × 100 ) × (0.3000 × 101 ))
= f l(0.09999000 × 101 )
= f l(0.9999 × 100 )
On a une légère perte de précision. Mais la multiplication et la division sont
simples
Exercice 7 Effectuer les opérations suivantes en arithmétique décimale tronquée
à 4 chiffres.
a) (0.4035 × 106 ) × (0.1978 × 10−1 )
b) (0.56789 × 104 ) ÷ (0.1234321 × 10−3 )
Pour l’addition et la soustraction il faut être particulièrement attentif. Ici il
y a lieu de transformer la mantisse du nombre ayant le plus petit exposant de
manière à opérer avec le plus grand exposant.
Exemple
6
En arithmétique décimale arrondie à 4 chiffres, pour x = 43 et y = 11 , calculons
x + y et x − y
x+y → f l(f l( 34 ) + f l( 11
6
))
= ¯ . . .))
f l(f l(1.3̄ . . .) + f l(0.545
= f l(1.333 + 0.5455)
= f l(1.333 + 0.546)
= f l(1.879) = 1.879
62
Comme x + y = alors l’erreur absolue commise est
33
62
e = |x + y − f l(x + y)| = − 1.879 = 0.0002121
33
d’où l’erreur relative
|x + y − f l(x + y)|
er = = 0.1129 10−3
x+y
Exercice 8 Effectuer les opérations x−y, x×y et x÷y en arithmétique décimale
arrondie à 4 chiffres pour x = 43 et y = 11
6
. Calculer l’erreur relative commise
dans chaque cas.

20
2.1.6 La perte de chiffres significatifs
Nous allons observer le phénomène à partir d’un exemple.
Considérons les deux réels suivants :
x = 0.3721448693
y = 0.3720214371

Calculons l’erreur relative commise lors de l’évaluation de x−y avec une arithmétique
décimale tronquée précise à 5 chiffres.
En effet, on a

x̃ = f l(x) = 0.37214
ỹ = f l(y) = 0.37202
x̃ − ỹ = 0.00012
Comme x − y = 0.0001234322 on a
|(x − y) − (x̃ − ỹ)| 0.0000034322 34322
= = = 3 10−2
|x − y| 0.0001234322 1234322
Sachant que
x̃ ỹ
≤ 0.5 10−4 et ≤ 0.5 10−4
x y
alors on observe une perte de précision.
Dans le calcul de f l(x − y) = (x − y)(1 + δ) nous n’avons que 2 chiffres
significatifs or on en espérait jusqu’à 5 dans notre machine.
Par ailleurs on peut corriger cette perte de précision en changeant la précision
de calcul, mais cela coûtera inutilement cher.

Théorème 2.1.3 Soit x et y deux réels flottants normalisés vérifiant

0 < y < x.
Si 2−p ≤ 1 − xy ≤ 2−q pour des entiers p et q, alors au plus p chiffres et au
moins q chiffres significatifs binaires sont perdus lors de la soustraction x − y.

Preuve
t
X bi
Soit x = 2e i
avec b1 6= 0 Si 0 < y < x alors
i=1
2

t t
y X bi X bi
2−p ≤ 1 − ≤ 2−q ⇔ 2−p x ≤ x − y ≤ 2−q x ⇔ 2e i+p
≤ x − y ≤ 2e
i+q
x i=1
2 i=1
2

On observe que pour 1 ≤ i ≤ t on a 1 + p ≤ i + p ≤ t + p et 1 + q ≤ i + q ≤ t + q.

Du fait de ce décalage seuls les i pour lesquels i + p ≤ t ou i + q ≤ t ne sont pas
perdus.
D’où la conclusion
Exemple p
On considère l’expression f (x) = x2 + 1 − 1. Déterninons les réels x pour

21
lesquels une perte de chiffres significatifs est inévitable lors de l’évaluation de
f (x) à l’aide d’une calculatrice électronique.
√
Soit x ∈ R I ∗ ; on observe que X = x2 + 1 > 1 = Y et la calculatrice est
binaire. Donc il y a perte de chiffres significatifs si on en perd au moins un, c’est
Y
à dire qu’on a : 1 − X ≤ 12

1 1 p √
1− √ ≤ ⇔ x2 + 1 ≤ 2 ⇔ |x| ≤ 3
x2 +1 2

L’évaluation de√f (x) génère une perte de chiffres significatifs pour tout x
vérifiant 0 < |x| ≤ 3.

2.1.7 Comment corriger la perte de précision dans une

soustraction
On a diverses techniques qui vont de la rationnalisation au développement de
Taylor en passant par d’éventuels changements de domaine en cas de périodicité.
La stratégie de toutes ces techniques consiste à mettre en œuvre des mécanismes
permettant d’éviter la soustraction.
Nous allons illustration ces mécanismes par des exemples.
Exemple : Rationnalisation
√
Evaluons au voisinage de 0, f (x) = x2 + 1 − 1.
Il suffit d’observer que

x2
f (x) = √
x2 + 1 + 1
Ce qui éclipse la soustraction.
Exemple : Déleveloppement de Taylor

Considérons f (x) = x − sin x au voisinage de 0.

On observe que limx→0 sinx x = 1. Ainsi la meilleure forme réduisant la perte
de précision lors de l’évaluation ded f (x) fait appel à un déleveloppement de
Taylor et on a
x3 x5 x7
f (x) = x − sin x = − + ...
3! 5! 7!
Lors du calcul de fonctions périodiques, les machines utilisent des change-
ments de domaine de manière à se ramener au domaine de manipulation qui leur
est assigné. Il peut en résulter une perte de chiffres significatifs dans d’éventuelle
soustraction.
Exemple : Changement de domaine

Considérons f (x) = sin x pour x ≥ 2π. Soit n ∈ N

I tel que 0 ≥ x − nπ ≥ 2π.
On utilisera alors f (x) = f (x − 2 nπ)

22
Exercice 9 Pour quelles valeurs de x ∈ R I , une perte de chiffres significatifs
est inévitable en utilisant une calculatrice électronique pour l’évaluation de f (x)
et comment peut-t-on contourner cette difficulté dans chacun des cas suivants :
1. f (x) = ex − e−3x
p
2. f (x) = x2 + 1 − x
√ √
3. f (x) = x + 3 − x

2.2 Notion d’erreur

Il est rare que le traitement numérique d’un problème par exemple la résolution
de Ax = b ne donne une solution non entachée d’erreurs.
Une partie importante de l’analyse numérique consiste à contenir les effets
des erreurs ainsi introduites. Ces erreurs ont diverses origines qui vont des incer-
titudes sur les données à la représentation machine de ces données et proviennent
de trois sources essentielles :
– Erreurs de modélisation ;
– Erreurs de troncature ;
– Erreurs de représentation sur ordinateur ;
En effet, mises à part les erreurs d’acquisition des données réelles du problème,
la représentation machine d’un nombre réel n’utilise qu’un nombre limité de
chiffres significatifs. On dit qu’on a une arithmétique à précision limitée et l’er-
reur commise est dite erreur d’arrondi.
La prise en compte de toutes ces altérations dans le traitement de la recherche
de solution du problème Ax = b nous amème en réalité à obtenir la solution du
problème perturbé (A + ∆A)y = (b + ∆b).
En général les perturbations ∆A et ∆b ne sont connues que par leurs majo-
rations. Et comme pour un algorithme particulier, l’analyse de l’accumulation
d’erreurs d’arrondi commises à chaque étape permet de majorer l’erreur finale,
on va se focaliser un peu sur les erreurs d’arrondi avant de regarder les problèmes
de conditionnement.

2.2.1 Erreurs de modélisation

Bien que ne faisant spécifiquement pas l’objet de ce cours, la première étape
de la résolution d’un problème et peut être la plus délicate, consiste en la
modélisation du phénomène observé. On identifie tous les facteurs internes et ex-
ternes qui influencent la description du phénomène ; dans le cas d’un phénomène
physique, on fait un inventaire des forces en présence. Cet effort conduit en
général à des systèmes d’équations complexes comprennant un grand nombre
de paramàtres inconnus. Il faut alors simplifier certaines composantes et négliger
les moins importantes. On fait alors une première erreur dite de modélisation.
Du fait des difficultés de mise en équation des phénomènes physiques, on in-
troduit un modèle qui décrit le mieux son influence tout en demeurant une
approximation de la réalité.

23
Exemple Le pendule pesant idéal
On considère un pendule pesant fait d’une boule de masse m suspendu en A
à un fil OA de masse négligeable et de longueur l et fixé en O. On se propose
de décrire
le mouvement de A. On travail dans un repère orthonormé direct
0 , i , j où le poids du point A est p~ = −mg~j. Soit θ la mesure de l’angle
~ ~
−→

orienté −~j , OA .
En négligeant le frottement de l’air sur la boule, le point A est soumis à trois
forces vérifiant d’après la loi de Newton, T~ + p~ + F~N = ~0 où :
– T~ est la tension du fil
– p~ est le poids du point A
−→
2
– F~N = m d dt OA est la force de Newton
2

Par projection orthogonale sur la tangente en A à la trajectoire de A on a :

d2 θ
ml 2 + mg sin θ = 0 d’où
dt
d2 θ g
+ sin θ = 0 (2.1)
dt2 l
Cette équation ne donne qu’une approximation du mouvement réel de la
boule. Malgré cela sa résolution est loin d’être évidente

2.2.2 Erreurs de Troncature

Elles constituent une catégorie très importante surtout liée à l’utilisation de
développement de Taylor. Nous reviendrons plus en détail sur ce type d’erreurs
quand le besoin se fera sentir dans le déroulement de ce cours.
Nous devons simplement faire remarquer qu’il est très important de ne pas
confondre les erreurs de troncature avec la troncature utilisée pour la représentation
des nombres sur un ordinateur.
Exemple Le pendule pesant en faible amplitude
Le mouvement décrit par l’équation 2.1 ne fait aucune hypothèse sur les valeurs
de θ. Il faut observer que pour de petites valeurs de θ on peut remplacer sin θ
par son développement de Taylor à l’ordre 2n + 1 au voisinage de 0 ; c’est à dire
θ3 n θ
2n+1
sin θ = θ − . . . + (−1) + O(θ2n+2 )
3! (2n + 1)!
Pour n = 0, l’approximation obtenue est une équation différentielle ordinaire
linéaire d’ordre 2 à coefficients constants que nous savons résoudre :

d2 θ g
+ θ=0 (2.2)
dt2 l
La description donnée par l’équation 2.2 est entachée d’une erreur de tron-
cature par rapport à celle issue de l’équation différentielle non linéaire 2.1 de la
sous section précédente.

Exercice 10 Donner une approximation cubique du pendule pesant idéal.

24
Exercice 11 On considère la fonction f définie par

f (x) = 1.01349e5x − 5.2621e3x − 0.017326e2x + 0.83892ex − 1.9125

Question 1
Donner l’évaluation optimale de f (0.84) que peut offrir votre calculatrice.
Question 2
On évaluera f (0.8) à l’aide d’une arithmétique décimale tronquée à 4 chiffres
sachant que pour x = 0.8, ex = 2.225 :
a) En calculant d’abord les expressions suivantes e2x , e3x , e4x et e5x .
Indiquer l’erreur absolue et l’erreur relative.
b) En utilisant la méthode de Horner.
Indiquer l’erreur absolue et l’erreur relative.

25
Chapitre 3

Résolution des systèmes

linéaires

3.1 Introduction
La résolution d’un système linéaire de la forme Au = b est au centre de la
plupart des traitements de problèmes numériques, provenant de divers domaines
parmi lesquels on peut citer le calcul de reseau électrique.
On va considèrer le système modèle :


 15x1 − 2x2 − 6x3 = 300
−2x1 + 12x2 − 4x3 − x4 = 0


 −6x 1 − 4x 2 + 19x 3 − 9x 4 = 0
− x2 − 9x3 + 21x4 = 0


La solution de ce système en arithmétique décimale arrondie à 4 chiffres est

x1 = 26.55 x2 = 9.354 x3 = 13.25 x4 = 6.126

Le but de ce chapitre est de passer en revue quelques techniques ayant fait

leur preuve sur des systèmes du genre, même si ceux ci ont des centaines d’in-
connues.
Ces techniques sont classées en deux catégories : les méthodes directes et les
méthodes itératives.
Nous étudirons d’abord les méthodes directes c’est à dire celles qui con-
duisent en un nombre fini d’opérations élémentaires à la solution exacte lorsque
l’arithmétique utilisée est à précision infinie.
Nous aborderons ensuite les méthodes itératives consistant à rechercher la
solution du système comme limite d’une suite vectorielle. Ici la solution est
toujours entachée d’erreurs.
Dans la pratique nous utiliserons une arithmétique flottante, il faut donc
prendre en compte l’influence de l’erreur d’arrondi.

26
3.2 Méthodes directes de résolution de systèmes
linéaires Au = b
3.2.1 Systèmes triangulaires
Ils sont de deux types selon que le matrice A est une matrice triangulaure
supérieure ou une matrice triangulaure inférieure. Dans ces cas on utilise de
simples techniques de substitution.

Substitution par la remontée

Le système est de la forme


 a11 x1 + a12 x2 + . . . ... + a1n xn = b1


 a22 x2 + . . . ... + a2n xn = b2
.. ..


 ..


 . . .
aii xi + ... + ain xn = bi

 .. .. ..
. . .




an−1,n−1 xn−1 + an−1,n xn = bn−1




ann xn = bn


Le système est inversible si pour tout i = 1, 2, . . . , n, aii 6= 0.

Ainsi la procédure de remontée démarre par
bn
xn =
ann
et par une reccurence regressive on a
 
n
1  X
xi = bi − aij xj  i = n − 1, n − 2, . . . , 1
aii j=i+1

D’où l’algorithme :
real array (aij )n×n , (bi )n , (xi )n
integer i, j, n
real sum
xn ← bn /ann
for i = n − 1 to 1 step -1 do
begin
sum ← bi
for j = i + 1 to n do
begin
sum ← sum − aij xj
end
xi ← sum/aii
end

27
Le nombre de couples (multiplication - addition) que necessite cet algorithme
2
est de l’ordre de n2

Substitution par la descente

Le système est de la forme


 a11 x1 = b1


 a21 x1 + a22 x2 = b2
.. ..


 ..


 . . .
ai1 x1 + ... + aii xi = bi

 .. .. ..
. . .




a x + ... + ... + an−1,n−1 xn−1 = bn−1


 n−1,1 1


an1 x1 + ... + ... + an,n−1 xn−1 + ann xn = bn

Le système est inversible si pour tout i = 1, 2, . . . , n, aii 6= 0.

Ainsi la procédure de descente démarre par
b1
x1 =
a11
et par une reccurence progressive on a
 
i−1
1  X
xi = bi − aij xj  i = 2, 3, . . . , n
aii j=1

D’où l’algorithme :

real array (aij )n×n , (bi )n , (xi )n

integer i, j, n
real sum
x1 ← b1 /a11
for i = 2 to n do
begin
sum ← bi
for j = 1 to i − 1 do
begin
sum ← sum − aij xj
end
xi ← sum/aii
end

Il necessite le même nombre d’opérations que l’algorithme précédent.

28
Exemple
Résoudre le système

6x1 = 12
3x1 +6x2 = −12
4x1 −2x2 +7x3 = 14
5x1 −3x2 +9x3 +21x4 = −2

On a
x1 = 12/6 = 2 x1 = 2
6x2 = −12 − 3x1 = −18 x2 = −18/6 = −3
⇔
7x3 = 14 − (4x1 − 2x2 ) 7x3 = 14 − 14 = 0
21x4 = −2 − (5x1 − 3x2 + 9x3 ) 21x4 = −2 − (5x1 − 3x2 + 9x3 )

x1 = 2 x1 = 2
x2 = −3 x2 = −3
⇔ ⇔
x3 = 0 x3 = 0
21x4 = −2 − 19 = −21 x4 = −21/21 = −1

Exercice 12 Résoudre le système

6x1 −2x2 +2x3 +4x4 = 16

−4x2 +2x3 +2x4 = −6
+2x3 −5x4 = −9
−3x4 = −3

3.2.2 Elimination naı̈ve de Gauss

Le processus d’élimination naı̈ve de Gauss comporte deux étapes :
– la triangularisation qui consiste à transformer le système initial en un
système triangulaire supérieur
– la subtitution par la remontée pour résoudre le système triangulaire supérieur
obtenu

3.2.3 Triangularisation
Elle est obtenue en appliquant des transformations élémentaires successives
de Gauss à la matrice initiale. Commençons d’abord par nous faire une idée
précise d’une transformation de Gauss.

Transformations de Gauss
Elles ont pour fonction d’annuler les composantes d’un vecteur à partir d’un
certain rang. Soit n, k ∈ NI ∗ avec k < n. On considère un vecteur x ∈ K I n tel
que
xT = [x1 , x2 , . . . , xk , xk+1 . . . , xn ] , xk 6= 0

29
En posant
xi
αi = pour i = k + 1, . . . , n
xk

αT = [0, 0, . . . , 0, αk+1 . . . , αn ]
et  
1 ... 0 0 0 ... 0
 .. .. .. .. .. .. .. 
 .
 . . . . . . 
 0 ... 1 0 0 ... 0 
T
 
Mk = I − αek = 
 0 ... 0 1 0 ... 0  
 0
 ... 0 −αk+1 1 ... 0  
 . .. .. .. . . .. 
 .. . . . . . 
0 ... 0 −αn 0 ... 1
on a
T
Mk x = x − xk α ⇒ (Mk x) = [x1 , x2 , . . . , xk , 0 . . . , 0]
T
Exercice 13 Soit α = [0, . . . , 0, αk+1 , . . . , αn ] . Montrer que
si Mk = I − αeTk alors Mk−1 = I + αeTk
Exercice 14 Soit M = Mk . . . M1 où Mi = I − α(i) eTi avec
T
α(i) = [0, . . . , 0, −li+1,i , . . . , −lni ]
Montrer que eTj α(i) = 0 pour j < i.
Montrer que M est une matrice triangulaire inférieure de diagonale unité.
Plus explicitement on a
 
1
 l21 1 
 
 l31 l32 1 
 ..
 
.. 
 . . 
M =  
l
 k1 l k2 lk3 . . . lkk−1 1 

 lk+1,1 lk+1,2 lk+13 . . . lk+1,k 1 
 
 . .
 .. ..


ln1 ln2 . . . lnk 0 ... 1
Montrer que
k
X
M −1 = (I + α(1) eT1 ) . . . (I + α(k) eTk ) = I + α(i) eTi
i=1

Exercice 15 Déterminer la transformation I 3×3 telle que

de Gauss M ∈ R
   
2 2
M 3  =  7 
4 8

30
Théorème 3.2.1 On suppose que A ∈ K I m×n et que pour k < min(m, n) on a
déterminé les transformations de Gauss M1 , . . . , Mk−1 ∈ KI m×m telles que
" #
(k−1) (k−1)
(k−1) A11 A12 k−1
A = Mk−1 . . . M1 A = (k−1)
0 A22 m −k+1
k−1 n−k+1
(k−1)
où A11 est une matrice triangulaire supérieure. Si
 
(k−1) (k−1)
akk . . . akn
(k−1)  . .. 
A22 = .. .


(k−1) (k−1)
amk ... amn

(k−1)
et akk 6= 0 alors les multiplicateurs
(k−1)
aik
lik = (k−1)
i = k + 1, . . . , m
akk
T
sont définis et en posant Mk = I −α(k) eTk avec α(k) = (0, . . . , 0, lk+1,k , . . . , lmk )
on a
" #
(k) (k)
A A k
A(k) = Mk A(k−1) = 11 12
(k)
0 A22 m −k
k n−k
(k)
avec A11 matrice triangulaire supérieure.

Preuve
Elle est une utilisation immédiate des propriétés de la transformation de Gauss.
(Vérification laissée en exercice)

Théorème 3.2.2 Soit A ∈ K I n×n une matrice non singulière. On suppose que
l’on a déterminé les transformations de Gauss succeccives M1 , . . . , Mn−1 telles
que
U = A(n−1) = Mn−1 A(n−2) = Mn−1 . . . M1 A
soit une matrice triangulaire supérieure. Alors on a

Ax = b ⇔ U x = (Mn−1 . . . M1 ) b

On peut alors utiliser l’algorithme de substitution par la remontée

En résumé l’algorithme de triangularisation ou d’élimination de Gauss peut

s’écrire :

31
real array (aij )n×n , (bi )n
integer i, j, k, n
for k = to n − 1 do
begin
for i = k + 1 to n do
begin
for j = k to n do
begin
aij ← aij − (aik /akk )akj
end
bi ← bi − (aik /akk )bk
end
end

Remarque :

Dans cet algorithme nous observons que :

1) le multiplicateur aakk
ik
ne depend pas de j donc peut sortir de la boucle j
et on obtient

real array (aij )n×n , (bi )n

integer i, j, k, n
real xmult
for k = to n − 1 do
begin
for i = k + 1 to n do
begin
xmult ← aik /akk
aik ← xmult
for j = k to n do
begin
aij ← aij − (xmult)akj
end
bi ← bi − (xmult)bk
end
end

2) les termes aij pour i > j qui devaient être mis à 0 n’ont pas été mis à jour,
opération complètement superflue. L’espace occupé par ces élements pourra être
utilisé dans l’algorithme de factorisation LU

3.2.4 Factorisation LU
Théorème 3.2.3 Une matrice régulière A ∈ K I n×n possède une factorisation
A = LU où L est triangulaire inférieure à diagonale unité et U est triangulaire

32
supérieure, si et seulement si toutes les sous matrices principales A(k) de A sont
régulière.

Preuve

Si A = LU , nous pouvons avoir une décomposition par bloc des trois matrices

k A11 A12 L11 0 U11 U12 k
= .
(n − k) A21 A22 L21 L22 0 U22 (n − k)
k (n − k) k (n − k) k (n − k)

où A11 est la sous matrice principale d’ordre k de A c’est à dire A(k) . Le produit
par bloc donne

A11 = L11 U11 et det A11 = det L11 det U11

L11 est triangulaire inférieure à diagonale unité donc det L11 = 1.

U étantQune matrice triangulaire supérieure régulière on a
k
det U11 = i=1 U (i, i) 6= 0.
Ainsi
det A(k) = det A11 = det U11 6= 0

Réciproquement si pour tout k det A(k) 6= 0 on va montrer par reccurence

que l’algorithme naı̈f de Gauss s’applique.
En effet pour k = 1 la sous matrice est le scalaire non nul A11 . Donc l’etape
1 de l’algorthme naı̈f de Gauss est réalisable.
Supposons que nous ayons réalisé l’étape k − 1 de l’algorithme naı̈ de Gauss.
Alors on
A = A1 = L1 L2 . . . Lk−1 A(k) = RA(k)
où R = L1 L2 . . . Lk−1 est une matrice triangulaire inférieure à diagonale unité.
Posons B = A(k) Par une decomposition par bloc de dimension k et n − k on a

A11 A12 R11 0 B11 B12
= .
A21 A22 R21 I B21 B22

et A11 = R11 B11 avec

det B11 = det A11 = det A(k) puisque det R11 = 1

Etant donnée la structure de B = A(k) alors B11 est est une matrice trian-
(k)
gulaire supérieure avec tous les termes diagonaux non nuls. En particulier Akk
eme
est non nuls et peut être choisi comme le k pivot.

Théorème 3.2.4 Si une matrice régulière A ∈ K I n×n possède une factorisation

A = LU où L est triangulaire inférieure à diagonale unité et U est triangulaire
supérieure, alors cette factorisation est unique.

33
Preuve

Supposons qu’il existe deux décompositions A = L1 U1 = L2 U2 où toutes les

matrices sont régulières. Alors X = L−1 2 L1 = U2 U1
−1
est à la fois triangulaire
inférieure à diagonale unité et triangulaire supérieure. X est nécessairement une
matrice unité. Ainsi L1 = L2 et U1 = U2 d’où l’unicité de la décomposition.
Exemple
Résoudre par l’élimination naı̈ve de Gauss puis déduire la factorisation de la
matrice A du système
    
6 −2 2 4 x1 16
 12 −8 6 10   x2  =  26 
   

 3 −13 9 3   x3   −19 
−6 4 1 −18 x4 −34

On part du tableau :
 
E1 : 6 −2 2 4 16
E2 : 
 12 −8 6 10 26 

E3 :  3 −13 9 3 −19 
E4 : −6 4 1 −18 −34

A l’étape 1 on a
 
E1 : 6 −2 2 4 16
E2 ← E2 − 12
6 E1 : 
 0 −4 2 2 −6 
E3 ← E3 − 36 E1 :  0 −12 8 1 −27 
E4 ← E4 − −6
6 E1 : 0 2 3 −14 −18

A l’étape 2 on a
 
E1 : 6 −2 2 4 16
E2 : 
 0 −4 2 2 −6 
E3 ← E3 − −12

−4 E2 :  0 0 2 −5 −9 
2
E4 ← E4 − −4 E2 : 0 0 4 −13 −21

A l’étape 3 on a
 
E1 : 6 −2 2 4 16
E2 : 
 0 −4 2 2 −6 

E3 :  0 0 2 −5 −9 
E4 ← E4 − 42 E3 : 0 0 0 −3 −3

D’où par la remontée on a

x4 = −3/(−3) = 1 x4 = −3/(−3) = 1
2x3 = −9 − (−5x4 ) = −4 x3 = −4/(2) = −2
⇔
−4x2 = −6 − (2x3 + 2x4 ) −4x2 = −6 − (−4 + 2) = −4
6x1 = 16 − (−2x2 + 2x3 + 4x4 ) 6x1 = 16 − (−2x2 + 2x3 + 4x4 )

34
x4 = −3/(−3) = 1 x4 = −3/(−3) = 1
x3 = −4/(2) = −2 x3 = −4/(2) = −2
⇔ ⇔
x2 = −4/(−4) = 1 x2 = −4/(−4) = 1
6x1 = 16 − (−2 − 4 + 4) = 18 x1 = 18/6 = 3
On a également la factorisation A = LU avec
   
6 −2 2 4 1 0 0 0
 0 −4 2
 et L =  12
2  1 0 0 

U = 
 0 0 2 −5  
2 3 1 0 
0 0 0 −3 −1 − 21 2 1

Exercice 16 En utilisant l’élimination naı̈ve de Gauss, calculer la factorisation

LU de A
0 13
 
  1 0
3 0 3  0 1 3 −1 
a) A =  0 −1 3  b) A =   3 −3 0

6 
1 −3 0
0 2 4 −6

Exercice 17 On considère
 
2 2 1
A= 1 1 1 
3 2 1

a) Montrer que A n’admet pas une factorisation LU avec L une matrice

triangulaire inférieure à diagonale unité et U une matrice triangulaire supérieure
b)Permuter les lignes pour que la factorisation soit possible.

I n×n les matrices

Exercice 18 On considère dans R

N (y, k) = I + yeTk avec In

y ∈R

appelée transformation de Gauss-Jordan.

−1
1) En supposant que N (y, k) existe, donner en une expression.
n
2) Etant donné x ∈ R I n pour que
I , à quelle condition peut on trouvé y ∈ R

N (y, k)x = ek

3)Donner un algorithme utilisant les transformations de Gauss Jordan pour

surcharger A par A−1 . Donner des conditions sur A pour que cet algorithme
marche.

Exercice 19 Soit A ∈ RI n×n définie par


 1 si i = j ou j = n
aij = −1 si i > j
0 autrement


Montrer que A admet une décomposition LU avec |lij | ≤ 1 et unn = 2n−1 .

35
Exercice 20 Montrer que si AT ∈ R I n×n est à diagonale dominante alors A
admet une décomposition LU avec |lij | ≤ 1.

Exercice 21 Résoudre par élimination de Gauss les systèmes :


 3x1 +4x2 +3x3 = 10
(a) x1 +5x2 −x3 =7
6x1 +3x2 +7x3 = 15



 3x1 +2x2 −5x3 = 0
(b) 2x1 −3x2 x3 = 0
x1 +4x2 −x3 = 4

    
1 −1 2 1 x1 1
 3 2 1 4   x2   1 
(c) 
 5
 = 
8 6 3   x3   1 
4 2 5 3 x4 −1

3.2.5 Elimination de Gauss avec stratégie de pivot

Il s’agit ici de modifier l’approche naı̈ve afin de prendre en compte une
stratégie de choix des pivots. Il faut par exemple commencer par résoudre le cas
(k−1)
d’inopérabilité où l’étape k est bloqué par Akk = 0. On peut avoir satisfaction,
soit à travers une permutation de la ligne k avec une ligne en dessous, c’est à
dire une permutation d’équations ; soit à travers une permutation d’inconnues.
Nous allons regarder essentiellement une stratégie dite de pivot partiel
normalisé. Cette stratégie prend en compte des critères de stabilité du proces-
sus de triangularisation.
étape k = 1
On normalise chaque équation du système, en recherchant tout simplement
le vecteur de normalisation s(1) = [s1 , s2 , . . . , sn ] avec

si = max |aij |
1≤j≤n

On regarde ensuite les rapports

|ai1 |
i = 1, 2, . . . , n
si

On choisit l’indice j correspondant à la première occurence de la plus grande

valeur de cet ensemble. On procede alors à la permutation des lignes (équations)
indexées 1 et j. L’étape se termine par l’application de la transformation de
Gauss M1 . En résumé on a
" #
(1) (1)
A A
A(0) = A → A g (0) = P A(0) → A(1) = M P A(0) =
1 1 1
11 12
(1)
0 A22

36
(1) (0)
avec A11 = A11 ∈ K I le pivot 1.
g
Pour l’étape 2 on applique l’étape 1 à la sous matrice
 
(1) (1)
a22 . . . a2n
(1)  . ..  (n−1)×(n−1)
 ..
A22 =  ∈K
.  I
(1) (1)
an2 . . . ann

. On construit le vecteur de normalisation s(2) = [s2 , . . . , sn ] avec

(1)
si = max aij
2≤j≤n

On regarde ensuite les rapports

 (1) 
 ai2 
i = 2, . . . , n
 si 

On choisit l’indice j correspondant à la première occurence de la plus grande

valeur de cet ensemble. On procède alors à la permutation des lignes (équations)
I (n−1)×(n−1) cette permutation. On obtient
f2 ∈ K
indexées 2 et j. Soit P
 
(1) (1)
 a22 . . . a2n 
g g

A22 = P2 A22 =  ... ..

(1) (1)
I (n−1)×(n−1)
∈K
g  
.
f
 
(1) (1)
an2 . . . ann
g g

. En posant P2 = diag(1, P
f2 ) on a
" #
(2) (2)
(2) (1) A11 A12
A = M 2 P2 A = (2)
0 A22
(2) (2)
où A11 ∈ K I 2×2 est triangulaire supérieure et A22 ∈ K I (n−2)×(n−2) ; M2 étant
une transformation de Gauss.
En supposant ce processus réalisé jusqu’à l’étape k − 1, c’est à dire exis-
tence de transformation de Gauss M1 , . . . , Mk−1 ∈ K I n×n et de permutations
n×n
P1 , . . . , Pk−1 ∈ K
I telles que
" #
(k−1) (k−1)
A A
A(k−1) = Mk−1 Pk−1 . . . M1 P1 A = 11 12
(k−1)
0 A22
(k−1) (k−1)
avec A11 I (k−1)×(k−1) triangulaire supérieure et A22
∈K ∈K I (n−k+1)×(n−k+1)
étape k
On construit le vecteur de normalisation s(k) = [sk , . . . , sn ] avec
(k−1)
si = max aij
k≤j≤n

37
On regarde ensuite les rapports
 (k−1) 
 aik 
i = k, . . . , n
 si 

On choisit l’indice j correspondant à la première occurence de la plus grande

valeur de cet ensemble. On procède alors à la permutation des lignes (équations)
fk ∈ K
indexées k et j. Soit P I (n−k+1)×(n−k+1) cette permutation. On obtient
 
(k−1) (k−1)
 a22 . . . a2n
g g

(k) fk A(k−1) =  .. .. I (n−k+1)×(n−k+1)
A22 = P ∈K
g 
22  . .
 
(k−1) (k−1)
an2 . . . ann
g g

. En posant Pk = diag(Ik−1 , P
fk ) on a
" #
(k) (k)
(k) (k−1) A11 A12
A = M k Pk A = (k)
0 A22

I (k−1)×(k−1) et Mk = In − α(k) eTk avec

Ik−1 ∈ K
h iT
α(k) = 0, . . . , e
lk+1,k . . . , e
ln,k

où
(k−1)
a
g
li,k = ik
e i = k + 1, . . . , n
(k−1)
akk
g

En résumé, l’élimination de Gauss avec une stratégie de pivot partiel nor-

malisé consiste à construire des transformations de Gauss M1 , . . . , Mn−1 et des
permutations P1 , . . . , Pn−1 telles que

I n×n
Mn−1 Pn−1 . . . M1 P1 A = U ∈ K

où U est une matrice triangulaire supérieure.

Posons P = Pn−1 . . . P1 alors
−1
L = P (Mn−1 Pn−1 . . . M1 P1 )

est une matrice triangulaire inférieure est l’on a la factorisation

P A = LU

Remarque :

Au lieu de se lancer dans des permutations effectives et laborieuses des lignes,

dans la pratique, on se contentera d’associer aux équations, un vecteur d’index

38
T T
I n que l’on initialisera à ` = [1, 2, . . . , n] . Les permutations
` = [`1 , . . . , `n ] ∈ N
de lignes se résumeront à des permutations des éléments de `.
Par ailleurs en pratique les vecteurs de normalisation s ne change pas sensi-
blement.
On a alors l’algorithme :

real array (aij )n×n ,(bi )n , (`i )n

real array (si )n
integer i, k, n
real r , rmax, smax , xmult
for i = 1 to n do
begin
`←i
smax ← 0
for j = 1 to n do
begin
smax ← max (smax, |aij |)
end
si ← smax
end
for k = 1 to n − 1 do
begin
rmax ← 0
for i = k to n do
begin
r ← |aì ,k /sì |
if (r > rmax) then
begin
rmax ← r
j←i
end
end
`j ↔ `k
for i = k + 1 to n do
begin
xmult ← aì ,k /a`k ,k
aì ,k ← xmult
for j = k + 1 to n do
begin
aì j ← aì j − (xmult)a`k j
end
bì ← bì − (xmult)b`k
end
end

39
Théorème 3.2.5 Dans l’algorithme de Gauss appliqué à une matrice A ∈
3
I n×n , la phase de triangularisation a un coût de l’ordre de n3
K

Remarque :
(k−1)
Il faut observer que dans le processus de triangularisation, plus un pivot akk
(k−1)
aik
est petit en valeur absolue, plus les coefficients multiplicateurs (k−1) sont
akk
grands ce qui amplifira l’effet des erreurs d’arrondi inhérentes à l’utilisation
d’une arithmétique à précision finie.
Ainsi l’objectif d’une stratégie de pivot est de réduire le plus possible les
(k−1)
aik
facteurs d’amplification que sont les coefficients (k−1) . Il faudra tout au moins
akk
assurer la condition
(k−1)
aik
(k−1)
≤ 1pour i > k
akk
Cette condition est réalisée avec la stratégie de pivot partiel normalisé.
On peut aussi envisager une stratégie dite de pivot total normalisé qui est
plus stable du point de vue numérique.
Remarque :
L’inversion de matrice peut se faire directement en utilisant la méthode de
Gauss.
Exemple
Déterminons l’inverse de la matrice
 
4 −1 7
A =  −1 4 0 
−7 0 4

Il s’agit de résoudre l’équation AX = I avec

   
1 0 0 x1 y1 z1
I =  0 1 0  et X =  x2 y2 z2 
0 0 1 x3 y3 z3

On part du tableau :
 
E1 : 4 −1 7 1 0 0
E2 :  −1 4 0 0 1 0 
E3 : −7 0 4 0 0 1

A l’étape
1 les coefficients normalisés sont dans l’ordre
4 1 7
, , On procède à une permutation de ligne :
7 4 4
 
E1 ← E3 : −7 0 4 0 0 1
E2 :  −1 4 0 0 1 0 
E3 ← E1 : 4 −1 7 1 0 0

40
On a alors les facteurs
multiplicateurs
−1 1 4 4
= , =− et on a
−7 7 −7 7
 
E1 : −7 0 4 0 0 1
E2 ← E2 − 17 E1 :  0 4 − 47 0 1 − 17 
E3 ← E3 + 47 E1 : 0 −1 657 1 0 4
7

A l’étape2 les coefficients normalisés sont dans l’ordre

4 7
= 1,
4 65
On a alors le facteur multiplicateur
−1 1
= − et on a
4 4
 
E1 : −7 0 4 0 0 1
E2 :  0 4 − 74 0 1 − 71 
E3 ← E3 + 14 E2 : 0 0 64
7 1 1
4
15
28

Intervient alors un processus de remontée :

 64  7 7
 7 x3 = 1  x3 = 1 64 = 64
64 1 1 7 7
7 y3 = 4 ⇔ y3 = 4 64 = 256
 64 15 15 7 15
z = z3 = 28 64 = 256

7 3 28

= 0 + 47 x3 = 74 64
7 1 1
 
 4x2 = 16  x2 = 64
4 4 7 65 65
4y2 = 1 + 7 y3 = 1 + 7 256 = 64 ⇔ y2 = 256
1 4 1 4 15 7 7
4z2 = − 7 + 7 z3 = − 7 + 7 256 = − 64 z2 = − 256
 

7 7 1
 
 −7x1 = 0 − 4x3 = −4 64 = − 16  x1 = 16
7 7 1
−7y1 = 0 − 4y3 = −4 256 = − 64 ⇔ y1 = 64
15 49 7
−7z1 = 1 − 4z3 = 1 − 4 256 = 64 z1 = − 64
 

D’où
1 1 7
 
16 64 − 64
A−1 =  1
64
65
256
7 
− 256
7 7 15
64 256 256

Exercice 22 On considére la matrice

 
2 −1 0
A =  −1 2 −1 
0 −1 2

Calculer l’inverse A−1 de A en utilisant la méthode de Gauss.

41
3.2.6 Système linéaires spéciaux : factorisations de Cholesky,
LDM t ...
T
Factorisation L-D-M
Théorème 3.2.6 Si toutes les sous matrices principale de la matrice régulière
A∈K I n×n sont régulières alors il existe deux matrices triangulaires inférieures
à diagonale unité L , M et une matrice diagonale D = diag(d1 , . . . , dn ) telles
que A = LDM T

Preuve

D’après le théorème 3.2.3, A admet une décomposition A = LU avec L une

matrice triangulaire inférieure à diagonale unité et U une matrice triangulaire
supérieure.
En définissant D = diag(d1 , . . . , dn ) telle que di = uii pour tout i = 1, . . . , n,
on a D régulière et la matrice M définie par M T = D−1 U est triangulaire
inférieure à diagonale unité. Ainsi

A = LU = LD(D−1 U ) = LDM T

La décomposition A = LDM T est obtenue de façon plus interessante en

regardant directement terme à terme l’équation matricielle.
En effet on a pour tout k
k−1
X
akk = lkp dp mkp + dk
p=1

et pour tout i > k

Pk−1
aik = p=1 lip dp mkp + lik dk
Pk−1
aki = p=1 lkp dp mip + dk mik

Ainsi
Pour tout k = 1, . . . , n
k−1
X
dk = akk − lkp dp mkp
p=1

et pour i = k + 1, . . . , n
Pk−1
lik = aik − p=1 lip dp mkp /dk
Pk−1
mik = aki − p=1 lkp dp mip /dk

n3
L’algorithme associé a un coût de l’ordre de 3 opérations.

42
Exercice 23 Ecrire l’algorithme de la décomposition A = LDM T en sur-
chargeant aij par lij si i > j et par mij si i < j.

Exemple
On considère la matrice  
2 −1 2
A= 2 −3 3 
6 −1 8
1. Trouver la factorisation A = LDV où L est triangulaire inférieure, V
triangulaire supérieure à diagonales unité et D diagonale.
T
2. Utiliser cette décomposition pour résoudre Ax = b avec b = [−2, −5, 0]
1 ) On commence par une factorisation LU classique en utilisant le processus
naı̈f de triangularisation de Gauss.
On part du tableau
   
E1 : 2 −1 2 1 0 0
E2 :  2 −3 3  et on construit  22 1 0 
6
E3 : 6 −1 8 2 0 1

Etape 1 :
   
E1 : 2 −1 2 1 0 0
E2 ← E2 − 22 E1 :  0 −2 1  et on construit  1 1 0 
2
E3 ← E3 − 62 E1 : 0 2 2 3 −2 1

Etape 2 :
   
E1 : 2 −1 2 1 0 0
E2 :  0 −2 1  et on construit  1 1 0 
2
E3 ← E3 − −2 E2 : 0 0 3 3 −1 1

D’où    
1 0 0 2 −1 2
L= 1 1 0  et U =  0 −2 1 
3 −1 1 0 0 3
Soit D = diag(2 , −2 , 3). Calculons la matrice V telle que U = DV . Puisque
D est inversible, on a V = D−1 U . Ainsi

1 − 21
 
1
V = 0 1 − 21 
0 0 1
2 ) Résolution de Ax = b se ramène celle de LDV x = b. On procède donc
par la séquence :
Ly = b → Dz = y → V x = z

43
T
L’équation Ly = b, avec b = [−2, −5, 0] , donne

y1 = −2
y2 = −5 − y1 = −5 + 2 = −3
y3 = 0 − 3y1 + y2 = 6 − 3 = 3
T
et y = [−2, −3, 3]
L’équation Dz = y donne z = D−1 y. Comme D = diag(2 , −2 , 3) alors on
a T T
−2 −3 3 3
z= , , = −1, , 1
2 −2 3 2
Enfin l’équation V x = z donne :

x3 = 1
x2 = 32 + 12 x3 = 23 + 12 = 2
x1 = −1 + 12 x2 − x3 = −1 + 21 (2) − 1 = −1
T
D’où x = [−1, 2, 1]
Remarque :
L’utilisation de la factorisation de A dans la résolution du système linéaire
Ax = b peut être rentable lorsque l’on opère pour toute une série de b.

Exercice 24 On considère
   
−2 1 −2 1
A =  −4 3 −3  et b =  4 
2 2 4 4

1. Trouver la factorisation A = LDU où L est triangulaire inférieure, U

triangulaire supérieure à diagonales unité et D diagonale.
2. Utiliser cette décomposition pour résoudre Ax = b

T
Factorisation L-D-L
Théorème 3.2.7 Si A = LDM T est la décomposition LDM T d’une matrice
régulière symétrique A alors L = M et on a A = LDLT .

Preuve
−1
La matrice M −1 A(M T ) = M −1 LD est symétrique et triangulaire inférieieure.
Elle est donc diagonale. Puisque D est régulière alors M −1 L est aussi diagonale.
Mais M −1 L est triangulaire à diagonale unité. On a donc M −1 L = I.
Remarque :
L’agorithme de factorisation A = LDLT d’une matrice symétrique A ∈ K I n×n a
3
un coût de l’ordre de n6

44
Factorisation de Cholesky
Théorème 3.2.8 Si une matrice A ∈ R I n×n est symétrique définie positive
I n×n à termes diagonaux
alors il existe une matrice triangulaire inférieure G ∈ R
T
positifs telle que A = GG .

Preuve

D’après le théorème 3.2.7 il existe une matrice triangulaire inférieure à diag-

onale unité L ∈ R I n×n et une matrice diagonale D = diag(d1 , . . . , dn ) telle que
T
A = LDL .
Comme A est definie positive et L régulière alors pour tout k = 1, . . . , n il
existe xk ∈ RI n tel que Lxk = ek et

o < (Axk |xk ) = (DLxk |Lxk ) = (Dek |ek ) = dk

√ √
La matrice G = Ldiag( d1 , . . . , dn ) est réelle, triangulaire inférieure à termes
diagonaux positifs et on a GGT = LDLT = A.
Pour calculer concrètement la factorisation de Cholesky il est plus interessant
de regarder directement les termes de l’équation matricielle A = GGT .
Ainsi pour i ≥ k on a
Xk
aik = gip gkp
p=1

après arrangement on a
r
Pk−1 2
gkk = akk − p=1 gkp
Pk−1
gik = aik − p=1 gip gkp /gkk

L’algorithme de décomposition de Cholesky pour une matrice réelle, symétrique,

definie positive A sera donné ici en calculant la matrice triangulaire inférieure
puis en surchargeant les termes aij par gij pour i ≥ k
On a :

real array (aij )n×n

integer i, j, k, n
real sum
for k = 1 to n do
begin
for j = 1 to k − 1 do
begin
akk ← akk − akj ∗ akj
end
√
akk ← akk
for i = k + 1 to n do
begin

45
for j = 1 to k − 1 do
begin
aik ← aik − aij ∗ akj
end
aik ← aik /akk
end
end

Exemple
Calculer la factorisation A = LDLT puis la factorisation de Cholesky A = GGT
de la matrice  
6 2 1 −1
 2 4 1 0 
A= 
 1 1 4 −1 
1 0 −1 3
Factorisation LDLT
Commençons par une factorisation A = LU classique avec L à diagonale
unité, en utilisant le processus naı̈f de triangularisation de Gauss.
On part du tableau
   
E1 : 6 2 1 −1 1 0 0 0
 2
E2 :  2 4 1 0  et via  61 1 0 0 

E3 :  1 1 4 −1  
6 0 1 0 
−1
E4 : −1 0 −1 3 6 0 0 1

Etape 1 :
   
E1 : 6 2 1 −1 1 0 0 0
E2 ← E2 − 26 E1 : 
 0
10
3
2
3
1
3
 
 et via 
1
3 1 0 0 

E3 ← E3 − 16 E1 :  0 2
3
23
6 − 56   1
6
2
10 1 0 
E4 ← E4 − −1
6 E1 : 0 1
3 − 56 17
6 − 61 1
10 0 1

Etape 2 :
   
E1 : 6 2 1 −1 1 0 0 0
10 2 1 1
E2 : 
 0 3 3 3
 
 et via  3 1 0 0 
2 37 9 1 1

E3 ← E3 − 10 E2 :  0 0 10 − 10  
6 5 1 0 
1 27 84 −9
E4 ← E4 − 10 E2 : 0 0 − 30 30 − 61 1
10 37 1

Etape 3 :
   
E1 : 6 2 1 −1 1 0 0 0
10 2 1 1
E2 : 
 0 3 3 3
 
 et via  3 1 0 0 

37 9 1 1
E3 :  0 0 10 − 10  
6 5 1 0 
E4 ← E4 − −9 E
37 3 : 0 0 0 191
74 − 16 1
10
9
− 37 1

46
Ainsi
   
1 0 0 0 6 2 1 −1
1 10 2 1

3 1 0 0   0
3 3 3

L= 1 1
 et U =  37 9


6 5 1 0   0 0 10 − 10 
1 1 9 191
−6 10 − 37 1 0 0 0 74

10 37 191
Pour D = diag 6 , , , on obtient A = LDLT .
3 10 74
Factorisation LDLT √
T
On peut déduire r de!Cholesky A = GG avec G = L D où
r la factorisation
√ √
r
10 37 191
D = diag 6, , ,
3 10 74
On a donc
 √ 
6 0 0 0
 √6
q
10
0 0 

 √3 3

G= q q 
6 1 10 37
0

 √6 5q 3 q 10 q
 

6 1 10 9 37 191
− 6 10 3 − 37 10 74

Exercice 25 Reprendre l’exercice de l’exemple ci-dessus en travaillant dans

une arithmétique flottante décimale tronquée à 4 chiffres.

Exercice 26 Calculer la factorisation A = LDLT puis la factorisation de

Cholesky A = GGT de la matrice
 
4 1 −1 0
 1 3 −1 0 
A=  −1 −1

5 2 
0 0 2 4

Exercice 27 Soient α ∈ R I n−1×n−1 et w, ν ∈ R

I , β ∈R I n−1 . Montrer que si

α wT

A=
ν β
est définie positive alors il en est de même pour β − νwT /α

I n×n est symétrique, définie positive alors si

Exercice 28 Montrer que si A ∈ R
i 6= j, on a
(a +a )
|aij | ≤ ii 2 jj
√
|aij | ≤ aii ajj

47
3.2.7 Analyse d’erreur d’arrondi
Dans cette introduction nous n’allons pas nous focaliser sur l’analyse théorique
de l’effet de l’arithmétique finie qui fait largement état du conditionnement de la
matrice en jeu. Nous observerons ce phénomène à travers un exemple illustratif
laissant le reste en exercice.
Exemple
En utilisant une arithmétique décimale arrondie à 4 chiffres, résoudre par la
méthode de Gauss avec une stratégie de pivot partiel pondéré, le système


 15x1 − 2x2 − 6x3 = 300
−2x1 + 12x2 − 4x3 − x4 = 0


 −6x 1 − 4x 2 + 19x 3 − 9x 4 = 0
− x2 − 9x3 + 21x4 = 0


Déduire de ce processus une factorisation de A

1. Méthode de Gauss
On part du tableau :
 
E1 : 15 −2 −6 0 300
E2  −2 12 −4 −1
:  0 

E3 :  −6 −4 19 −9 0 
E4 : 0 −1 −9 21 0

A
l’étape 1 les coefficients normalisés sont dans
l’ordre
15 2 6
= 1, = 0.1666667 , = 0.6666667 , 0
15 12 9
On
a alors les facteurs multiplicateurs
−2 −6 0
= −0.1333333 , = −04 , =0
15 15 15
Puis on obtient :
 
E1 : 15 −2 −6 0 300
E2 ← E2 + 0.1333E1 :   0 11.73 −4.533 −1 39.49 

E3 ← E3 + 0.4E1 :  0 −4.8 16.6 −9 120 
E4 ← E4 + 0E1 : 0 −1 −9 21 0

A
l’étape 2 les coefficients normalisés sont dans l’ordre

11.73 4.8 1
= 1, = 0.2891566 , = 0.0476190
11.73 16.6 21
On
a alors les facteurs multiplicateurs
−4.8 −1
= −0.4092072 , = −0.0852515
11.73 11.73
Puis on obtient :
 
E1 : 15 −2 −6 0 300
E2 :  0 11.73 −4.533 −1 39.49 
E3 ← E3 + 0.4092E2 :  0 0 14.74 −9.409 136.2 
E4 ← E4 + 0.08525E2 : 0 0 −9.386 20.91 3.367

48
A
l’étape 3 les coefficients normalisés
sont dans l’ordre
14.74 9.386
= 1, = 0.4488761
14.74 20.91
On a alors le facteur multiplicateur
−9.386
= −0.6367707
14.74
Puis on obtient :
 
E1 : 15 −2 −6 0 300
 0 11.73 −4.533
E2 :  −1 39.49 
E3 :  0 0 14.74 −9.409 136.2 
E4 ← E4 + 0.6368E3 : 0 0 0 14.92 90.10
Par la remontée on a :
x4 = 90.10/14.92 = 6.039 x4 = 90.10/14.92 = 6.039
14.74x3 = 136.2 − (−9.409x4 ) = 193 x3 = 193/14.74 = 13.09
⇔
11.73x2 = 39.49 − (−4.533x3 − 1x4 ) 11.73x2 = 39.49 − (−4.533 ∗ 13.09 − 6.039) = 104.9
15x1 = 300 − (−2x2 − 62x3 ) 15x1 = 300 − (−2x2 − 6x3 )

x4 = 90.10/14.92 = 6.039 x4 = 90.10/14.92 = 6.039

x3 = 193/14.74 = 13.09 x3 = 193/14.74 = 13.09
⇔ ⇔
x2 = 104.9/11.73 = 8.943 x2 = 104.9/11.73 = 8.943
15x1 = 300 − (−2 ∗ 8.943 − 6 ∗ 13.09) = 396.4 x1 = 396.4/15 = 26.43
2. Factorisation :

– D’une part, il n’y a eu aucune permutation de ligne, donc on a la fac-

torisation A = LU avec
   
15 −2 −6 0 1 0 0 0
 0 11.73 −4.533 −1  et L =  −0.1333 1 0 0 

U = 
 0 0 14.74 −9.409   −0.4 −0.4092 1 0 
0 0 0 14.92 0 −0.08525 −0.6368 1

– D’autre part on a une factorisation A = LDM T avec

D = diag(15 , 11.73 , 14.74 , 14.92),
ce qui donne :
D−1 = diag(0.6667 , 0.08525 , 0.06784 , 0.06702)
puis on a
 
1 −0.1333 −0.4 0
0 1 −0.3864 −0.08525
M T = D−1 U = 
 

 0 0 1 −0.6383 
0 0 0 1
d’où  
1 0 0 0
 −0.1333 1 0 0 
M = ≈L
 −0.4 −0.3864 1 0 
0 −0.08525 −0.6383 1

49
ce qui est attendu car A symétrique.
– Comme les termes de D = diag(15 , 11.73 , 14.74 , 14.92)√sont positifs,
on a une factorisation de cholesky A = GGT avec G = L D où
√ √ √ √ √
D = diag 15 , 11.73 , 14.74 , 14.92 = diag(3.873 , 3.425 , 3.839 , 3.863)

Ainsi on a
 
3.873 0 0 0
 −0.5163 3.425 0 0 
G=
 −1.549

−1.402 3.839 0 
0 −0.2920 −2.445 3.863

Exercice 29 Soit u la borne supérieure de l’erreur d’arrondi dans la représentation

d’un nombre réel en arithmétique flottante. C’est à dire que pour α ∈ R I , sa
représentation
f l(α) = α(1 + ε)avec |ε| ≤ u
I n . Montrer par reccurence sur n que
Soit x, y ∈ R

f l(xT y) − xT y ≤ nu xT |y| + O(u2 )

sachant que si xT = (x1 , . . . , xn ), xT = (|x1 | , . . . , |xn |)

I m×n , B ∈ R
Exercice 30 1) Montrer que si A ∈ R I n×p on a

f l(AB) = AB + E avec |E| ≤ nu |A| |B| + O(u2 )

2) Montrer que si B est régulière alors

kf l(AB) − ABkF
≤ nuκF (B) + O(u2 )
kABkF

où
m X
X n
kAkF = a2ij pour tout I m×n
A(aij ) ∈ R
i=1 j=1

50
3.3 Méthodes itératives de résolution de systèmes
linéaires Au = b
3.3.1 Contexte général
Nous aborderons la résolution du système linéaire Ax = b avec un esprit
complètement différent de celui qui a prévalu lors de la mise en œuvre des
méthodes dites directes.
Il s’agit ici de construire une suite de vecteurs x(0) , x(1) , . . . , x(k) , . . . ap-
prochant la solution x du système. Le processus numérique est conçu de manière
à assurer la convergence de la suite vers la solution. Ce processus peut s’arrête
dès qu’une précision suffisante préalablement fixée est atteinte.
De façon générale un algorithme itératif pour le système linéaire Ax = b est
conçu en se donnant une matrice régulière Q ∈ K I n×n , puis, partant d’un vecteur
arbitraire x , d’engendrer la suite de vecteurs x(1) , . . . , x(k) , . . . par l’équation
(0)

reccurente :
Qx(k) = (Q − A)x(k−1) + b k = 1, 2, . . .
En supposant que la suite x(k) k converge vers x∗ , du fait de la continuité

des applications linéaires, le passage à la limite quand k → ∞ dans l’équation

Qx(k) = (Q − A)x(k−1) + b k = 1, 2, . . .

donne Qx∗ = (Q − A)x∗ + b ⇔ Ax∗ = b. Cette limite est donc la solution du

système initial.
L’algorithme général d’une méthode itérative est le suivant :

array (A)n×n ,(b)n ,(c)n , (y)n , (x)n

integer k, kmax
begin
x ← x(0)
for k = 1 to kmax do
begin
y←x
c ← (Q − A)x + b
Resolution de Qx = c
Ecrire k, x
if kx − yk < then
begin
Ecrire ”Convergence”
stop
end
end
Ecrire ”Nombre maximal d’iterations atteint”
end

51
Remarque :

Le choix de la matrice régulière Q doit prendre en compte le fait que le

système reccurent

Qx(k) = (Q − A)x(k−1) + b k = 1, 2, . . .

doit être facile à résoudre en x(k) lorsque x(k−1) est connu. Ce choix doit être
fait pour assurer une convergence éventuellement rapide du processus.

Théorème 3.3.1 Soit b ∈ R I n et A, Q ∈ RI n×n deux matrices régulières. Si le

rayon spectrale de I − Q A satisfait à la condition ρ(I − Q−1 A) < 1 alors
−1

la suite de vecteurs x(k) définie par Qx(k) = (Q − A)x(k−1) + b converge vers

x = A−1 b pour toute valeur initiale x(0) .

Preuve

Posons e(k) = x(k) − x l’erreur dans l’itération k. Des équations

Ax b
(k)
Qx = (Q − A)x(k−1) + b
on a
e(k) = x(k) − x = (I − Q−1 A)x(k−1) + Q−1 b − x
= (I − Q−1 A)x(k−1) + Q−1 Ax − x
= (I − Q−1 A)x(k−1) − (I − Q−1 A)x
= (I − Q−1 A)(x(k−1) − x)
= (I − Q−1 A)e(k−1)
k
= (I − Q−1 A) e(0)
Puisque ρ(I − Q−1 A) < 1, on sait que la suite matricielle
k
(I − Q−1 A) → 0

Par définition la matrice I −Q−1 A est dite la matrice itérative de la méthode.

Exemple
On considère le système linéaire
1 1 1
2 x2 + 3 x2 = 63
1 1 1
3 x2 + 4 x2 = 168

Soit A = (aij ) la matrice du système et b le second membre.

Soit D = diag(a11 , a22 ) la matrice diagonale associée.
On se propose de résoudre le système par la méthode itérative définie par
1
D−1 b − AX (k)

X (k+1) = X (k) + 10 pour k ≥ 0
X (0) = (0 , 0)

1. Calculer X (1) et X (2)

52
2. Déterminer la matrice itérative G de cette méthode puis calculer ses
valeurs propres.
3. Que peut-t-on dire de la converge de la méthode ?
On observe que :
1 1 1

A= 2 3 et b = 63
1 1 1
3 4 168

1 1
On a donc I = diag(1 , 1), Q = D = diag , et Q−1 = D−1 = diag(2 , 4)
2 4
1. Calcul de X (1) et X (2)
T
– X (0) = (0 , 0) alors
T
AX (0) = (0 , 0)
1 1 T

b − AX (0) = 63 , 168
1 T
D−1 2

b − AX (0) = 63 , 42
1 −1 1 1 T

10 D b − AX (0) = 315 , 420
1 T
X + 10 D−1
1 1
(0)

b − AX (0) = 315 , 420
1 1 T
X (1) = 315 , 420
T
1 1
– X (1) = , alors
315 420
1 5
T
AX (1) = 420 , 3024
17 13 T
b − AX (1) = 1260 , 3024.
13 T
D−1 17

b − AX (1) = 630 , 756
1 −1 17 13 T

10 D b − AX (1) = 6300 , 7560
1 −1 37 31 T

X (1) + 10 D b − AX (1) = 6300 , 7560
37 31 T
X (2) = 6300 , 7560
2. Détermination de la matrice itérative G de la méthode et ses valeurs pro-
pres.
9 1

1 −1 1 −1 − 15
G=I− Q A=I− D A= 10
2 9
10 10 − 15 10
Ainsi les valeurs propres de G sont données par
2 ( √ √ )
9 2 9 2 9 2
det(G − λI) = 0 ⇔ −λ − 2 =0 ⇔ λ∈ − , +
10 15 10 15 10 15

3. Pour la converge de la méthode on observe que le rayon spectral

( √ √ ) √
9 2 9 2 9 2
ρ(G) = max |λ| = max − , + = + <1
λ∈sp(A) 10 15 10 15 10 15

La méthode est convergente.

53
Exercice 31 Reprendre l’exercice de l’exemple ci-dessus sur le système linéaire

7x1 − 6x2 = 3
−8x1 + 9x2 = −4

3.3.2 Méthodes itératives standard : Jacobi, Gauss Sei-

del, méthodes de relaxation
Considérons le système Ax = b sous la forme explicite
n
X
aij xj = bi 1≤i≤n
j=1

Nous allons construire explicitement la matrice Q en décomposant la matrice

A de façons simples dans les cas standards.

Méthode de Jacobi
Lorsque pour tout i on a aii 6= 0, en résolvant l’équation numero i pour
l’inconnue xi , on établit la relation de reccurence :
n
(k) (k−1)
X
aii xi =− aij xj + bi 1≤i≤n
j=1
j6=i

ou
 
n
(k) 1  X (k−1)
xi = − aij xj + bi  1≤i≤n

aii j=1
j6=i

On a ainsi formulé point par point la méthode de Jacobi

Dans la suite de cette section on décomposera la matrice A sous la forme
A = D − CL − CU où D est diagonale, CL triangulaire inférieure à diagonale
nulle et CU triangulaire supérieure à diagonale nulle avec

D = diag(A) = (aii )
CL = (−aij )i>j
CU = (−aij )i<j

En formulation matricielle, la méthode de Jacobi s’écrit

Dx(k) = (CL + CU )x(k−1) + b

Ainsi Q = D et la matrice itérative de la méthode de Jacobi est

G = I − D−1 A = D−1 (CL + CU )

L’algorithme de la méthode de Jacobi peut s’écrire :

54
real array (A)n×n ,(b)n ,(c)n , (y)n , (x)n
real sum, diag
integer i, j, k, kmax , n
begin
n ← taille(A)
x ← x(0)
for k = 1 to kmax do
begin
y←x
for i = 1 to n do
begin
sum ← bi
diag ← Aii
if |diag| < δ then
begin
Ecrire ”Element Diagonal trop petit”
return
end
for j = 1 to n do
begin
if j 6= i then
begin
sum ← sum − Aij yj
end
end
xi ← sum/diag
end
Ecrire k, x
if kx − yk < then
begin
Ecrire ”Convergence”
stop
end
end
Ecrire ”Nombre maximal d’iterations atteint”
end

Exemple
On considère le système linéaire
1 1 1
2 x1 + 3 x2 = 63
1 1 1
3 x1 + 4 x2 = 168

T
1. En posant X (k) = xk1 , xk2 la solution à la k ieme itération de la méthode
T
de Jacobi et partant de X (0) = (0 , 0) , calculer X (1) et X (2) .

55
2. Déterminer la matrice itérative G de cette méthode puis calculer ses
valeurs propres.
3. Que peut-t-on dire de la converge de la méthode ?
Pour le système on a :
1 1 1

A= 2 3 et b = 63
1 1 1
3 4 168

Dans la méthode de Jacobi on a :

1 1
I = diag(1 , 1), Q = D = diag , et Q−1 = D−1 = diag(2 , 4)
2 4
1. Calcul de X (1) et X (2) .
Connaissant X (k) on a explicitement :
1 k+1 1 k 1
2 x1 + 3 x2 = 63
1 k 1 k+1 1
x
3 1 + x
4 2 = 168

c’est à dire :
1 k+1 1 1 k
xk+1 1 1 k

2 x1 = 63 − 3 x2 ⇔ 1 = 2 63 − 3 x2
1 k+1
4 x2 = 1 1 k
168 − 3 x1 xk+1
2 = 4 1 1 k
168 − 3 x1

T
Pour X (0) = (0 , 0) on a :
1 2 2
1
x1 = 2 63 − 0 = 63 (1)
1 1 4 ⇔ X = 63
1
x2 = 4 168 − 0 = 168 42

et
1 1 1 1 1

x21

= 2 63 − 3 42 = 126 ⇔ X (2) = 126
1 1 2 −7. 1
x22 = 4 168 − 3 63 = 1512 − 216
2. Calcul de la matrice itérative G de la méthode de Jacobi et ses valeurs
propres.
0 − 23

−1
G=I −D A=
− 43 0
Ainsi
( √ √ )
8 2 2 2 2 2
det(G − λI) = 0 ⇔ (0 − λ) − 2 = 0 ⇔ λ ∈ − ,
3 3 3

3. Le ra yon spectral de la matrice itérative G

( √ √ ) √
2 2 2 2 2 2
ρ(G) = max |λ| = max − , = <1
λ∈sp(A) 3 3 3

Donc la méthode de Jacobi converge.

56
Exercice 32 Reprendre les questions de l’exemple ci-dessus pour chacun des
systèmes linéaires suivants :
1.
7x1 − 6x2 = 3
−8x1 + 9x2 = −4
2.
7x1 − 12x2 = −5
−8x1 + 9x2 = 1
3. 
 5x1 − x2 = 7
−x1 + 3x2 − x3 = −8
− x2 + 2x3 = 4


Méthode de Gauss-Seidel
Intuitement on espère améliorer l’approximation en utilisant les inconnues
x1 , . . . , xj−1 déjà mis à jour lors de l’évaluation de la j ieme inconnue dans le
processus de Jacobi. On établit alors la relation de reccurence :
i−1 n
(k) (k) (k−1)
X X
aii xi =− aij xj − aij xj + bi 1≤i≤n
j=1 j=i−1

ou
 
i−1 n
(k) 1  X (k)
X (k−1)
xi = − aij xj − aij xj + bi  1≤i≤n
aii j=1 j=i+1

On a ainsi la méthode de Gauss-Seidel dont la formulation matricielle est :

Dx(k) = CL x(k) + CU x(k−1) + b

c’est à dire
(D − CL )x(k) = CU x(k−1) + b
et la matrice Q = D − CL est la matrice triangulaire inférieure de A.
L’algorithme de la méthode de Gauss-Seidel s’écrire alors :

real array (A)n×n ,(b)n ,(c)n , (y)n , (x)n

real sum, diag
integer i, j, k, kmax , n
begin
n ← taille(A)
x ← x(0)
for k = 1 to kmax do
begin
y←x

57
for i = 1 to n do
begin
sum ← bi
diag ← Aii
if |diag| < δ then
begin
Ecrire ”Element Diagonal trop petit”
return
end
for j = 1 to i − 1 do
begin
sum ← sum − Aij xj
end
for j = i + 1 to n do
begin
sum ← sum − Aij xj
end
xi ← sum/diag
end
Ecrire k, x
if kx − yk < then
begin
Ecrire ”Convergence”
stop
end
end
Ecrire ”Nombre maximal d’iterations atteint”
end

Exemple
On considère le système linéaire
1 1 1
2 x1 + 3 x2 = 63
1 1 1
3 x1 + 4 x2 = 168

T
1. En posant X (k) = xk1 , xk2 la solution à la k ieme itération de la méthode
T
de Gauss-Seidel et partant de X (0) = (0 , 0) , calculer X (1) et X (2) .
2. Déterminer la matrice itérative G de cette méthode puis calculer ses
valeurs propres.
3. Que peut-t-on dire de la converge de la méthode ?
Pour la méthode de Gauss-Seidel on a :
1 1 1 1

2 3 2 0 63
A= 1 1 Q= 1 1 et b = 1
3 4 3 4 168

58
T
Explicitement, pour X (k) = xk1 , xk2 connu, on a
1 k+1 1 k 1
2 x1 + 3 x2 = 63
1 k+1 1 k+1 1
3 x1 + 4 x2 = 168
qui donne le système triangulaire inférieur
1 k+1 1
k+1
− 13 xk2 1 1 k

2 x1 = 63
⇔
x1 = 2 63 − 3 x2
1 k+1 1 k+1 1
xk+1 1 1 k+1

x
3 1 + x
4 2 = 168 2 = 4 168 − 3 x1
T
En partant de X (0) = (0 , 0) on a
1 2
− 13 x02 2
1 1
x1 = 2 63 x1 =
1 ⇔ 63
1 1 2 ⇔ X (1) = 63
1
− 13 x11

x12 = 4 168 x12 = 4 168 − 3 63 − 54

1 1 1 1 1 1 25

x21 x21

= 2 63 − 3 x2 ⇔
= 2 63 + 3 54 ⇔ X (2) = 567
1 1 2 1 1 25 119
x22 = 4 168 − 3 x1 x22 = 4 168 − 3 567
− 3402

Méthode de Surrelaxation : SOR

Ici on introduit un facteur de relaxation ω sur deux itérations successives de
Gauss-Seidel. Ainsi on pourra écrire
(k) (k) (k−1)
xSOR = ωx + (1 − ω)x
Gauss-Seidel Gauss-Seidel
On a donc
 
i−1 n
(k) ω X (k)
X (k−1) (k−1)
xi = − aij xj − aij xj + bi  + (1 − ω)xi 1≤i≤n
aii j=1 j=i+1

La formulation matricielle s’écrit :

(D + ωCL )x(k) = [ωCU + (1 − ω)D]x(k−1) + ωb

Exercice 33 Ecrire l’algorithme de la méthode de sur-relexation (SOR)

Exercice 34 On considère le système linéaire

1 1 1
2 x1 + 3 x2 = 63
1 1 1
3 x 1 + 4 x2 = 168
T
1. En posant X (k) = xk1 , xk2 la solution à la k ieme itération de la méthode
11
de sur-relaxation avec comme coefficient de relaxation ω = et partant
10
T
de X (0) = (0 , 0) , calculer X (1) et X (2) .
2. Déterminer la matrice itérative G de cette méthode puis calculer ses valeurs
propres.
3. Que peut-t-on dire de la converge de la méthode ?

59
3.3.3 Convergence des méthodes
Théorème 3.3.2 Si la matrice A ∈ R I n×n est à diagonale strictement dom-
inante alors les méthodes de Jacobi et de Gauss-Seidel appliquées au système
Ax = b sont convergentes pour toute solution initiale x(0)

Preuve
La démontration sera faite en exercice.
Remarque :
Le théorème 3.3.2 donne une condition suffisante de convergence qui n’est absol-
ument pas nécessaire. On a d’autres critères de convergence ; par exemple dans
le cas des matrices hermitiennes.

Exercice 35 Soit Dx(k) = (D − A)x(k−1) + b l’itération de Jacobi appliquée

à l’équation Ax = b où A est à diagonale strictement dominante. Soit G =
I − D−1 A
Montrer kGk∞ < 1. En déduire la convergence de la méthode de Jacobi.

Exercice 36 Montrer que la méthode de Jacobi converge pour le système Ax =

I 2×2 symétrique définie positive.
b pour A ∈ R

Exercice 37 Comparer les rayons spectraux ρ(I − Q−1 −1

J A) et ρ(I − QG A) re-
spectifs des méthodes de Jacobi et de Gauss-Seidel pour la matrice
 
4 −1 −1
A =  −1 4 −1 
−1 −1 4

3.3.4 Effet de l’arithmétique finie

On va simplement mettre en évidence l’influence de l’arithmétique à travers
un exemple illustratif.
Exemple
Considérons à cet effet le système :

 3x1 − x2 + x3 = 1
3x1 + 6x2 + 2x3 = 0
3x1 + 3x2 + 7x3 = 4


En utilisant une arithmétique décimale tronquée à 4 chiffres et l’une des

T
méthodes itératives suivantes, calculons X (1) et X (2) pour X (0) = (0 , 0 , 0)
1. Méthode de Jacobi
2. Méthode de Gauss-Seidel
3. Méthode de sur-relaxation avec un coefficient de relaxation ω = 1.25
T
En posant X (k) = xk1 , xk2 , xk3

60
1. Pour la méthode de Jacobi, connaissant X (k) , on a explicitement :
 k+1
= 31 1 + xk2 − xk3

 x1
xk+1 = 16 −3xk1 − 2xk3
 2k+1
x3 = 17 4 − 3xk1 − 3xk2
On obtient en arithmétique finie
 k+1
 x1 = 0.3333 1 + xk2 − xk3
xk+1 = 0.1666 −3xk1 − 2xk3
 2k+1
x3 = 0.1428 4 − 3xk1 − 3xk2
T
Pour X (0) = (0 , 0 , 0) , on a :
 1
 x1 = 0.3333(1 + 0 − 0) = 0.3333
x1 = 0.1666(−3 ∗ 0 − 2 ∗ 0) = 0
 21
x3 = 0.1428(4 − 3 ∗ 0 − 3 ∗ 0) = 0.5712
m 
0.3333
X (1) =  0 
0.5712
et  2
 x1 = 0.3333(1 + 0 − 0.5712) = 0.1429
x2 = 0.1666(−3 ∗ 0.3333 − 2 ∗ 0.5712) = −0.3566
 22
x3 = 0.1428(4 − 3 ∗ 0.3333 − 3 ∗ 0) = 0.4285
m 
0.1429
X (2) =  −0.3566 
0.4285
2. La méthode de Gauss-Seidel, connaissant X (k) , s’explicite sous la forme
 k+1
= 13 1 + xk2 − xk3

 x1
xk+1 = 16 −3xk+1 − 2xk3
 2k+1 1
= 7 4 − 3x1 − 3xk+1
1 k+1

x3 2

On obtient en arithmétique finie

 k+1
 x1 = 0.3333 1 + xk2 − xk3
k+1
x = 0.1666 −3xk+1 − 2xk3
 2k+1 1
0.1428 4 − 3x1 − 3xk+1
k+1

x3 = 2
T
Pour X (0) = (0 , 0 , 0) , on a :
 1
 x1 = 0.3333(1 + 0 − 0) = 0.3333
x1 = 0.1666(−3 ∗ 0.3333 − 2 ∗ 0) = −0.1665
 12
x3 = 0.1428(4 − 3 ∗ 0.3333 − 3 ∗ (−0.1665)) = 0.4998
m 
0.3333
X (1) =  −0.1665 
0.4998

61
et
 2
 x1 = 0.3333(1 − 0.1665 − 0.4998) = 0.1113
x2 = 0.1666(−3 ∗ 0.1113 − 2 ∗ 0.4998) = −0.2220
 22
x3 = 0.1428(4 − 3 ∗ 0.1113 − 3 ∗ (−0.2220)) = 0.6186
m 
0.1113
X (2) =  −0.2220 
0.6186

3. La méthode de sur-relaxation avec un coefficient de relaxation ω = 1.25,

connaissant X (k) , peut s’écrire :
 k+1
= 1.25

 x1 3 1 + xk2 − xk3 − 0.25xk1
xk+1
2 = 1.25 6 −3xk+1
1 − 2xk3 − 0.25x k
2
 k+1 1.25 k+1 k+1

x3 = 7 4 − 3x1 − 3x2 − 0.25xk3

On obtient en arithmétique finie

 k+1
 x1 = 0.4166 1 + xk2 − xk3 − 0.25xk1
xk+1 = 0.2083 −3xk+1 − 2xk3 − 0.25x k
2
 2k+1 1
k+1 k+1

x3 = 0.1785 4 − 3x1 − 3x2 − 0.25xk3

T
Pour X (0) = (0 , 0 , 0) , on a :
 1
 x1 = 0.4166(1 + 0 − 0) − 0.25 ∗ 0 = 0.4166
x1 = 0.2083(−3 ∗ 0.4166 − 2 ∗ 0) − 0.25 ∗ 0 = −0.2601
 21
x3 = 0.1785(4 − 3 ∗ 0.4166 − 3 ∗ (−0.2601)) − 0.25 ∗ 0 = 0.6302
m 
0.4166
X (1) =  −0.2601 
0.6302

et
 2
 x1 = 0.4166(1 + (−0.2601) − 0.6302) − 0.25 ∗ 0.4166 = −0.0583
x2 = 0.2083(−3 ∗ (−0.0583) − 2 ∗ 0.6302) − 0.25 ∗ (−0.2601) = −0.1612
 22
x3 = 0.1785(4 − 3 ∗ (−0.0583) − 3 ∗ (−0.1612)) − 0.25 ∗ 0.6302 = 0.6739
m 
−0.0583
X (2) =  −0.1612 
0.6739

Exercice 38 Reprendre les questions de l’exemple ci-dessus avec une arithmétique

décimale arrondie à 4 chiffres.

62
3.4 Accélération de convergence
La relaxation est une technique beaucoup plus générale notamment utilisée
pour accélérer la convergence d’une méthode itérative. Dans sa version simple,
elle se formule par :

Qy (k) = (Q − A)x(k−1) + b

k = 1, 2, . . .
x(k) = ωy (k) + (1 − ω)x(k−1)

où ω est le coefficient de relaxation.

Remarque :
La méthode SOR est une relaxation de la méthode de Gauss-Seidel.
Nous pouvons définir une méthode de Jacobi relaxée par :

x(k) = x(k−1) + ω b − Ax(k−1) k = 1, 2, . . .

On a aussi l’habitude de mettre en oeuvre des techniques plus générale

de type polynômial en particilier la méthode de Tchebychev. La méthode de
Tchebychev comme celle de gradient conjugué qui sont d’une grande éfficacité
sont au delà des ambitions de ce cours d’introduction.

63
Chapitre 4

Résolution numérique des

équations différentielles
ordinaires

L’analyse qualitative des solutions d’une équation ou d’un système différentiel

ordinaire n’est pas à l’ordre du jour dans ce cours (amateur, voir [8]). Nous al-
lons juste nous restreindre à la résolution numérique des problèmes de valeurs
initiales associés à des équations différentielles ordinaires. Un accent particulier
sera mis sur le traitement en arithmétique finie avec les algorithmes proposés.
Loin de recourir à la mémorisation des algorithmes, l’objectif ici est surtout
de les utiliser concrètement et de pouvoir comparer leur efficacité à celle des
méthodes de Taylor.

4.1 Equation différentielle et classification

Définition 4.1.1 Toute équation impliquant des dérivées d’une ou plusieurs
fonctions d’une ou plusieurs variables est dite équation différentielle.

Remarque :
Dans cette définition, on exclut les identités comme

d2 x d dv du
(e ) = ex ou (uv) = u +v
dx2 dx dx dx

64
Exemple
2
d2 y dy
dx2 + xy dx =0

d4 x 2

dt4 + 5t ddt2x + 3x = cos t

∂v ∂v
∂t + ∂x =0

∂2u ∂2u ∂2u

∂x2 + ∂y 2 + ∂z 2 =0
Remarque :
Lorsque la fonction inconnue est à une variable, on dit qu’on a une équation
différentielle ordinaire (EDO). Dans le cas contraire on a une équation aux
dérivées partielles (EDP).
Remarque :
L’ordre le plus élevé de dérivation intervenant dans une équation différentielle
est appelé l’ordre de l’équation différentielle.

4.1.1 Equation différentielle et solution

Plus formellement
Définition 4.1.2 Pour n ∈ N I ∗ , une équation différentielle ordinaire d’ordre n
est une équation de la forme

dy d2 y dn−1 y dn y

F t, y, , 2 , . . . , n−1 , n = 0 (4.1)
dt dt dt dt
où
I × E n+1 → E
F :R
I p , p ≥ 1).
avec E est un espace vectoriel normé ( par exemple E = R
Remarque :
I p avec p > 1 on parlera de système d’équations différentielles.
Lorsque E = R

Définition 4.1.3 Une fonction f : I ⊂ R I → E est dite solution explicite de

l’équation 4.1 sur I si f est n fois dérivable sur I et

x 7→ F x, f (x), f 0 (x), . . . f (n−1) (x), f (n) (x),

est définie en tout point de I avec

F x, f (x), f 0 (x), . . . f (n−1) (x), f (n) (x), = 0 ∀x ∈ I

Définition 4.1.4 Une relation g(x, y) = 0 où g est une fonction R I × E → E,

est dite solution implicite de l’équation 4.1 sur I si cette relation définie une
solution explicite de 4.1 sur I

65
Exemple
d2 y
f (x) = 2 cos x − 3 sin x est solution de dx 2 + y = 0 sur RI car f est deux fois
00
I avec f 0 (x) = −2 sin x − 3 cos x et f (x) = −2 cos x + 3 sin x d’où
dérivable sur R
00
f (x) + f (x) = (−2 cos x + 3 sin x) + (2 cos x − 3 sin x) = 0 ∀x ∈ R
I

Exemple
dy
La relation x2 +y 2 −25 = 0 est une solution implicite x+y dx = 0 sur I =]−5, 5[
2 2
En effet
√ sur I la relation x + y − 25 = 0 définit la fonction f donnée par
f (x) = 25 − x2 qui est dérivable sur I avec
x
f 0 (x) = − √ et x + f (x)f 0 (x) = x − x = 0 ∀x ∈ I
25 − x2
dy
Ainsi f est solution de x + y dx = 0 sur I d’où x2 + y 2 − 25 = 0 est une solution
dy
implicite de x + y dx = 0.
Remarque :
dy
De la relation x2 + y 2 + 25 = 0 on a par dérivation 2x + 2y dx = 0 c’est à dire
dy
x + y dx = 0 mais pour tout x ∈ R I

x2 + y 2 + 25 = 0 ⇒ 0 ≤ x2 + y 2 = −25 < 0
dy
ce qui est absurde. Ainsi x2 + y 2 + 25 = 0 n’est pas une solution de x + y dx =0
Nous laissons quelques exemples à l’appréciation du lecteur.

Equation a pour Solution

0 t
x (t) − x(t) = e I 7→ x(t) = tet + cet
t∈R (4.2)
00
y + 9y = 0 t∈R
I 7→ y(t) = c1 sin 3t + c2 cos 3t (4.3)
1
u0 + =0 u2 (t) + t = c sur ] − ∞, c[ (4.4)
2u
Remarque :
Dans les exemples énumérés ci-dessus nous pouvons constater que la solution
d’une équation différentielle ordinaire est connue à une constante réelle c près.
Nous avons donc une indétermination sur la solution. Pour connaı̂tre précisément
la solution, il est nécessaire d’adjoindre des conditions exprimant certaines con-
traintes sur la solution recherchée.

4.1.2 Problème de Cauchy

Dans notre cours nous nous interesserons essentiellement aux problèmes dits
de valeur initiale ou problème de Cauchy.
Exemple
Trouvons la solution f de l’équation différentielle
dy
= −2x
dx

66
tel que f (1) = 4.
Nous pouvons vérifier que y(x) = −x2 + c est une famille à un paramètre de
solution de cette EDO. Déterminons cette constante c pour que y(1) = 4. Nous
obtenons l’equation −1 + c = 4 d’où c = 5. La solution f sera alors donnée par
f (x) = −x2 + 5
Le problème que nous venons de résoudre se note simplement
dy
dx = −2x
y(1) = 4
Exemple
De façon similaire on peut trouver la solution du problème du second ordre
suivant  2
d y
 dx 2 + y = 0

y(1) = 3
 0
y (1) = −4
qui est aussi problème de valeur initiale possedant une solution unique.
En effet la solution générale de cette EDO est connue et y(t) = a sin(t) + b cos(t).
Donc les contraintes donnent le système :

y(1) = 3 a sin(1) + b cos(1) = 3
⇔
y 0 (1) = −4 a cos(1) − b sin(1) = −4

qui est de Cramer car son déterminant D = − sin2 (1) − cos2 (1) = −1. On a donc
I 2.
une solution unique (a , b) ∈ R
Remarque :
Nous pouvons associer un autre type de contraintes à l’équation
d2 y
+y =0
dx2
que nous notons
d2 y

 dx2+y =0
y(0) = 1
y( π2 ) = 5


Il s’agit ici d’un problème de valeur aux limites qui a aussi une solution unique
f (x) = cos(x) + 5 sin(x).
En effet les contraintes donnent

y(0) = 1 a sin(0) + b cos(0) = 1 b=1
⇔ ⇔
y( π2 ) = 5 a sin( π2 ) + b cos( π2 ) = 5 a=5
Nous devons cependant noter qu’un problème de valeurs aux limites n’est
pas à prendre à la légere car par exemple, le problème
 2
d y
 dx 2 + y = 0

y(0) = 1

y(π) = 5

67
n’a pas de solution car les contraintes conduisent au système :

y(0) = 1 a sin(0) + b cos(0) = 1 b=1
⇔ ⇔
y(π) = 5 a sin(π) + b cos(π) = 5 −b = 5

ce qui est absurde.

4.1.3 Problème de Cauchy d’ordre 1

Dans toute la suite nous traiterons de problème d’ordre 1 qu’on définira
formellement comme suit :

I p, p ≥
Définition 4.1.5 E est un espace vectoriel normé (par exemple E = R
1). Considérons l’équation différentielle ordinaire d’ordre 1

dy
= f (x, y) (4.5)
dx
où f est une fonction continue en (x, y) ∈ D, D un domaine de R I × E. Soit
(x0 , y0 ) ∈ D un point donné.
Le problème au valeur intiale ou problème de Cauchy associé à l’équation
(4.5) est dit problème de Cauchy d’ordre 1 et consiste en la recherche d’une
solution φ de l’équation différentielle (4.5) définie sur un intervalle I de R
I
contenant x0 à valeurs dans E et vérifiant

φ(x0 ) = y0

On note très simplement

dy

dx= f (x, y)
y(x0 ) = y0

Exemple
On peut lister quelques problèmes de Cauchy que nous manipulerons par la
suite.
1. Un problème linéaire
du


 = u + 1 pour 0 ≤ t ≤ 1
dt


u(0) = 0

2. Un problème non linéaire

du u u 2


 = 1+ + pour 1 ≤ t ≤ 2
dt t t


u(1) = 0

68
3. Un système linéaire

d u v π
= pour 0 ≤ t ≤


 dt v −u 2



u 0


(0) =


v 1


D’une façon générale un problème de Cauchy d’ordre n > 1 peut se ramener

à un problème de Cauchy d’ordre 1.
Exemple
Le problème d’ordre 2 :
2

 d u = −u pour 0 ≤ t ≤ π

dt2



 2


u(0) = 0




 du (0) = 1



dt
est équivalent au problème d’ordre 1 :

d u v π
= pour 0 ≤ t ≤


 dt v −u 2



u 0


(0) =


v 1


du
moyennant la variable intermédiaire v = .
dt
Exercice 39 Pour chacun des problèmes de Cauchy suivants, déterminer le
problème de Cauchy d’ordre 1 associé
1.  2
d u du
+ + u2 + 2t = 0 pour 0 ≤ t ≤ 1


dt2



 dt


u(0) = 0





 du
(0) = 0.21


dt
2.
3 2

3d u 2d u du
t − t + 3t − 4u = 5t3 ln t + 9t3 pour 1 ≤ t ≤ 2




 dt 3 dt 2 dt



u(1) = 0




du d2 u


(1) = 1 (1) = 3


dt dt2

69
3.
d2 u

du dv
= −2 +u+3 + v + ln t pour 0 ≤ t ≤ 1


dt2 dt dt







d2 v du dv


= 5 + 2u + t − 3v − sin t


dt2

 dt dt
du


u(0) = 1 (0) = 2





 dt





 v(0) = 3 dv
(0) = 4

dt

4.2 Résolution d’un problème de Cauchy

Soit E est un espace vectoriel normé (par exemple E = R I p , p ≥ 1). Soit Ω
un domaine borné de E et D = [a , b] × Ω̄. Soit f est une fonction continue sur
D. Soit (x0 , y0 ) ∈ D un point donné.
Il est hors du champs de ce cours l’étude de l’existence et l’unicité du
problème de Cauchy : dy
dx = f (x, y) (4.6)
y(x0 ) = y0
Cependant il n’est pas superflu avant de se lancer dans de fastidueux calcul
d’être sûr que ce ne sera pas peine perdue.

4.2.1 Résultats d’existence et unicité

Le cadre théorique de justification de l’existence et de l’unicité d’une solution
locale d’un problème de Cauchy d’ordre 1 s’appuie sur le Théorème de Cauchy
Picard qui peut s’énoncer :

Proposition 4.2.1 (Cauchy-Picard) Si la fonction f est continue et bornée

sur D, si de plus f vérifie la condition de lipschitz relativement à la seconde
variable c’est à dire il existe un réel k > 0 tel que

kf (x , y1 ) − f (x , y2 )k ≤ kky1 − y2 k pour tout (x , y1 ) ∈ D et tout (x , y2 ) ∈ D

alors le problème de Cauchy (4.6), admet une et une seule solution φ définie un
intervalle fermé I ⊂ [a , b] contenant x0

Preuve
Cette démonstration, que l’on peut trouver dans un cours de calcul différentiel
en Licence, sort du cadre de la présente note de cours. (Pour une approche
élémentaire on peut se reférer au livre de Shepley L. Ross [7]).
Remarque :
La proposition a des déclinaisons plus faciles d’emploi comme la suivante, énoncée
dans le cas E = R I et dont nous ne donnerons pas non plus la démonstration.

70
Proposition 4.2.2 Si f est continue sur D = [a , b]×IR et si sa dérivée partielle
fy par rapport à sa seconde variable est continue sur D alors pour tout y0 ∈ R I
le problème de Cauchy
dy
dt = f (t , y) pour a ≤ t ≤ b
y(a) = y0
admet une solution unique φ définie sur [a , b].

4.2.2 Résolution numérique des EDO et l’intégration

Les techniques de résolution numérique des EDO sont pour la plupart liées
aux méthodes d’approximation des intégrales définies. Ce lien apparaı̂t dans la
proposition suivante :
Proposition 4.2.3 Soit f une fonction réelle définie et continue sur un do-
maine D ∈ RI 2 . Soit φ une fonction réelle définie et continue sur un intervalle
α ≤ x ≤ β tel que (x, φ(x)) ∈ D pour tout x ∈ [α, β]. Soit x0 un réel tel que
α < x0 < β.
Alors φ est une solution du problème de Cauchy
dy
dx = f (x, y) x ∈ [α, β] (4.7)
y(x0 ) = y0
si et seulement si
Z x
φ(x) = y0 + f (t, φ(t))dt ∀x ∈ [α, β] (4.8)
x0

Preuve
dy
Si φ est solution de l’équation différentielle dx = f (x, y) sur [α, β] alors pour
tout x ∈ [α, β] on a
dφ(x)
= f (x, φ(x))
dx
par intégration on a
Z x
φ(x) − φ(x0 ) = f (t, φ(t))dt
x0

Comme φ(x0 ) = y0 on obtient (4.8)

Réciproquement si φ satisfait l’équation intégrale (4.8) alors par dérivation
on obtient
dφ(x)
= f (x, φ(x)) ∀x ∈ [α, β]
dx
Donc φ satisfait le problème de Cauchy ( 4.7).
Remarque :
Considérons x une solution dy dt = f (t, y) sur ]a, b[. Soit t ∈]a, b[ et h ∈ R
I tel que
t + h ∈]a, b[. Alors on
Z t+h
x(t + h) − x(t) = f (s, x(s))ds
t

71
L’approximation de l’intégrale détermine la méthode.
Exemple

1. La méthode de rectangle
Z t+h
f (s, x(s))ds = hf (t, x(t)) + E
t

conduit à la méthode d’Euler

x(t + h) = x(t) + hf (t, x(t))

2. La méthode de rectangle
Z t+h
f (s, x(s))ds = hf (t + h, x(t + h)) + E
t

conduit à la méthode d’Euler implicite

x(t + h) = x(t) + hf (t + h, x(t + h))

dont le traitement final nécessite la résolution d’une équation où x(t + h)

est l’inconnue. Cette tache peut se réveler hardue
3. La méthode de trapèze
Z t+h
h
f (s, x(s))ds = [f (t, x(t)) + f (t + h, x(t + h))] + E
t 2
donne également une méthode implicite
h
x(t + h) = x(t) + [f (t, x(t)) + f (t + h, x(t + h))]
2
On peut contourner la difficulté en prenant une valeur prédicte de x(t +
h) par la méthode d’Euler dans f (t + h, x(t + h)). On obtient alors une
méthode de Runge Kutta d’ordre 2
h
x(t + h) = x(t) + [f (t, x(t)) + f (t + h, x(t) + hf (t, x(t)))]
2

4.3 Méthodes de Taylor

L’approximation de la solution du problème de valeur initiale semble être
une application naturelle du Théorème de Taylor. Ici on développe la fonction
inconnue en polynôme de Taylor. La forme la plus élémentaire conduit à la
méthode d’Euler qui est d’une importance tant théorique que pratique.
Dans cette cette section comme dans les suivantes, la fonction f est réelle
sauf indication contraire ; les résultats énoncés sont trivialement extensibles au
cas vectoriel.

72
4.3.1 Méthode d’Euler
Le but de la méthode d’Euler est de trouver une approximation du problème
de Cauchy dy
dt = f (t, y) t ∈ [a, b]
y(a) = α

Théorème 4.3.1 Soit ti = a + ih pour i = 0 . . . N un réseau régulier de pas

h = b−a
N . Soit yi une approximation par la méthode d’Euler de y(ti ), solution
du problème de Cauchy. La méthode d’Euler est donnée par :

y0 = α
yi+1 = yi + hf (ti , yi ) ∀i = 0, . . . , N − 1

h2 (2)
Si la solution est deux fois dérivable sur [a, b] alors l’erreur locale est 2 y (ξ)
avec ξ ∈]ti , ti+1 [

Preuve
Il découle immédiatement de la formule de Taylor - Lagrange de la page 7 pour
l’ordre n=1.
L’algorithme de la méthode d’Euler est le suivant :

real array (yi )n , (xi )n

function externe f
real h, t0 , tf in , y0
integer k
begin
h ← (tf in − t0 )/(n − 1)
t ← t0
u ← y0
x1 ← t
y1 ← u
Ecrire t, u
for k = 2 to n do
begin
u ← u + hf (t , u)
t←t+h
xk ← t
yk ← u
Ecrire t, u
end
end

real function f (t , y)
real t, y
integer k
begin

73
A definir
end

Exemple
Considérons le problème de Cauchy
0
y = y − t2 + 1 0≤t≤2
y(0) = 0.5
10
Déterminer une approximation (yi )i=0 de la solution sur un réseau régulier par
la méthode d’Euler.
On observe que N = 10 donc h = N2 = 0.2 et ti = 0.2i pour i = 1, . . . , 10.
On a alors l’expression

yi+1 = yi + h(yi − t2i + 1) ∀i = 0, . . . , N − 1

d’où la formule de récurrence

y0 = α
yi+1 = 1.2yi − 0.008i2 + 0.2 ∀i = 0, . . . , 9

C’est une équation aux différences que nous savons pas toujours résoudre ex-
plicitement même dans un cas linéaire comme ici.
Remarque :
En pratique il n’est nullement nécessaire d’établir une relation de réccurence
avant de faire le calcul numérique de la solution en arithmétique finie qui est un
objectif de ce cours.
Exemple
En utilisant une arithmétique décimale arrondie à 4 chiffres, calculer sur une
grille uniforme de pas h = 0.2 par la méthode d’Euler la solution du problème
de Cauchy :
dy y y 2


 = 1+ + pour 1 ≤ t ≤ 2
dt t t


y(1) = 0
y y 2
Soit f (t , y) = 1 + + .
t t
En notant ŷ(t) l’approximation de y(t), une itération de la méthode d’Euler
se traduit par : pour h = 0.2, t et ŷ(t) connus,

ŷ(t + h) = ŷ(t) + hf (t , ŷ(t))

Ainsi on a la séquence :

74
1. Pour t = 1, ŷ(t) = 0 on a
ŷ(t)
t =0
f (t , ŷ(t)) = 1
hf (t , ŷ(t)) = 0.2
ŷ(t) + f (t , ŷ(t)) = 0 + 0.2 = 0.2
t + h = 1.2
ŷ(t + h) = 0.2
2. Pour t = 1.2, ŷ(t) = 0.2 on a
ŷ(t)
t = 0.1667
f (t , ŷ(t)) = 1 + (0.1667) + (0.1667)2 = 1.195
hf (t , ŷ(t)) = 0.2 ∗ 1.195 = 0.239
ŷ(t) + f (t , ŷ(t)) = 0.2 + 0.239 = 0.439
t + h = 1.4
ŷ(t + h) = 0.439
3. Pour t = 1.4, ŷ(t) = 0.439 on a
ŷ(t)
t = 0.3136
f (t , ŷ(t)) = 1 + (0.3136) + (0.3136)2 = 1.412
hf (t , ŷ(t)) = 0.2 ∗ 1.412 = 0.2824
ŷ(t) + f (t , ŷ(t)) = 0.439 + 0.2824 = 0.7214
t + h = 1.6
ŷ(t + h) = 0.7214
4. Pour t = 1.6, ŷ(t) = 0.7214 on a
ŷ(t)
t = 0.4509
f (t , ŷ(t)) = 1 + (0.4509) + (0.4509)2 = 1.654
hf (t , ŷ(t)) = 0.2 ∗ 1.654 = 0.3308
ŷ(t) + f (t , ŷ(t)) = 0.7214 + 0.3308 = 1.052
t + h = 1.8
ŷ(t + h) = 1.052
5. Pour t = 1.8, ŷ(t) = 1.052 on a
ŷ(t)
t = 0.5844
f (t , ŷ(t)) = 1 + (0.5844) + (0.5844)2 = 1.926
hf (t , ŷ(t)) = 0.2 ∗ 1.926 = 0.3852
ŷ(t) + f (t , ŷ(t)) = 1.052 + 0.385 = 1.437
t+h=2
ŷ(t + h) = 1.437
En résumé on a
t 1 1.2 1.4 1.6 1.8 2
ŷ(t) 0 0.2 0.439 0.7214 1.052 1.437

75
4.3.2 Majoration d’erreur de la méthode d’Euler
Cette majoration importante au plan théorique sera donnée à titre indicatif
et ne fera pas l’objet d’une démonstration dans ce cours.

Théorème 4.3.2 Soit y(t) l’unique solution du problème de Cauchy

dy
dt = f (t, y) t ∈ [a, b]
y(a) = α
N
et (yi )i=0 l’approximation générée par la méthode d’Euler sur un réseau régulier,
pour un N positif donné.
Si la fonction f est continue pour tout t ∈ [a, b] et pour tout y ∈ R I , et s’il
existe deux réels positifs A et B tels que
∂f
k (t, y(t))k ≤ A et ky”(t)k ≤ B
∂t
alors pour tout i on a
hB A(ti −a)
|y(ti ) − yi | ≤ e −1
2A
Preuve
(Les fans d’analyse amusez vous ! On peut simplement consulter [4] )
Remarque :
Ce resultat théorique statue sur la convergence de la méthode d’Euler.

4.3.3 Méthode de Taylor d’ordre n

Théorème 4.3.3 Soit ti = t0 + ih pour i = 0 . . . N un réseau régulier de pas
t −t
h = f inN 0 . Soit yi une approximation par la méthode de Taylor d’ordre n
de y(ti ), solution du problème de Cauchy. La méthode de Taylor d’ordre n est
donnée par :

y0 = α
yi+1 = yi + hT (n) (ti , yi ) ∀i = 0, . . . , N − 1
avec
h 0 h2 hn−1 (n−1)
T (n) (ti , yi ) = f (ti , yi ) + f (ti , yi ) + f (2) (ti , yi ) + . . . + f (ti , yi )
2! 3! n!
hn+1 (n+1)
Si la solution est n+1 fois dérivable sur [a, b] alors l’erreur locale est (n+1)! y (ξ)
avec ξ ∈]ti , ti+1 [.

Preuve
Application immédiate de Taylor-Lagrange (voir page 7)

76
Exemple
Considérons le problème de Cauchy
0
y = y − t2 + 1 0≤t≤2
y(0) = 0.5
10
Déterminer une approximation (yi )i=0 de la solution sur un réseau régulier :
1) par la méthode de Taylor d’ordre 2.
2) par la méthode de Taylor d’ordre 4.

Dans cet exemple, on a f (t, y(t)) = y(t) − t2 + 1. Ainsi

0
f (t, y(t)) = d
dt (y − t2 + 1) = y 0 − 2t = y(t) − t2 − 2t + 1
00
f (t, y(t)) = d
dt (y − t2 − 2t + 1) = y 0 − 2t − 2 = y(t) − t2 − 2t − 1
000
f (t, y(t)) = d
dt (y − t2 − 2t − 1) = y 0 − 2t − 2 = y(t) − t2 − 2t − 1

Donc pour N = 10 on a h = 0.2 et ti = 0.2i, puis

0
T (2) = f (ti , yi) + h2 f (ti , yi ) = yi − t2i + 1 + h2 (yi − t2i − 2t + 1)
= 1 + h2 (yi − t2i + 1) − hti
= 1.1yi − 0.044i2 − 0.04i + 1.1

et
0 2 00 3 000
T (4) = f (ti , yi ) + h2 f (ti , yi ) + h6 f (ti , yi ) + h24 f (ti , yi )
= yi − t2i + 1 + h2 (yi − t2i − 2t + 1)
2 3
+ h6 (yi − t2i − 2t − 1) + h24 (yi − t2i − 2t − 1)
2 3 2
h2 h3
= 1 + h2 + h6 + h24 (yi − t2i ) − 1 + h3 + h12 hti + 1 + h
2 − 6 − 24
= 1.107yi − 0.04428i2 − 0.0428i + 1.093

On déduit ainsi les formules de réccurence :

Pour Taylor d’ordre 2

y0 = 0.5
yi+1 = 1.22yi − 0.0088i2 − 0.08i + 0.22 ∀i = 0 . . . 9

Pour Taylor d’ordre 4

y0 = 0.5
yi+1 = 1.221yi − 0.008856i2 − 0.0856i + 0.218 ∀i = 0 . . . 9

Comme on l’a déjà signalé plus haut, cette démarche est fastidueuse et ne
donne pas les résultats numériques.
L’algorithme de la méthode de Taylor d’ordre 3 par exemple, peut s’écrire :

77
real array (yi )n , (xi )n
function externe f , f1 , f2
real h, t0 , tf in , y0
integer k
begin
h ← (tf in − t0 )/(n − 1)
t ← t0
u ← y0
x1 ← t
y1 ← u
Ecrire t, u
for k = 2 to n do
begin
u1 ← f (t , u)
u2 ← f1 (t , u , u1 )
u3 ← f2 (t , u , u1 , u2 )
u ← u + h u1 + h u22 + u3

6 h
t←t+h
xk ← t
yk ← u
Ecrire t, u
end
end

real function f (t , y)
real t, y
integer k
begin
A definir
end

real function f1 (t , y , y1 )
real t, y, y1
integer k
begin
A definir
end

real function f2 (t , y , y1 , y2 )
real t, y, y1 , y2
integer k
begin
A definir
end

78
Donnons une illustration numérique des méthodes de Taylor en reprenant
l’exemple traité par la méthode d’Euler.
Exemple
On considère le problème de Cauchy
dy y y 2


 = 1+ + pour 1 ≤ t ≤ 2
dt t t


y(1) = 0

Calculer sur une grille uniforme de pas h = 0.2 la solution du problème de

Cauchy
1. En utilisant une arithmétique décimale arrondie à 4 chiffres et la méthode
de Taylor à l’ordre 2
2. En utilisant une arithmétique décimale tronquée à 4 chiffres et la méthode
de Taylor à l’ordre 3
y y 2
Soit f (t , y) = 1 + + .
t t
d d2
On va avoir besoin des fonctions f1 = f et f2 = 2 f Ainsi
dt dt
0
y y y
f1 (t , y , y 0 ) = − 2 1+2
t t t
1 0 y y
= y − 1+2
t t t 2
00 0
0
0 00 y y y y y y
f2 (t , y , y , y ) = −2 2 +2 3 1+2 +2 − 2
t t t t t t
1 00 1 0 y
h i y 1 0 y 2
= y −2 y − 1+2 +2 2 y −
t t t t t t

En notant ŷ(t) l’approximation de y(t),

1. une itération de la méthode de Taylor d’ordre 2 se traduit par :
pour h = 0.2, t et ŷ(t) connus,

ŷ(t + h) = ŷ(t) + hy 0 (t) + 21 h2 y 00 (t)

avec
y 0 (t) = f (t , ŷ(t))
y 00 (t) = f1 (t , ŷ(t) , y 0 (t))

Ainsi en arithmétique arrondie à 4 chiffres on a la séquence :

0
(a) Pour t = 1, ŷ(t) = 1 = 0 on a

79
ŷ(t)
y = t =0
y1 = 1 + y + y2 = 1
1
y2 = t (y1 − y)(1 + 2y)
= (1/1) ∗ (1 − 0) ∗ (1 + 2 ∗ 0) = 1
ŷ(t + h) = ŷ(t) + h(y1 + 0.5hy2 )
= 0 + 0.2 ∗ (1 + 0.1 ∗ 1.) = 0.22
t+h = 1.2
ŷ(t + h) = 0.22
(b) Pour t = 1.2, ŷ(t) = 0.22 on a
ŷ(t)
y = t = 0.22
1.2 = 0.1833
y1 = 1 + y(1 + y) = 1 + 0.1833 ∗ (1 + 0.1833) = 1.217
1
y2 = t (y1 − y)(1 + 2y)
= (1/1.2) ∗ (1.217 − 0.1833) ∗ (1 + 2 ∗ 0.1833) = 1.178
ŷ(t + h) = ŷ(t) + h(y1 + 0.5hy2 )
= 0.22 + 0.2 ∗ (1.178 + 0.1 ∗ 1.178) = 0.4792
t+h = 1.4
ŷ(t + h) = 0.4792
(c) Pour t = 1.4, ŷ(t) = 0.4792 on a
ŷ(t)
y = t = 0.4792
1.4 = 0.3423
y1 = 1 + y(1 + y) = 1 + 0.3423 ∗ (1 + 0.3423) = 1.459
1
y2 = t (y1 − y)(1 + 2y)
= (1/1.4) ∗ (1.459 − 0.3423) ∗ (1 + 2 ∗ 0.3423) = 1.344
ŷ(t + h) = ŷ(t) + h(y1 + 0.5hy2 ) = 0.4792 + 0.2 ∗ (1.459 + 0.1 ∗ 1.344) = 0.7978
t+h = 1.6
ŷ(t + h) = 0.7978
(d) Pour t = 1.6, ŷ(t) = 0.7978 on a

y = ŷ(t) 0.7978
t = 1.6 = 0.4986
y1 = 1 + y(1 + y) = 1 + 0.4986 ∗ (1 + 0.4986) = 1.747
y2 = 1t (y1 − y)(1 + 2y)
= (1/1.6) ∗ (1.747 − 0.4986) ∗ (1 + 2 ∗ 0.4986) = 1.558
ŷ(t) + h(y1 + 0.5hy2 )
= 0.7978 + 0.2 ∗ (1.747 + 0.1 ∗ 1.558) = 1.178
t + h = 1.8
ŷ(t + h) = 1.178
(e) Pour t = 1.8, ŷ(t) = 1.178 on a

80
ŷ(t)
y = t = 1.178/1.8 = 0.6544
y1 = 1 + y(1 + y) = 1 + 0.6544 ∗ (1 + 0.6544) = 2.082
1
y2 = t (y1 − y)(1 + 2y)
= (1/1.8) ∗ (2.082 − 0.6544) ∗ (1 + 2 ∗ 0.6544) = 1.832
ŷ(t) + h(y1 + 0.5hy2 )
= 1.178 + 0.2 ∗ (2.082 + 0.1 ∗ 1.832) = 1.631
t+h = 2
ŷ(t + h) = 1.631
On résume dans la table
t 1 1.2 1.4 1.6 1.8 2
ŷ(t) 0 0.22 0.4792 0.7978 1.178 1.631
2. une itération de la méthode de Taylor d’ordre 3 se traduit par :
pour h = 0.2, t et ŷ(t) connus,

ŷ(t + h) = ŷ(t) + hy 0 (t) + 21 h2 y 00 (t) + 16 h3 y 000 (t)

avec
y 0 (t) = f (t , ŷ(t))
y 00 (t) = f1 (t , ŷ(t) , y 0 (t))
y 000 (t) = f2 (t , ŷ(t) , y 0 (t) , y 00 (t))

Ainsi en arithmétique tronquée à 4 chiffres on a la séquence :

(a) Pour t = 1, ŷ(t) = 0 on a
ŷ(t)
y = t = 0/1 = 0
y1 = 1 + y(1 + y) = 1 + 0 ∗ (1 + 0) = 1
1
t1 = t = 1/1 = 1
u = 1 + 2y = 1 + 2 ∗ 0 = 1
u1 = y1 − y = 1 − 0 = 1
y2 = t1 u1 u = 1 ∗ 1 ∗ 1 = 1
y3 = t1 (y2 − 2t1 u1 )u + 2(t21 )(u21 )
= 1 ∗ (1 − 2 ∗ 1 ∗ 1) ∗ 1 + 2 ∗ (12 ) ∗ (12 )
= 1
ŷ(t + h) = ŷ(t) + h(y1 + h(0.5y2 + 0.1666hy3 ))
= 0 + 0.2 ∗ (1 + 0.2 ∗ (0.5 ∗ 1 + 0.03332 ∗ 1))
= 0.2212
t+h = 1.2
ŷ(t + h) = 0.2212
(b) Pour t = 1.2, ŷ(t) = 0.2212 on a

81
ŷ(t)
y = t = 0.2212/1.2 = 0.1843
y1 = 1 + y(1 + y) = 1 + 0.1843 ∗ (1 + 0.1843) = 1.218
1
t1 = t = 1/1.2 = 0.8333
u = 1 + 2y = 1 + 2 ∗ 0.1843 = 1.368
u1 = y1 − y = 1.218 − 0.1843 = 1.034
y2 = t1 u1 u = 0.8333 ∗ 1.034 ∗ 1.368 = 1.178
y3 = t1 (y2 − 2t1 u1 )u + 2(t21 )(u21 )
= 0.8333 ∗ (1 − 2 ∗ 0.8333 ∗ 1.034) ∗ 1.368 + 2 ∗ (0.83332 ) ∗ (1.0342 )
= 0.66
ŷ(t + h) = ŷ(t) + h(y1 + h(0.5y2 + 0.1666hy3 ))
= 0.2212 + 0.2 ∗ (1.218 + 0.2 ∗ (0.5 ∗ 1.178 + 0.03332 ∗ 0.66))
= 0.4892
t+h = 1.4
ŷ(t + h) = 0.4892
(c) Pour t = 1.4, ŷ(t) = 0.4892 on a
ŷ(t)
y = t = 0.4892/1.4 = 0.3494
y1 = 1 + y(1 + y) = 1 + 0.3494 ∗ (1 + 0.3494) = 1.471
1
t1 = t = 1/1.4 = 0.7142
u = 1 + 2y = 1 + 2 ∗ 0.3494 = 1.698
u1 = y1 − y = 1.471 − 0.3494 = 1.122
y2 = t1 u1 u = 0.7142 ∗ 1.122 ∗ 1.698 = 1.360
y3 = t1 (y2 − 2t1 u1 )u + 2(t21 )(u21 )
= 0.7142 ∗ (1.360 − 2 ∗ 0.7142 ∗ 1.122) ∗ 1.698 + 2 ∗ (0.71422 ) ∗ (1.1222 )
= 0.99
ŷ(t + h) = ŷ(t) + h(y1 + h(0.5y2 + 0.1666hy3 ))
= 0.4892 + 0.2 ∗ (1.471 + 0.2 ∗ (0.5 ∗ 1.360 + 0.03332 ∗ 0.99))
= 0.8118
t+h = 1.6
ŷ(t + h) = 0.8118
(d) Pour t = 1.6, ŷ(t) = 0.8118 on a

82
ŷ(t)
y = t = 0.8118/1.6 = 0.5073
y1 = 1 + y(1 + y) = 1 + 0.5073 ∗ (1 + 0.5073) = 1.764
1
t1 = t = 1/1.6 = 0.625
u = 1 + 2y = 1 + 2 ∗ 0.5073 = 2.014
u1 = y1 − y = 1.764 − 0.5073 = 1.257
y2 = t1 u1 u = 0.625 ∗ 1.257 ∗ 2.014 = 1.582
y3 = t1 (y2 − 2t1 u1 )u + 2(t21 )(u21 )
= 0.625 ∗ (1.582 − 2 ∗ 0.625 ∗ 1.257) ∗ 2.014 + 2 ∗ (0.6252 ) ∗ (1.2572 )
= 1.247
ŷ(t + h) = ŷ(t) + h(y1 + h(0.5y2 + 0.1666hy3 ))
= 0.8118 + 0.2 ∗ (1.764 + 0.2 ∗ (0.5 ∗ 1.582 + 0.03332 ∗ 1.247))
= 1.197
t+h = 1.8
ŷ(t + h) = 1.197
(e) Pour t = 1.8, ŷ(t) = 1.197 on a
ŷ(t)
y = t = 1.197/1.8 = 0.665
y1 = 1 + y(1 + y) = 1 + 0.665 ∗ (1 + 0.665) = 2.107
1
t1 = t = 1/1.8 = 0.5555
u = 1 + 2y = 1 + 2 ∗ 0.665 = 2.33
u1 = y1 − y = 2.107 − 0.665 = 1.442
y2 = t1 u1 u = 0.5555 ∗ 1.442 ∗ 2.33 = 1.866
y3 = t1 (y2 − 2t1 u1 )u + 2(t21 )(u21 )
= 0.5555 ∗ (1.866 − 2 ∗ 0.5555 ∗ 1.442) ∗ 2.33 + 2 ∗ (0.55552 ) ∗ (1.4422 )
= 1.623
ŷ(t + h) = ŷ(t) + h(y1 + h(0.5y2 + 0.1666hy3 ))
= 1.197 + 0.2 ∗ (2.107 + 0.2 ∗ (0.5 ∗ 1.866 + 0.03332 ∗ 1.623))
= 1.657
t+h = 2
ŷ(t + h) = 1.657
On la résume dans la table
t 1 1.2 1.4 1.6 1.8 2
ŷ(t) 0 0.2212 0.4892 0.8118 1.197 1.657

Exercice 40 Reprendre le traitement du problème de l’exemple ci-dessus En

utilisant une arithmétique décimale arrondie à 4 chiffres et la méthode de Taylor
à l’ordre 3

4.4 Méthodes de Runge Kutta

Pour l’élaboration des méthodes de Runge Kutta nous avons encore besoin
de l’outil essentiel qu’est la formule de Taylor. Dans ce contexte sa version à
deux variables sera d’une grande utilité.

83
4.4.1 Rappel de la formule de Taylor à deux variables
Théorème 4.4.1 Soit f une fonction réelle telle que f ainsi que ses dérivées
partielles jusqu’à l’ordre n + 1 sont continues sur

D = {(t, y)| a≤t≤b c ≤ y ≤ d}

Soit (t, y) , (t + h, y + k) ∈ D, alors il existe θ ∈]0, 1[ tel que

h i
f (t + h, y + k) = f (t, y) + h ∂f ∂t (t, y) + k ∂f
∂y (t, y)
h 2 2 2
i
+ 1
2 h2 ∂∂t2f (t, y) + 2hk ∂t∂y
∂ f
(t, y) + k 2 ∂∂yf2 (t, y)

+ ...
hP n
i
1 n ∂ f
+ n! j=0 Cnj hn−j k j ∂tn−j ∂y j (t, y)

hP n+1
i
1 n+1 j ∂ f
+ (n+1)! j=0 Cn+1 hn−j+1 k j ∂tn−j+1 ∂y j (t + θh, y + θk)

Preuve
Soient (t , y) ∈ D et (h , k) ∈ R I 2 fixés tel que (t + h , y + k) ∈ D.
Il suffit d’appliquer le théorème de Taylor Lagrange de la page 7 à la fonction
réelle g définie sur [0 , 1] par

s ∈ [0 , 1] 7→ g(s) = f (t + sh , y + sk)

qui est n + 1 fois continûment dérivable sur [0 , 1] comme composée de deux

fonctions n + 1 fois continûment différentiables.
Il existe donc θ ∈]0 , 1[ tel que :
n
X 1 dl g 1 dn+1 g
f (t + h , y + k) = g(1) = g(0) + (0) + (θ)
l! dsl (n + 1)! dsn+1
l=1

Et on montre en posant l’opérateur différentiel linéaire

∂ ∂
D = h +k
∂t ∂y
que
dg ∂f ∂f
(s) = h (t + sh , y + sk) + k (t + sh , y + sk) = Df (t + sh , y + sk)
dt ∂t ∂y
et pour tout entier naturel l ≤ n + 1 on a

dl g
(s) = Dl f (t + sh , y + sk)
dtl

84
Ainsi
dl g
(0) = Dl f (t , y)
dtl
dn+1 g
(θ) = Dn+1 f (t + θh , y + θk)
dtn+1
Grâce au théorème de Schwarz on peut appliquer la formule du binôme de
Newton au développement des puissance Dl de D
Remarque :
En condensé la formule de Taylor peut s’écrire :
n
X 1
f (t + h , y + k) = f (t , y) + Dl f (t , y) + Dn+1 f (t + θh , y + θk)
(n + 1)!
l=1

4.4.2 Méthodes de Runge Kutta d’ordre 2

Considérons le problème de Cauchy suivant
dy
dt = f (t, y) t ∈ [t0 , t0 + T ] avec T > 0
y(t0 ) = s
Le but des méthodes de Runge Kutta est de trouver une approximation du
problème sans se soumettre à la contrainte d’une évaluation des dérivées de la
fonction f tout en maintenant un niveau de précision équivalent à celui d’une
méthode de Taylor de même ordre.
La démarche de construction sera montrée dans le cas des méthodes d’ordre
2, en nous servant du développement de Taylor d’ordre 2 de la solution du
problème d Cauchy en un point ti d’un réseau de pas h.
Il existe ξi ∈]ti , ti+1 [ tel que
0
h2 00 h3 000
y(ti+1 ) = y(ti ) + hyh (ti , y(ti )) + 2 y (ti , y(ti )) + 6 y (ξ i , y(ξi )) i
h ∂f h ∂f
= y(ti ) + h f (ti , y(ti )) + (t
2 ∂t i , y(ti )) + 2 f ∂y (ti , y(ti ))
h3 000
+ 6 y (ξi , y(ξi ))
Cherchons l’expression entre crochets à l’aide de
af (ti , y(ti )) + bf (ti + α, y(ti ) + β)
Considérons un développement d’ordre 1 de f (ti + α, y(ti ) + β) au point
(ti , y(ti )) on a
∂f ∂f
f (ti + α, y(ti ) + β) = f (ti , y(ti )) + α (ti , y(ti )) + β (ti , y(ti )) + . . .
∂t ∂y
On obtient alors
af (ti , y(ti )) + bf (ti + α, y(ti ) + β)
= (a + b)f (ti , y(ti )) + bα ∂f ∂t (ti , y(t
∂f
i )) + bβ ∂y (ti , y(ti )) + . . .

≡ f (ti , y(ti )) + h2 ∂f h ∂f
∂t (ti , y(ti )) + 2 f ∂y (ti , y(ti )) + . . .

85
Par identification on a
a+b = 1
h
bα = 2
h
bβ = 2 f (ti , y(ti ))

En supposant b 6= 0 alors on a une famille à un paramètre

a = 1−b
h
α = 2b
h
β = 2b f (ti , y(ti ))

On peut citer quelques cas particuliers importants :

– b = 1 on a la méthode du point milieu qui est un Euler amélioré que
l’on peut écrire sous la forme :


 y0 donné en t0
Pour tout i




K1 = hf (ti , yi )


 K2 = hf (ti + 12 h , yi + 12 K1 )
yi+1 = yi + K2




ti+1 = ti + h


Ici l’erreur locale est en O(h3 )

– a = b = 12 on a la méthode Euler modifié que l’on peut écrire sous la
forme : 

 y0 donné en t0
 Pour tout i



K1 = hf (ti , yi )


 K2 = hf (ti + h , yi + K1 )
1
y i+1 = yi + 2 [K1 + K2 ]




ti+1 = ti + h


avec une erreur locale est en O(h3 )

– a = 14 et b = 34 on a la méthode Heun que l’on peut écrire sous la
forme :


 y0 donné en t0
Pour tout i




K1 = hf (ti , yi )


 K2 = hf (ti + 23 h , yi + 23 K1 )
yi+1 = yi + 14 [K1 + 3K2 ]




ti+1 = ti + h


avec une erreur locale est en O(h3 )

4.4.3 Méthodes de Runge Kutta d’ordre 4

Bien que l’objectif de ce cours n’est pas de mémoriser les algorithmes etudiés,
nous ne citons, ici dans la grande famille des méthodes de Runge Kutta d’ordre

86
4, que la plus utilisée dite méthode de Runge Kutta d’ordre 4 classique
et qui mérite un tel effort. Elle peut s’écrire sous la forme :


 y0 donné en t0
Pour tout i




K1 = hf (ti , yi )



K2 = hf (ti + 21 h , yi + 12 K1 )



 K3 = hf (ti + 21 h , yi + 12 K2 )
K4 = hf (ti + h , yi + K3 )



 1
y i+1 = yi + 6 [K1 + 2K2 + 2K3 + K4 ]




ti+1 = ti


avec une erreur locale est en O(h5 )

Remarque :
A chaque itération d’une méthode de Runge Kutta d’ordre 4, on procède à
quatre évaluations de la fonction f en vue d’avoir la précision du même ordre
que celle de la méthode de Taylor d’ordre 4 c’est à dire une erreur locale en
O(h5 ).
En comparaison la méthode de Taylor d’ordre 4 nécessite par itération, une
évaluation de f , celle de chacune des fonctions f (1) , f (2) et f (3) en passant
d’abord l’exercice de leur évaluation formelle, ce qui est très peu recommandé
dans un contexte numérique. Cette évaluation formelle est souvent source d’er-
reur.
L’algorithme de la méthode de Runge Kutta d’ordre 4 classique est :
real array (yi )n , (xi )n
function externe f
real h, t0 , tf in , y0 , K1 , K2 , K3 , K4
integer k
begin
h ← (tf in − t0 )/(n − 1)
t ← t0
u ← y0
x1 ← t
y1 ← u
Ecrire t, u
for k = 2 to n do
begin
K1 ← hf (t , u)
u1 ← u + 0.5K1
t1 ← t + 0.5h
K2 ← hf (t1 , u1 )
u1 ← u + 0.5K2
t1 ← t + 0.5h
K3 ← hf (t1 , u1 )
u1 ← u + K3
t1 ← t + h

87
K4 ← hf (t1 , u1 )
u ← u + 16 (K1 + 2K2 + 2K3 + K4 )
t←t+h
xk ← t
yk ← u
Ecrire t, u
end
end

real function f (t , y)
real t, y
integer k
begin
A definir
end

Remarque :
En théorie la structure d’une méthode de Runge Kutta d’ordre supérieur à 2
peut être très élaborée aussi bien implicite qu’explicite (amateur voir [4] !). Dans
la classe des structures explicites, les plus simples sont diagonales comme pour
le cas de la méthode de Runge Kutta d’ordre 4 classique. En général la structure
explicite est triangulaire inférieure.
Exemple
Pour le problème de Cauchy
dy


 = f (t , y) pour a ≤ t ≤ b
dt


y(a) = α
on peut envisager les traitement par l’une des deux méthodes de Runge Kutta
d’ordre 3 suivantes :
1. t et ȳ(t) connus, une itération est donnée par

ȳ(t + h) = ȳ(t) + 91 (2K1 + 3K2 + 4K3 )




K1 = hf (t, ȳ(t))

(M1)

 K2 = hf (t + 21 h, ȳ(t) + 12 K1 )
K3 = hf (t + 34 h, ȳ(t) + 34 K2 )


La structure de la méthode est diagonale : pour tout i > 1 seul Ki−1 est
utilisé pour calculer Ki .
2. t et ȳ(t) connus, une itération est donnée par
1


 ȳ(t + h) = ȳ(t) + 30 (5K1 + 9K2 + 16K3 )
K1 = hf (t, ȳ(t))

(M2)

 K2 = hf (t + 31 h, ȳ(t) + 13 K1 )
K3 = hf (t + 43 h, ȳ(t) − 163
K1 + 15
16 K2 )


88
La structure de la méthode est triangulaire inférieure : pour tout i > 1 on
utilise K1 à Ki−1 pour calculer Ki .
Remarque :
Outre l’utilisation concrète d’un algorithme de Runge Kutta nous devons être
capable d’analyser sa précision c’est à dire le comparer à la méthode de Taylor
du même ordre.
Exemple
Considérons le problème de Cauchy
dy y y 2


 = 1+ + pour 1 ≤ t ≤ 2
dt t t


y(1) = 0

1. En utilisant une arithmétique décimale arrondie à 4 chiffres calculer une

approximation de la solution sur une grille uniforme de pas h = 0.2 par la
méthode RK3 définie pour une EDO y 0 = f (t , y) par :
t et ȳ(t) connus, une itération est donnée par
1

 ȳ(t + h) = ȳ(t) + 9 (2K1 + 3K2 + 4K3 )

K1 = hf (t, ȳ(t))

(RK3)

 K2 = hf (t + 21 h, ȳ(t) + 12 K1 )
K3 = hf (t + 43 h, ȳ(t) + 34 K2 )


2. Montrons que la méthode RK3 restitue la solution donnée par la méthode

de Taylor d’ordre 3.
y y 2
Soit f (t , h) = 1 + +
t t
1. Calcul de solution
Une itération de la méthode RK3 se traduit par :
pour h = 0.2, t et ŷ(t) donnés

ȳ(t + h) = ȳ(t) + 91 (2K1 + 3K2 + 4K3 )

avec
K1 = hf (t, ȳ(t))
K2 = hf (t + 12 h, ȳ(t) + 12 K1 )
K3 = hf (t + 34 h, ȳ(t) + 34 K2 )

Ainsi en arithmétique décimale arrondie à 4 chiffres on a la séquence :

(a) Pour t = 1, ŷ(t) = 0 on a

89
t1 = t=1
y1 = ŷ(t) = 0
y1
u = t1 = 0/1 = 0
f (t , y1 ) = 1 + u(1 + u) = 1 + 0 ∗ (1 + 0) = 1
K1 = hf (t , y1 ) = 0.2 ∗ 1 = 0.2
t1 = t + 0.5h = 1 + 0.1 = 1.1
y1 = ŷ(t) + 0.5K1 = 0 + 0.5 ∗ 0.2 = 0.1
y1
u = t1 = 0.1/1.1 = 0.09091
f (t , y1 ) = 1 + u(1 + u) = 1 + 0.09091 ∗ (1 + 0.09091) = 1.099
K2 = hf (t , y1 ) = 0.2 ∗ 1.099 = 0.2198
t1 = t + 0.75h = 1 + 0.15 = 1.15
y1 = ŷ(t) + 0.5K2 = 0 + 0.75 ∗ 0.2198 = 0.1649
y1
u = t1 = 0.1649/1.15 = 0.1434
f (t , y1 ) = 1 + u(1 + u) = 1 + 0.1434 ∗ (1 + 0.1434) = 1.164
K3 = hf (t , y1 ) = 0.2 ∗ 1.164 = 0.2328
K = 2K1 + 3K2 + 4K3 = 2 ∗ 0.2 + 3 ∗ 0.2198 + 4 ∗ 0.2328 = 1.99
ŷ(t + h) = ŷ(t) + 19 K = 0 + (1/9) ∗ 1.99 = 0.2211
t+h = 1.2
ŷ(t + h) = 0.2211
(b) Pour t = 1.2, ŷ(t) = 0.2211 on a

t1 = t = 1.2
y1 = ŷ(t) = 0.2211
y1
u = t1 = 0.2211/1.2 = 0.1843
f (t , y1 ) = 1 + u(1 + u) = 1 + 0.1843 ∗ (1 + 0.1843) = 1.218
K1 = hf (t , y1 ) = 0.2 ∗ 1.218 = 0.2436
t1 = t + 0.5h = 1.2 + 0.1 = 1.3
y1 = ŷ(t) + 0.5K1 = 0.2211 + 0.5 ∗ 0.2436 = 0.3429
y1
u = t1 = 0.3429/1.3 = 0.2638
f (t , y1 ) = 1 + u(1 + u) = 1 + 0.2638 ∗ (1 + 0.2638) = 1.333
K2 = hf (t , y1 ) = 0.2 ∗ 1.333 = 0.2666
t1 = t + 0.75h = 1.2 + 0.15 = 1.35
y1 = ŷ(t) + 0.5K2 = 0.2211 + 0.75 ∗ 0.2666 = 0.4211
y1
u = t1 = 0.4211/1.35 = 0.3119
f (t , y1 ) = 1 + u(1 + u) = 1 + 0.3119 ∗ (1 + 0.3119) = 1.409
K3 = hf (t , y1 ) = 0.2 ∗ 1.409 = 0.2818
K = 2K1 + 3K2 + 4K3 = 2 ∗ 0.2436 + 3 ∗ 0.2666 + 4 ∗ 0.2818 = 2.414
ŷ(t + h) = ŷ(t) + 19 K = 0.2211 + (1/9) ∗ 2.414 = 0.4893
t+h = 1.4
ŷ(t + h) = 0.4893
(c) Pour t = 1.4, ŷ(t) = 0.4893 on a

90
t1 = t = 1.4
y1 = ŷ(t) = 0.4893
y1
u = t1 = 0.4893/1.4 = 0.3495
f (t , y1 ) = 1 + u(1 + u) = 1 + 0.3495 ∗ (1 + 0.3495) = 1.472
K1 = hf (t , y1 ) = 0.2 ∗ 1.472 = 0.2944
t1 = t + 0.5h = 1.4 + 0.1 = 1.5
y1 = ŷ(t) + 0.5K1 = 0.4893 + 0.5 ∗ 0.2944 = 0.6365
y1
u = t1 = 0.6365/1.5 = 0.4243
f (t , y1 ) = 1 + u(1 + u) = 1 + 0.4243 ∗ (1 + 0.4243) = 1.604
K2 = hf (t , y1 ) = 0.2 ∗ 1.604 = 0.3208
t1 = t + 0.75h = 1.4 + 0.15 = 1.55
y1 = ŷ(t) + 0.5K2 = 0.4893 + 0.75 ∗ 0.3208 = 0.7299
y1
u = t1 = 0.7299/1.55 = 0.4709
f (t , y1 ) = 1 + u(1 + u) = 1 + 0.4709 ∗ (1 + 0.4709) = 1.693
K3 = hf (t , y1 ) = 0.2 ∗ 1.693 = 0.3386
K = 2K1 + 3K2 + 4K3 = 2 ∗ 0.2944 + 3 ∗ 0.3208 + 4 ∗ 0.3386 = 2.905
ŷ(t + h) = ŷ(t) + 19 K = 0.4893 + (1/9) ∗ 2.905 = 0.8121
t+h = 1.6
ŷ(t + h) = 0.8121
(d) Pour t = 1.6, ŷ(t) = 0.8121 on a

t1 = t = 1.6
y1 = ŷ(t) = 0.8121
y1
u = t1 = 0.8121/1.6 = 0.5076
f (t , y1 ) = 1 + u(1 + u) = 1 + 0.5076 ∗ (1 + 0.5076) = 1.766
K1 = hf (t , y1 ) = 0.2 ∗ 1.766 = 0.3532
t1 = t + 0.5h = 1.6 + 0.1 = 1.7
y1 = ŷ(t) + 0.5K1 = 0.8121 + 0.5 ∗ 0.3532 = 0.9887
y1
u = t1 = 0.9887/1.7 = 0.5816
f (t , y1 ) = 1 + u(1 + u) = 1 + 0.5816 ∗ (1 + 0.5816) = 1.92
K2 = hf (t , y1 ) = 0.2 ∗ 1.92 = 0.384
t1 = t + 0.75h = 1.6 + 0.15 = 1.75
y1 = ŷ(t) + 0.5K2 = 0.8121 + 0.75 ∗ 0.384 = 1.100
y1
u = t1 = 1.100/1.75 = 0.6286
f (t , y1 ) = 1 + u(1 + u) = 1 + 0.6286 ∗ (1 + 0.6286) = 2.024
K3 = hf (t , y1 ) = 0.2 ∗ 2.024 = 0.4048
K = 2K1 + 3K2 + 4K3 = 2 ∗ 0.3532 + 3 ∗ 0.384 + 4 ∗ 0.4048 = 3.477
ŷ(t + h) = ŷ(t) + 19 K = 0.8121 + (1/9) ∗ 3.477 = 1.198
t+h = 1.8
ŷ(t + h) = 1.198
(e) Pour t = 1.8, ŷ(t) = 1.198 on a

91
t1 = t = 1.8
y1 = ŷ(t) = 1.198
y1
u = t1 = 1.198/1.8 = 0.6656
f (t , y1 ) = 1 + u(1 + u) = 1 + 0.6656 ∗ (1 + 0.6656) = 2.109
K1 = hf (t , y1 ) = 0.2 ∗ 2.109 = 0.4218
t1 = t + 0.5h = 1.8 + 0.1 = 1.9
y1 = ŷ(t) + 0.5K1 = 1.198 + 0.5 ∗ 0.4218 = 1.409
y1
u = t1 = 1.409/1.9 = 0.7416
f (t , y1 ) = 1 + u(1 + u) = 1 + 0.7416 ∗ (1 + 0.7416) = 2.292
K2 = hf (t , y1 ) = 0.2 ∗ 2.292 = 0.4584
t1 = t + 0.75h = 1.8 + 0.15 = 1.95
y1 = ŷ(t) + 0.5K2 = 1.198 + 0.75 ∗ 0.4584 = 1.542
y1
u = t1 = 1.542/1.95 = 0.7908
f (t , y1 ) = 1 + u(1 + u) = 1 + 0.7908 ∗ (1 + 0.7908) = 2.416
K3 = hf (t , y1 ) = 0.2 ∗ 2.416 = 0.4832
K = 2K1 + 3K2 + 4K3 = 2 ∗ 0.4218 + 3 ∗ 0.4584 + 4 ∗ 0.4832 = 4.152
ŷ(t + h) = ŷ(t) + 19 K = 1.198 + (1/9) ∗ 4.152 = 1.659
t+h = 2
ŷ(t + h) = 1.659
En résumé on a
t 1 1.2 1.4 1.6 1.8 2
ŷ(t) 0 0.2211 0.4893 0.8121 1.198 1.659
2. Analyse de précision
Pour l’EDO y 0 = f (t , y) et connaissant à la date t la solution y(t) :
La méthode de Taylor d’ordre 3 s’écrit
1 1
ytaylor (t + h) = y(t) + hy 0 (t) + hy 00 (t) + hy 000 (t)
2 6
avec  0
 y (t) = f (t , y(t)) ≡ f
y 00 (t) = ft + f fy (4.9)
 000
y (t) = ftt + 2f fty + f 2 fyy + (ft + f fy )fy
La méthode RK3 s’écrit
1
ŷ(t + h) = y(t) + (2K1 + 3K2 + 4K3 )
9
avec 
 K1 = hf (t , y(t)) ≡ hf
= hf t + 21 h , y(t) + 12 K1

K2
K3 = hf t + 34 h , y(t) + 34 K2


Le développement de Taylor avec y(t) ≡ y

f (t + α , y + β) = f (t , y) + (αft + βfy )
+ 21 α2 ftt + 2αβfty + β 2 fyy + . . .

92
f t + 12 h , y(t) + 12 K1 = f (t , y) + 12 hft + 12 K1 fy

+ 12 41 h2 ftt + 24 hK1 fty + 14 K12 fyy + . . .

K1 = hf ⇒ hK1 = h2 f et K12 = h2 f 2

f t + 12 h , y(t) + 21 K1 = f (t , y) + 12 hft + 21 hf fy

+ 12 14 h2 ftt + 24 h2 f fty + 14 h2 f 2 fyy + . . .

= f (t , y) + 12 h(ft + f fy )
+ 18 h2 ftt + 2f fty + f 2 fyy + O(h3 )

K2 = hf (t , y) + 12 h2 (ft + f fy )
+ 18 h3 ftt + 2f fty + f 2 fyy + O(h4 )

f t + 34 h , y(t) + 34 K2 = f (t , y) + 43 hft + 34 K2 fy

+ 12 16
9 2 18 9

h ftt + 16 hK2 fty + 16 K22 fyy + . . .

Avec f (t , y) ≡ f on a
K2 = hf + 12 h2 (ft + f fy )
1 3 2
4
+ 8 h ftt + 2f fty + f fyy + O(h )
⇓
hK2 = h2 f + 12 h3 (ft + f fy ) + O(h4 )
h2 K2 = h3 f + O(h4 )
K22 = h2 f 2 + h3 ft + f 2 fy + O(h4 )
hK22 = h3 f 2 + O(h4 )

f t + 43 h , y(t) + 34 K2 3
+ 34 K2 fy

= f+ 4 hft
1 9 2
+ 18 9 2

+ 2 16 h ftt 16 hK2 fty + 16 K2 fyy + ...

hf + 34 h2 ft + 34 hK2 fy

K3 =
1 9 3 18 2 9 2

+ 2 16 h ftt + 16 h K2 fty + 16 hK2 fyy + . ..
= hf + 34 h2 ft + 34 h2f + 12 h3 (ft + f fy ) fy
1 9 3 18 3 9 3 2
+ 2 16 h ftt + 16 h f fty + 16 h f fyy + . . .
3 2 3 3
= hf + 4 h (ft + f fy ) + 8 h(ft + f fy )fy
9 3 2 4
+ 32 h ftt + 2f fty + f fyy + O(h )

= (2 + 3 + 4)hf + 23 + 12
2
2K1 + 3K2 + 4K3 4 h (ft + f fy )
12 3
+ 8 h (f t 3+ f fy )fy
3 36 2
4
+ 8 + 32 h ftt + 2f fty + f fyy + O(h )
9 2
= 9hf + 2 h (ft + f fy )
+ 32 h3 (ft + f fy )fy + 32 h3 ftt + 2f fty + f 2 fyy + O(h4 )

= 9hf + 92 h2 (ft + f fy )
+ 32 h3 ftt + 2f fty + f 2 fyy + (ft + f fy )fy + O(h4 )

93
En utilisant l’équation 4.9 on obtient

2K1 + 3K2 + 4K3 = 9hf + 92 h2 (ft + f fy )

+ 32 h3 (ft + f fy )fy + 32 h3 ftt + 2f fty + f 2 fyy + O(h4 )

= 9hy 0 (t) + 29 h2 y 00 (t) + 32 h3 y 000 (t) + O(h4 )

Ainsi par la méthode RK3 on a

ŷ(t + h) = y(t) + 91 (2K1 + 3K2 + 4K3 )

y(t) + 91 9hy 0 (t) + 29 h2 y 00 (t) + 32 h3 y 000 (t) + O(h4 )

=
= y(t) + hy 0 (t) + 21 h2 y 00 (t) + 61 h3 y 000 (t) + O(h4 )
= ytaylor (t + h) + O(h4 )

Exercice 41 Reprendre les questions de l’exemple ci-dessus en utilisant une

arithmétique décimale tronquée à 4 chiffres et une méthode RK3 définie pour
une EDO y 0 = f (t , y) par :
t et ȳ(t) connus, une itération est donnée par
1


 ȳ(t + h) = ȳ(t) + 30 (5K1 + 9K2 + 16K3 )
K1 = hf (t, ȳ(t))

(RK3)

 K2 = hf (t + 13 h, ȳ(t) + 13 K1 )
K3 = hf (t + 43 h, ȳ(t) − 16
3
K1 + 15
16 K2 )


4.5 Méthodes multipas

Considérons à nouveau le problème de Cauchy suivant
dy
dt = f (t, y) t ∈ [a, b]
y(a) = s

L’objectif ici est, lors de l’approximation de la solution de ce problème sur un

réseau t0 = a < t1 < t2 < . . . < ti < ti+1 < . . . < tN = b, d’utiliser les
informations déjà obtenues.
En intégrant sur [ti , ti+1 ] ⊂ [a, b] on a
Z ti+1
y(ti+1 ) = y(ti ) + f (t , y(t))dt
ti

Si P (t) est un interpolant de f (t , y(t)) par certains des points de l’ensemble

(t0 , y0 ), (t1 , y1 ), . . . , (ti , yi )

où yi est une approximation de y(ti ), alors

Z ti+1
y(ti+1 ) ≈ yi + P (t)dt
ti

Supposons que pour calculer yi+1 on ait besoin de connaı̂tre yi−k+1 , yi−k+2 , . . . , yi−1 , yi
les k valeurs précédentes, on dit qu’on a une méthode multipas ou plus précisement

94
une méthode à k pas. Pour les méthodes de Runge Kutta et de Taylor on n’a
eu besoin que de yi ; ce sont des méthodes à un pas.
On a deux classes de méthodes multipas selon que l’interpolation de f (t , y(t))
utilise ou non le point (ti+1 , yi+1 ) qui est, pour l’heure, inconnu.
– Les méthodes explicites n’utilisent pas le point (ti+1 , yi+1 )
– Les méthodes implicites par contre qui utilisent le point (ti+1 , yi+1 )
Pour les méthodes à k pas, l’interploant P est un polynôme de degré inférieur
k et
 P
k
j=1 cj fi−j+1 explicite
Z ti+1 

P (t)dt = avec fl = f (tl , yl )∀l
ti  Pk c f

implicite
j=1 j i−j+2

où c1 , c2 , . . . , ck sont déterminés pour que la formule soit exacte pour tout
polynôme de degré inférieur k
En observant que sur un réseau régulier de pas h et que
Z ti+1 Z 1
P (t)dt = h P (ti + hs)ds
ti 0

alors on peut déterminer c1 , c2 , . . . , ck pour que pour tout polynôme F de degré

inférieur k on ait

Z 1  c1 F (0) + c2 F (−1) + . . . + ck F (1 − k) explicite
F (s)ds =
0 
c1 F (1) + c2 F (0) + . . . + ck F (2 − k) implicite

Ainsi pour une méthode explicite c1 , c2 , . . . , ck sont solutions du système

linéaire
 R1 Pk
 1 ≡ 0 ds
 = j=1 cj

1
R1 Pk i−j
si−1 ds =


i ≡ 0 j=1 cj (1 − j) ∀i = 2, . . . , k

Pour une méthode implicite c1 , c2 , . . . , ck sont solutions du système linéaire

 R1 Pk
 1 ≡ 0 ds
 = j=1 cj

1
R1 Pk i−j
si−1 ds =


i ≡ 0 j=1 cj (2 − j) ∀i = 2, . . . , k

Nous citerons quelques exemples dans la familles des méthodes d’Adams.

Il s’agit des méthodes d’Adams-Bashford pour les explicites et des méthodes
d’Adams-Moulton pour les implicites.

4.5.1 Méthodes d’Adams-Bashford

Théorème 4.5.1 Soit ti = a + ih pour i = 0 . . . N un réseau régulier de pas
h = b−a
N . Soit yi une approximation par la méthode d’Adams-Bashford à 2

95
pas de y(ti ), solution du problème de Cauchy. La méthode d’Adams-Bashford
à 2 pas est donnée par :


 y0 = α0
y1 = α1

 ∀i = 1, . . . , N − 1

yi+1 = yi + 12 h [3f (ti , yi ) − f (ti−1 , yi−1 )]


5 3 (3)
Si la solution est 3 fois dérivable sur [a, b] alors l’erreur locale est 12 h y (ξ)
avec ξ ∈]ti−1 , ti+1 [.

Théorème 4.5.2 Soit ti = a + ih pour i = 0 . . . N un réseau régulier de pas

h = b−a
N . Soit yi une approximation par la méthode d’Adams-Bashford à 3
pas de y(ti ), solution du problème de Cauchy. La méthode d’Adams-Bashford
à 3 pas est donnée par :


 y0 = α0
 y1 = α1


y2 = α2
∀i = 2, . . . , N − 1



 1
yi+1 = yi + 12 h [23f (ti , yi ) − 16f (ti−1 , yi−1 ) + 5f (ti−2 , yi−2 )]


3 4 (4)
Si la solution est 4 fois dérivable sur [a, b] alors l’erreur locale est 8h y (ξ)
avec ξ ∈]ti−2 , ti+1 [.

L’algorithme de la méthode d’Adams-Bashford à 3 pas est :

real array (yi )n , (xi )n

function externe f
real h, t0 , tf in , Km0 , Km1 , Km2
integer k
begin
h ← (tf in − t0 )/(n − 1)
t ← t0
u ← α0
x1 ← t
y1 ← u
Ecrire t, u
Km2 ← hf (t , u)
t←t+h
u ← α1
x2 ← t
y2 ← u
Ecrire t, u
Km1 ← hf (t , u)
t←t+h
u ← α2

96
x3 ← t
y3 ← u
Ecrire t, u
for k = 4 to n do
begin
Km0 ← hf (t , u)
1
u ← u + 12 (23Km0 − 16Km1 + 5Km2 )
Km2 ← Km1
Km1 ← Km0
t←t+h
xk ← t
yk ← u
Ecrire t, u
end
end

real function f (t , y)
real t, y
integer k
begin
A definir
end

Exemple
Considérons le problème de Cauchy
dy y y 2


 = 1+ + pour 1 ≤ t ≤ 2
dt t t


y(1) = 0
En utilisant une arithmétique décimale arrondie à 4 chiffres calculer une ap-
proximation de y(1.6) sur une grille uniforme de pas h = 0.2 par la méthode
d’Adams-Bashford à 3 pas en amorçant par deux itérations de la méthode RK3
définie pour une EDO y 0 = f (t , y) par :
t et ȳ(t) connus, une itération est donnée par
1

 ȳ(t + h) = ȳ(t) + 9 (2K1 + 3K2 + 4K3 )

K1 = hf (t, ȳ(t))

(RK3)

 K2 = hf (t + 21 h, ȳ(t) + 12 K1 )
K3 = hf (t + 43 h, ȳ(t) + 34 K2 )


y y 2
Soit f (t , h) = 1 + +
t t
Une itération de la méthode d’Adams-Bashford à 3 pas se traduit par :
pour h = 0.2, t, ŷ(t − 2h), ŷ(t − h) et ŷ(t) donnés
1
ŷ(t+h) = ŷ(t)+ h [23f (t , ŷ(t)) − 16f (t − h , ŷ(t − h)) + 5f (t − 2h , ŷ(t − 2h))]
12

97
Ainsi pour calculer ŷ(1.6) l’approximation de y(1.6) par la la méthode d’Adams-
Bashford à 3 pas, on a besoin de :

 ŷ(1.6 − 3h) = ŷ(1) = y(1) donné
ŷ(1.6 − 2h) = ŷ(1.2) à calculer
ŷ(1.6 − h) = ŷ(1.4) à calculer


L’armorceur RK3 va donc nous fournir ŷ(1.2) et ŷ(1.4) avant la mise en œuvre
de la d’Adams-Bashford à 3 pas.
Nous raccourcissons ici le développement en empruntant ces résultats à un
exemple précédent où nous avons travaillé avec la même machine arithmétique
et le même algorithme numérique. On a donc ŷ(1.2) = 0.2211 et ŷ(1.4) = 0.4893
1. Calcul de solution approchée ŷ(1.6)
Pour t = 1.4
t1 = = t − 2h = 1.4 − 0.4 = 1
y1 = ŷ(t − 2h) = ŷ(1) = 0
u = yt11 = 0/1 = 0
Km2 = f (t1 , y1 ) = 1 + u ∗ (1 + u) = 1 + 0 ∗ (1 + 0) = 1
t1 = = t − h = 1.4 − 0.2 = 1.2
y1 = ŷ(t − h) = ŷ(1.2) = 0.2211
u = yt11 = 0.2211/1.2 = 0.1843
Km1 = f (t1 , y1 ) = 1 + u ∗ (1 + u) = 1 + 0.1843 ∗ (1 + 0.1843) = 1.218
t1 = = t = 1.4
y1 = ŷ(t) = ŷ(1.4) = 0.4893
u = yt11 = 0.4893/1.4 = 0.3495
Km0 = f (t1 , y1 ) = 1 + u ∗ (1 + u) = 1 + 0.3495 ∗ (1 + 0.3495) = 1.472
K = h(23Km0 − 16Km1 + 5Km2 ) = 0.2 ∗ (23 ∗ 1.472 − 16 ∗ 1.218 + 5 ∗ 1) = 3.874
1
ŷ(t + h) = ŷ(t) + 12 K = 0.4893 + (1/12) ∗ 3.874 = 0.8121
t + h = 1.6
ŷ(t + h) = 0.8121

On a ŷ(1.6) = 0.8121

Exercice 42 Reprendre l’exemple ci-dessus et achever le calcul de la solution

sur l’intervalle [1 , 2].

Exercice 43 Reprendre la question de l’exemple ci-dessus en utilisant une arithmétique

décimale tronquée à 4 chiffres et la méthode d’Adams-Bashford à 2 pas amorcée
par une itération de la méthode de Taylor d’ordre 2.

4.5.2 Méthodes d’Adams-Moulton

Théorème 4.5.3 Soit ti = a + ih pour i = 0 . . . N un réseau régulier de pas
h = b−a
N . Soit yi une approximation par la méthode d’Adams-Moulton à 2

98
pas de y(ti ), solution du problème de Cauchy. La méthode d’Adams-Moulton
à 2 pas est donnée par :


 y0 = α0
y1 = α1

 ∀i = 1, . . . , N − 1
 1
yi+1 = yi + 12 h [5f (ti+1 , yi+1 ) + 8f (ti , yi ) − f (ti−1 , yi−1 )]


1 4 (4)
Si la solution est 4 fois dérivable sur [a, b] alors l’erreur locale est − 24 h y (ξ)
avec ξ ∈]ti−1 , ti+1 [.

Théorème 4.5.4 Soit ti = a + ih pour i = 0 . . . N un réseau régulier de pas

h = b−aN . Soit yi une approximation par la méthode d’Adams-Moulton à 3
pas de y(ti ), solution du problème de Cauchy. La méthode d’Adams-Moulton
à 3 pas est donnée par :


 y0 = α0
 y1 = α1


y2 = α2
∀i = 2, . . . , N − 1



 1
yi+1 = yi + 24 h [9f (ti+1 , yi+1 ) + 19f (ti , yi ) − 5f (ti−1 , yi−1 ) + f (ti−2 , yi−2 )]


19 5 (5)
Si la solution est 5 fois dérivable sur [a, b] alors l’erreur locale est − 720 h y (ξ)
avec ξ ∈]ti−2 , ti+1 [.

Exemple
Considérons le problème de Cauchy
dy y y 2


 = 1+ + pour 1 ≤ t ≤ 2
dt t t


y(1) = 0

En utilisant une arithmétique décimale tronquée à 4 chiffres calculer une ap-

proximation de y(1.4) sur une grille uniforme de pas h = 0.2 par la méthode
d’Adams-Moulton à 2 pas amorcée par la méthode de Taylor d’ordre 3.
y y 2
Soit f (t , h) = 1 + +
t t
Une itération de la méthode d’Adams-Moulton à 2 pas se traduit par :
pour h = 0.2, t, ŷ(t − h) et ŷ(t) donnés
1
ŷ(t + h) = ŷ(t) + 12 h [5f (t + h , ŷ(t + h)) + 8(t , ŷ(t)) − f (t − h , ŷ(t − h))]
m
5
0 = −ŷ(t + h) + 12 hf (t + h , ŷ(t + h))
1
+ ŷ(t) + 12 h [+8f (t , ŷ(t)) − f (t − h , ŷ(t − h))]

99
Ceci équivaut à l’équation
5
0 = −ŷ(t + h) + 12 hf (t + h , ŷ(t + h))
1
+ ŷ(t) + 12 h [+8Km0 − Km1 ]
m
5h 1 2 5h 1
0 = 12 (t+h)2 ŷ (t + h) + 12 (t+h) − 1 ŷ(t + h)
5h h
+ 12 + ŷ(t) + 12 [+8Km0 − Km1 ]

où Km0 = f (t , ŷ(t)) et Km1 = f (t − h , ŷ(t − h))

Ainsi pour calculer ŷ(1.4) l’approximation de y(1.4) par la la méthode d’Adams-
Moulton à 2 pas, on a besoin de :

ŷ(1.4 − 2h) = ŷ(1) = y(1) donné
ŷ(1.4 − h) = ŷ(1.2) à calculer par Taylor d’ordre 3
Nous empruntons ŷ(1.2) = 0.2212 déjà calculé dans l’exemple de la sec-
tion 4.3.3. Le reste de calcul va se dérouler selon la séquence :
1. Pour t = 1.2 et ŷ(t) = 0.2212
t1 = = t − h = 1.2 − 0.2 = 1
y1 = ŷ(t − h) = ŷ(1) = 0
u = yt11 = 0/1 = 0
Km1 = f (t1 , y1 ) = 1 + u ∗ (1 + u) = 1 + 0 ∗ (1 + 0) = 1
t1 = = t = 1.2
y1 = ŷ(t) = ŷ(1.2) = 0.2212
u = yt11 = 0.2212/1.2 = 0.1843
Km0 = f (t1 , y1 ) = 1 + u ∗ (1 + u) = 1 + 0.1843 ∗ (1 + 0.1843) = 1.218
K = 8Km0 − Km1 = 8 ∗ 1.218 − 1 = 8.744
Ici commence le traitement de l’équation ay 2 + by + c = 0
1
c1 = ŷ(t) + 12 hK = 0.2212 + 0.01666 ∗ 8.744 = 0.3668
5
c = 12 + c1 = 0.4166 + 0.3668 = 0.7834
tph = t + h = 1.2 + 0.2 = 1.4
5 1
b1 = 12 tph = 0.4166/1.4 = 0.2975
b1
a = tph = 0.2975/1.4 = 0.2125
b = b1 − 1 = 0.2975 − 1 = −0.703
∆ = b2 − 4ac = (−0.703)2 − 4 ∗ 0.2125 ∗ 0.7834 = −0.1716
Problème numériquement non solvable !
ŷ(t + h) = ???
Remarque :
En pratique une méthode multipas ne s’utilise pas seule car elle n’est pas auto-
amorçable. Il est nécessaire de faire appel à une méthode à un pas de précision
équivalente pour démarrer le processus. Les méthodes de Runge Kutta sont
souvent utilisées à cette fin.
Remarque :
Pour contourner les difficultés de résolution d’équation souvent non linéaire, une
méthode multipas implicite est souvent utilisée en association avec une méthode
explicite ; on obtient une méthode prédicteur-correcteur.

100
4.6 Méthodes prédicteur-correcteur
Dans une méthode prédicteur-correcteur on utilise une combinaison d’une
méthode explicite, le prédicteur et d’une méthode implicite de même ordre , le
correcteur.
La méthode explicite prédit une approximation que vient corriger la méthode
implicite. On peut citer les exemples suivantes.
Exemple
Dans ce exemple on utilise la méthode d’Euler comme prédicteur et la méthode
d’Euler implicite est la correctrice. Cela donne :

 y0 = α0

∀i = 0, . . . , N − 1


 ȳi+1 = yi + hf (ti , yi )
yi+1 = yi + hf (ti+1 , ȳi+1 )


qui est en fait une méthode de Runge Kutta d’ordre 2.

Exemple
En combinant la méthode d’Adams-Bashford à 3 pas à la méthode d’Adams-
Moulton à 2 pas, on obtient :


 y0 = α 0
y 1 = α1




y2 = α 2


 ∀i = 2, . . . , N − 1
1
ȳ i+1 = yi + 12 h [23f (ti , yi ) − 16f (ti−1 , yi−1 ) + 5f (ti−2 , yi−2 )]



 1
yi+1 = yi + 12 h [5f (ti+1 , ȳi+1 ) + 8f (ti , yi ) − f (ti−1 , yi−1 )]


Exercice 44 Analyser la précision de la méthode de l’exemple ci-dessus en

indiquant l’ordre de la méthode de Taylor qu’elle restitue. On pourra à cet effet
utiliser un développement de Taylor d’ordre 4.

4.7 Notion de stabilité de solution des ODE

Nous ouvrons une petite parenthèse pour soulever une question qui peut
être source d’embarras pour le numéricien si un minimum d’analyse théorique
ne soutend son travail ; c’est la question de la stabilité. Elle revêt deux aspects :
– le premier aspect est celui lié à la nature et au paramétrage de la méthode
numérique utilisée, avec comme corollaire, l’amplification des erreurs d’ar-
rondis
– le deuxième aspect est inhérent au comportement explosif propre à la
solution recherchée.
Nous ne regarderons ici que le dernier aspect.
Consiérons le problème de valeur initiale
dy
dt = f (t, y) avec t ∈ [a, b]
y(a) = s

101
La solution exacte y(t) dépend de la valeur initiale s donc nous pouvons l’ecrire
sous la forme y(t , s). Ainsi l’équation différentielle donne une famille de solu-
tions dépendant du paramètre s. Comme tout calcul sur ordinateur est toujours
entaché d’une certaine erreur d’arrondis, on est endroit de se demander ce que
adviendrait de la solution exacte si la valeur initiale est affectée d’une erreur
d’arrondis ?
Par un développement en serie de Taylor si s devient s + δs, on a

∂ 1 ∂2
y(t , s + δs) = y(t , s) + δs y(t , s) + δs2 2 y(t , s) + . . .
∂s 2 ∂s
ainsi
∂ 1 ∂2 ∂
|y(t , s + δs) − y(t , s)| = δs y(t , s) + δs2 2 y(t , s) + . . . ≈ δs y(t , s)
∂s 2 ∂s ∂s

Comme y(t , s) est solution de l’équation différentielle alors on a

∂ ∂ ∂ ∂
y(t , s) = f (t , y(t, s)) ⇒ y(t , s) = f (t , y(t, s))
∂t ∂s ∂t ∂s
Or
∂ ∂ ∂t ∂
f (t , y(t, s)) = fy (t , y(t, s)) y(t, s)+ft (t , y(t, s)) = fy (t , y(t, s)) y(t, s)
∂s ∂s ∂s ∂s
Donc en posant
Z t
∂
u(t) = y(t, s) , q(t) = fy (t , y(t, s)) et Q(t) = q(r)dr
∂s a

on obtient l’équation différentielle

0
u = qu dont la solution est u(t) = ceQ(t)

Supposons qu’il existe γ > 0 tel que 0 < γ < q(t) pour tout t > a. Alors on
a Z t Z t
Q(t) = q(r)dr > γdr = γ(t − a)
a a
puis
|y(t , s + δs) − y(t , s)| ≈ |δs u(t)| > |δs c| eγ(t−a) → ∞
lorsque t → ∞. On observe une divergence des solutions
Supposons qu’il existe γ > 0 tel que q(t) < −γ < 0 pour tout t > a. Alors
on a Z t Z t
Q(t) = q(r)dr < − γdr = −γ(t − a)
a a
puis
|y(t , s + δs) − y(t , s)| ≈ |δs u(t)| < |δs c| e−γ(t−a) → 0

102
lorsque t → ∞. On a convergence des deux solutions
Exemple
Considérons le problème de Cauchy x0 = x avec x(a) = s. On a sur ]a, ∞[, une
famille de solutions x(t) = se(t−a) . Ces solutions divergent.
Exemple
Considérons le problème de Cauchy x0 = −x avec x(a) = s. On a sur ]a, ∞[,
une famille de solutions x(t) = se(a−t) . Ces solutions convergent.

4.8 Les systèmes d’équations différentielles

Dans la plupart des applications, on est amené à résoudre des systèmes
différentiels d’ordre 1 ou des équations qui s’y ramènent. Du fait de l’extension
de la formule de Taylor aux fonctions vectorielles, toutes les méthodes mises au
point pour les problèmes scalaires s’étendent naturellement au cas vectoriel en
travaillant composante par composante.
Par exemple la méthode de Runge Kutta d’ordre 4 appliquée au problème
vectoriel  du1
 dt = f1 (t , u1 , u2 , . . . , um )
du2

= f2 (t , u1 , u2 , . . . , um )



 dt

 ..
.




 dum
dt = fm (t , u1 , u2 , . . . , um )

 u 1 = α1
u = α2



 2



 .
..



um = αm
pourra s’écrire :


 u0,1 = α1
u0,2 = α2





 ..
.




u0,m = αm




Pour tout i



 
K = hfj (ti , ui,1 , . . . , ui,j , . . . , ui,m )

1j

 
 
∀j = 1 . . . m

 


 
= hfj (ti + 21 h , ui,1 + 12 K11 , . . . , ui,j + 12 K1j , . . . , ui,m + 12 K1m )

K2j




∀j = 1 . . . m
 


 

 K3j = hfj (ti + 12 h , ui,1 + 12 K21 , . . . , ui,j + 12 K2j , . . . , ui,m + 12 K2m )

 



∀j = 1 . . . m




K4j = hfj (ti + h , ui,1 + K21 , . . . , ui,j + K3j , . . . , ui,m + K3m )

 


 

∀j = 1 . . . m

 


 

ui+1,j = ui,j + 61 [K1j + 2K2j + 2K3j + K4j ]

 


 

∀j = 1 . . . m

 


 

ti+1 = ti + h
 

103
4.9 Travaux Pratiques
Cette section réunit des exercices de programmation dans un environnement
de travail comme le SciLab.

Exercice 45 Considérons le modèle ”Proie-Prédateur” de dynamique de popu-

lation où les grandeurs x1 (t) et x2 (t) représentent respectivement le nombre de
proies et le nombre de prédateurs à une date t :
 d
 dt x1 (t) = 3x1 (t) − 0.002x1 (t)x2 (t)
0≤t≤4
 d
dt x2 (t) = 0.0006x1 (t)x2 (t) − 0.5x2 (t)

1)Résoudre le problème de Cauchy pour

x1 (0) = 1000
x2 (0) = 500

2) existe - t-il une solution stable pour ce modèle de population ? Si c’est le

cas ce serait pour quelles valeur de x1 et x2 .

Exercice 46 Considérons le modèle de dynamique de population mettant en

compétition deux espèces sur la même ressource alimentaire où les grandeurs
x1 (t) et x2 (t) représentent le nombre de chaque espèce à une date t :
 d
 dt x1 (t) = x1 (t)[4 − 0.0003x1 (t) − 0.0004x2 (t)]
0≤t≤4
 d
dt x 1 (t) = x 2 (t)[2 − 0.0002x 1 (t) − 0.0001x 2 (t)]
1)Résoudre le problème de Cauchy pour

x1 (0) = 10000
x2 (0) = 10000

2) existe - t-il une solution stable pour se modèle de population ? Si c’est le

cas ce serait pour quelles valeur de x1 et x2 .

104
Chapitre 5

Rappel et complément
d’algèbre linéaire

La plupart des démontrations de ce rappel seront laissées en exercice.

5.1 Principales définitions d’algèbre linéaire

Dans ce cours, le corps des scalaires qu’on notera K
I désignera le corps des
nombres réels R
I ou le corps des nombres complexes CI.

5.1.1 Vecteurs
Définition
I N l’ensemble des vecteurs x de composantes x1 , x2 , . . . , xN ∈
On désignera par K
K
I .
I N on définit deux opérations
Sur K

addition z =x+y définie par zi = xi + yi 1 ≤ i ≤ N
multiplication par un scalaire z = λx définie par zi = λxi 1 ≤ i ≤ N

qui en font un espace vectoriel sur K

I .
En notation matricielle le vecteur x sera représenté par le vecteur colonne
 
x1
 x2 
x= .
 
 ..


xN

et on notera xt et x∗ les vecteurs lignes suivantes :

xt = (x1 , x2 , . . . , xN ) x∗ = (x̄1 , x̄2 , . . . , x̄N )

105
On remarque qu’en posant
     
1 0 0
 0   1   0 
e1 =  .  , e2 =  .  , . . . , eN = 
     
..
 ..   ..

  . 
0 0 1
N
X
x= xi ei
i=1
et l’ensemble des N vecteurs (e1 , e2 , . . . , eN ) est appelé base canonique S de
I N.
K
On dit qu’un sous ensemble ∅ =6 F ⊂E =K I N est sous espace vectoriel de
E si
∀x, y ∈ F x+y ∈F
∀x ∈ F, ∀λ ∈ K I λx ∈ F

Independance linéaire, base

On dit que k vecteurs v1 , v2 , . . . , vk sont linéairement indépendants si
k
X
λi vi = 0 (où λi ∈ K
I ) ⇒ λ1 = λ2 = . . . = λk = 0
i=1

On dit qu’un sous espace vectoriel F ⊂ E = K I est engendré par l’ensemble

des vecteurs (v1 , v2 , . . . , vk ) si
( k
)
X
F = x∈E | x= λi vi où λi ∈ K
I = Comb(v1 , v2 , . . . , vk )
i=1

On appelle base d’un espace vectoriel E toute suite de vecteurs linéairement

indépendants qui engendre E.

Produit scalaire ou hermitien

L’application (.|.) : E × E → KI définie par
PN
(x|y) = y t x = xt y = i=1 xi yi si K
I =R
I
PN
(x|y) = y ∗ x = x∗ y = i=1 xi ȳi si K
I =C
I
est appelé produit scalaire.
Au produit scalaire on associe la norme euclidienne définie par
v
uN
1 uX
kxk = (x|x) = t
2
2
|x | i
i=1

On a l’inégalité de Schwarz
|(x|y)| ≤ kxk2 kyk2

106
Orthogonalité
Soit E un espace vectoriel sur K
I muni d’un produit scalaire
Deux vecteurs x et y sont dits orthogonaux si (x|y) = 0. Dans ce cas on a
Théorème de Pythagore
2 2 2
kx + yk2 = kxk2 + kyk2

On dit qu’un vecteur x est orthogonal à une partie F de E et on note x ⊥ F

si pour tout y ∈ F on a (x|y) = 0
Soit F un sous espace vectoriel donné de E, on définit son sous espace
orthogonal F ⊥ par

F ⊥ = {x ∈ E tel que ∀y ∈ F (x|y) = 0}

Remarque :
F étant sous espace vectoriel de E, on a E = F ⊕F ⊥ . Un ensemble {v1 , v2 , . . . , vk }
de vecteurs de l’espace vectoriel E est dit orthonormal si

(vi |vj ) = δij , 1 ≤ i, j ≤ k

où δij est le symbole de Kronecker défini par

1 si i = j
δij =
0 si i 6= j

Exercice 47 Méthode d’orthogonalisation de Gram-Schmidt

I N . On
Soit une suite (f1 , f2 , . . . , fk ) de vecteurs linéairement indépendants de K
construit la suite (p1 , p2 , . . . , pk ) de vecteurs par :
(
p1 = f1
Pi−1 (p |f )
pi = fi − j=1 kpj ki2 pj ∀i ≤ k
j 2

Montrer que pour tout l < i ≤ k on a pl ⊥ pi et

Comb(f1 , f2 , . . . , fk ) = Comb(p1 , p2 , . . . , pk )

5.1.2 Matrices
Définition
Soient E et F deux espaces vectoriels sur le même corps KI .
Une application f de E dans F est dite linéaire si on a les deux propriétés
suivantes :
∀x, y ∈ E f (x + y) = f (x) + f (y)
∀x ∈ E, ∀λ ∈ K
I f (λx) = λf (x)
On notera L(E, F ) l’ensemble des applications linéaires de E dans F .
L(E) = L(E, E) est l’ensemble des endomorphismes de E

107
Dans L(E, F ) on définit deux opérations :

f + g par (f + g)(x) = f (x) + g(x) ∀x ∈ E
λf par (λf )(x) = λf (x) ∀x ∈ E, ∀λ ∈ K
I

L(E, F ) a alors une structure d’espace vectoriel sur K I

En rapportant E à une base E = (ei )1≤i≤N et F à une base F = (fj )1≤j≤M
une application linéaire f ∈ L(E, F ) est définie de façon unique par la donnée
des vecteurs :
 
a1j
XM  a2j 
f (ej ) = aij fi =  .  1≤j≤N
 
 .. 
i=1
aM j

Ainsi f est représenté par la matrice à M lignes et à N colonnes :

 
a11 a12 . . . a1N
 a21 a22 . . . a2N 
A= .  = (aij )
 
.. ..
 .. . . 
aM 1 aM 2 . . . aM N
PN
Pour x = j=1 xj ej on exprime y = f (x) dans la base F par
       
y1 a11 a12 a1N
 y2   a21   a22   a2N 
 = x1   + x2   + . . . + xN 
       
 .. .. .. .. 
 .   .   .   . 
yM aM 1 aM 2 aM N

Symboliquement on ecrit y = Ax
La matrice A représentant l’application linéaire f ∈ L(E, F ) relativement
N
aux bases (E, F) est la donnée d’un élément de K I M I M ×N
≡K
Soit
     
a11 a12 a1N
 a21   a22   a2N 
a1 =  . a2 =  . , . . . , aN =  . ∈K I M
     
 ..  ..  ..
 
  
aM 1 aM 2 aM N

on peut écrire A = (a1 , a2 , . . . , aN )

On appelle image de A, le sous espace vectoriel
n o
ImA = y ∈ K I M , ∃x ∈ K I N tel que Ax = y = Comb(a1 , a2 , . . . , aN )

rg(A) = dim(ImA) est le rang de la matrice A

108
Remarque :
L’équation Ax = b admet au moins une solution si et seulement si rg(A) = M
On appelle noyau de A, le sous espace vectoriel
n o
KerA = x ∈ K I N , tel que Ax = 0
Remarque :
L’équation Ax = b admet au plus une solution si KerA = {0}
I M ×N on a
Pour toute matrice A ∈ K
rg(A) + dim(KerA) = N

5.1.3 Operations sur les matrices

Espace vectoriel de matrices, transposition, adjonction
La structure d’espace vectoriel de L(E, F ) sur K I induit sur KI M ×N une
structure d’espace vectoriel sur KI , identique à sa structure naturelle.
Soit la matrice A ∈ R I M ×N , on appelle matrice transposée de A l’unique
matrice At ∈ RI N ×M définie par :
(Ax|y)M = x|At y N ∀x ∈ R I N , ∀y ∈ R
IM

I M ×N et B = (bij ) ∈ R
Exercice 48 Montrer que si A = (aij ) ∈ R I N ×M ,
B = At ⇔ ∀(i, j) bij = aji 1 ≤ i ≤ N, 1 ≤ j ≤ M
Soit la matrice A ∈ C I M ×N , on appelle matrice adjointe de A l’unique
∗ N ×M
matrice A ∈ C
I définie par :
(Ax|y)M = (x|A∗ y)N I N , ∀y ∈ C
∀x ∈ C IM
I M ×N et B = (bij ) ∈ C
Exercice 49 Montrer que si A = (aij ) ∈ C I N ×M ,
B = A∗ ⇔ ∀(i, j) bij = āji 1 ≤ i ≤ N, 1 ≤ j ≤ M
I M ×N alors on a
Exercice 50 Montrer que si A ∈ K
⊥ ⊥
KerA∗ = (ImA) et ImA∗ = (KerA)

Produit de matrices
Le produit matriciel correspond à la composition des applications linéaires.
I M ×P et B ∈ K
Pour A ∈ K I P ×N on définit la matrice produit
M ×N
C = A.B ∈ K
I par
P
X
cij = aik bkj , 1 ≤ i ≤ M, 1 ≤ j ≤ N
k=1

Remarque :
Le produit matriciel est associatif mais non commutatif.
t
I M ×P et B ∈ K
Exercice 51 Pour A ∈ K I P ×N , montrer que (A.B) = B t .At et
∗ ∗ ∗
(A.B) = B .A

109
Sous matrices
On appelle sous matrice d’une matrice A = (aij ) ∈ K I M ×N toute matrice de
la forme :  
ai1 j1 ai1 j2 . . . ai1 jq
 ai2 j1 ai2 j2 . . . ai2 jq 
 
 .. .. .. 
 . . . 
aip j1 aip j2 . . . aip jq
où les entiers ik et jl vérifient

1 ≤ i1 < i2 < . . . < ip ≤ M ; 1 ≤ j1 < j2 < . . . < jq ≤ N

Remarque :
Une sous matrice de A s’obtient par suppression de lignes ou de colonne de A.
Soient E et F deux espaces vectoriels de dimensions finies avec les décompositions
en sommes directes suivantes :

E = E1 ⊕ E2 ⊕ . . . ⊕ EN , F = F1 ⊕ F2 ⊕ . . . ⊕ FM

Soit A la matrice d’une application linéaire de E dans F . A ces décompositions

de E et F on associe une décomposition par blocs de A
 
 A11 A12 ... A1N 
 
 
 
 A21 A22 ... A2N 
A=  = (AIJ )
 
 
 .. .. .. 
 . . ... . 
 
 
AM 1 AM 2 ... AM N

Chaque AIJ est une sous matrice de type (mI , nJ ) représentant une application
linéaire de EJ dans FI
Soit A = (AIK ) et B = (BKJ ) deux matrices décomposées par blocs avec la
même décomposition en K. Alors on a le produit par bloc
X
A.B = (CIJ ) avec CIJ = AIK BKJ
K

Le produit par bloc s’étend au produit matrice vecteur qui est un cas parti-
culier du produit matriciel.

5.1.4 Matrices carrées

Lorsque A ∈ K I N ×N on dit que A est une matrice carrée ou une matrice
d’ordre N . Il faut observer que K I N ×N a une structure d’anneau unitaire dont
l’élément unité est I = (δij )

110
Matrice inversible
Une matrice A est dite inversible s’il existe une (et une seule) matrice notée
A−1 appelée inverse de A, telle que A.A−1 = A−1 .A = I. Dans le cas contraire
on dit que A est singulière.

I N ×N sont inversibles alors on a :

Exercice 52 Si A et B ∈ K
−1 −1 t −1 ∗
(A.B) = B −1 A−1 , (At ) = (A−1 ) , (A∗ ) = (A−1 )

Changement de base
I N Soit E 0 = (e01 , e02 , . . . , e0N ) une
Soit E = (e1 , e2 , . . . , eN ) une base de K
N
nouvelle base de K
I
Soit S la matrice carrée définie par

e0i = Sei 1≤i≤N

Comme S est de rang N alors S −1 existe et on a

ei = S −1 e0i 1≤i≤N

I N representé par x dans la base E et par x0 dans E 0 .

Soit un vecteur de K
Alors on a :
x = Sx0 ou x0 = S −1 x
S est dite matrice de passage de la base E à la base E 0
Soit f ∈ L(E, F ) représentée par la matrice A ∈ K I M ×N dans la base E de
E et la base F de F .
Soit S la matrice de passage de la base E à la base E 0 dans E.
Soit T la matrice de passage de la base F à la base F 0 dans F .
L’application linéaire f est alors représentée dans les bases E 0 et F 0 par

A0 = T −1 AS

La matrice A0 est dite équivalente à A

Lorsque E = F avec E = F et E 0 = F 0 alors on a T = S puis A0 = S −1 AS
est dite semblable à A

Matrices spéciales
I N ×N est dite :
Une matrice A ∈ K

symétrique si A est réelle et A = At

anti-symétrique si A est réelle et A = −At
hermitienne si A = A∗
anti-hermitienne si A = −A∗
orthogonale si A est réelle et AAt = At A = I

111
unitaire si AA∗ = A∗ A = I
normale si AA∗ = A∗ A
définie positive si (Ax|x) > 0 I N
0 6= x ∈ K
P
à diagonale dominante si |aii | > j6=i |aij | ∀i
permutation si A = (eσ(1) , . . . , eσ(N ) ) où σ une permutation
de (1, 2, . . . , N )
diagonale si aij = 0 lorsque i 6= j
tridiagonale si aij = 0 lorsque |i − j| > 1
bidiagonale supérieure si aij = 0 lorsque i > jou j > i + 1
triangulaire supérieure si aij = 0 lorsque i > j
hessenberg supérieure si aij = 0 lorsque i > j + 1

Déterminant d’une matrice

I N ×N le scalaire det(A)
On appelle déterminant d’une matrice carrée A ∈ K
défini par la relation

det(A) = a si I 1×1
A = (a) ∈ K
N
X j+1
det(A) = (−1) a1j det(A1j )
j=1

où A1j ∈ KI N −1×N −1 est la sous matrice de A obtenue par suppression de la

ligne 1 et la colonne j.

Exercice 53 Pour les matrices carrées d’ordre N montrer qu’on a

det(I) =1
det(At ) = det(A)
det(A∗ ) = det(A)
det(αA) = αN det(A) où α ∈ KI
det(AB) = det(A)det(B)
det(A−1 ) = det(A)
1
si A−1 existe

Remarque :
On montre que σ désignant une permutation de (1, 2, . . . , N ) et εσ sa signature,
on a
X Y N
detA) = εσ aσ(i)i
σ i=1

Valeurs propres
I N ×N .
Soit A ∈ K

112
On appelle valeur propre de A, le nombre λ ∈ C IN
I tel qu’il existe 0 6= x ∈ C
vérifiant Ax = λx. le vecteur x est alors appelé vecteur propre associé à λ.
On appelle élément propre de A un couple (λ, x) ∈ C I ×C I N où λ est valeur
propre et x vecteur propre associé.
On appelle spectre de A l’ensemble sp(A) de ses valeurs propres. Le rayon
spectral de A est le réel
ρ(A) = max |λ|
λ∈sp(A)

Remarque :
Si (λ, x) est un élément propre de A alors 0 6= x ∈ Ker(A − λI). Ainsi la
matrice A − λI est singulière. En particulier si 0 est valeur propre de A alors A
est singulière.
On appelle polynôme caractéristique d’une matrice carrée A l’application

pA : λ ∈ C
I 7→ pA (λ) = det(A − λI)

Les N racines complexes λi , 1 ≤ i ≤ N du polynôme caractéristique sont les

valeurs propres de A
On appelle trace de A = (aij ) le scalaire
N
X
tr(A) = aii
i=1

Remarque :
On montre que pour sp(A) = {λi , 1 ≤ i ≤ N } on a :
N
X N
Y
tr(A) = λi et det(A) = λi
i=1 i=1

Remarque :
I M ×N . On appelle valeur singulière de A la racine carrée d’une valeur
Soit A ∈ K
propre de la matrice hermitienne A∗ A

I M ×N . Montrer que A∗ A possède N valeurs propres

Exercice 54 Soit A ∈ K
positives.

Exercice 55 (Théorème de Gerschgorin) Soit A une matrice carrée d’or-

dre N . Soit Dk le disque de Gerschgorin de centre akk défini par
 
 XN 
Dk = z ∈ C I , |z − akk | ≤ |akj |
 
j=1, j6=k

Montrer que
N
[
λ ∈ sp(A) ⇒ λ ∈ Dk
k=1

113
5.1.5 Réduction de matrices
Théorème 5.1.1 (1) Etant donné une matrice carrée A, il existe une matrice
unitaire U telle que la matrice U −1 AU soit triangulaire.
(2) Etant donné une matrice normale A, il existe une matrice unitaire U
telle que la matrice U −1 AU soit diagonale.
(3) Etant donné une matrice symétrique A, il existe une matrice orthogonale
O telle que la matrice O−1 AO soit diagonale.

Preuve
La démontration est laissée en exercice.

Théorème 5.1.2 Si A une matrice carrée réelle, il existe deux matrices or-
thogonales U et V telles que

U t AV = diag(µi )

et si A une matrice carrée complexe, il existe deux matrices unitaires U et V

telles que
U ∗ AV = diag(µi )
où les nombres µi sont les valeurs singulières de A

Preuve
La démontration est laissée en exercice.

5.2 Normes vectorielles et matricielles

5.2.1 Normes vectorielles
Definition
Soit E un espace vectoriel sur K I . On appelle norme sur E toute application
kk : E → RI vérifiant les propriétés suivantes :

1) kxk = 0 ⇔ x = 0 et kxk ≥ 0 ∀x ∈ E
2) kαxk = |α| kxk ∀α ∈ KI , ∀x ∈ E
3) kx + yk ≤ kxk + kyk ∀x, y ∈ E

Exemple
I N nous utiliserons couramment les normes suivantes :
Sur E = K
PN
1) x 7→ kxk1 = q i=1 |xi |
PN 2
2) x 7→ kxk2 = i=1 |xi |
3) x 7→ kxk∞ = max1≤i≤N |xi |

114
I N . Montrer que pour tout réel p ≥ 1 l’application k.kp
Exercice 56 Soit E = K
définie par :
N
! p1
X p
kxkp = |xi |
i=1

est une norme sur E

Pour ce faire on utilisera la convexité de la fonction exponentielle pour mon-
trer le lemme de Young :

Lemme 5.2.1 Soit p, q > 1 vérifiant , p1 + 1q = 1.

Soit α, β ∈ R
I + , on a :
αp βq
αβ ≤ +
p q
Puis on en déduira l’inégalité de Holder :
N
X
|xi yi | ≤ kxkp kykq
i=1

I N.
Considerons une suite (x(n) ) d’éléments de E = K
On dit que la suite tend vers x ∈ E si

lim kx(n) − xk = 0
n→∞

Propriétés
On dit que deux normes k.k et k.k∗ sont équivalentes s’ il existe deux con-
stantes C et C 0 telles que

kxk ≤ Ckxk∗ et kxk∗ ≤ C 0 kxk

Théorème 5.2.1 Dans un espace vectoriel de dimension finie toutes les normes
sont équivalentes

5.2.2 Normes matricielles

Definition
I M ×N . Considérons deux normes vectorielles sur K
Soit une matrice A ∈ K I N
M
et K
I
On appelle norme matricielle induite par les normes vectorielles, le nombre :

kAk = max kAxk

kxk=1

Exemple
Si A = (aij )1≤i,j≤n alors on montre que :

115
1. La norme `∞
n
X
kAk∞ = max kAxk∞ = max |aij |
kxk∞ =1 1≤i≤n
j=1

2. La norme `1
n
X
kAk1 = max kAxk1 = max |aij |
kxk1 =1 1≤j≤n
i=1
3. La norme `2 p
kAk2 = max kAxk2 = ρ(AA∗ )
kxk2 =1

Remarque :
I N ×N on a :
Pour toute norme matricielle induite, pour la matrice unité I ∈ K
kIk = 1
Remarque :
Toute norme matricielle induite est une norme
Remarque :
Il existe des normes matricielles qui ne sont induites par aucune norme vecto-
rielle.
Remarque :
Soit A ∈ KI N ×N , on a ρ(A) ≤ kAk pour toute norme matricielle.
Exercice 57 Norme de Schur
I M ×N on appelle norme de Schur de A, le réel
Soit A ∈ K
  21
XM X N
2
kAk = S |aij | 
i=1 j=1

Montrer que k.kS est une norme sur KI M ×N √

N ×N
Montrer que pour I ∈ K
I on a : kIkS = N

Propriétés
Lemme 5.2.2 Pour toute norme matricielle induite, on a
kAxk ≤ kAkkxk I N
∀x ∈ K
On dit qu’elle est compatible.
Preuve
La démontration est laissée en exercice.
Proposition 5.2.1 Pour toute norme matricielle induite, pour deux matrices
A∈KI M ×N et B ∈ K
I N ×P on a
kABk ≤ kAkkBk et ρ(A) ≤ kAk
Preuve
La démontration est laissée en exercice.

116
5.2.3 Suites vectorielles, suites matricielles
On rappelle qu’une suite de matrices A(k) ∈ K I M ×N est assimilable à une
MN
suite de vecteurs de K
I . Donc sa convergence équivaut à la convergence de
M N suites scalaires et elle est indépendante de la norme choisie.
Exemple
Etude de la convergence de la suite (xn )n définie par
p
−n 1 2
xn = e cos n , n sin , n +n−n
n

On observe que pour tout n ≥ 1 on a

e−n cos n ≤ e−n et lim e−n = 0

n→∞

donc lim e−n cos n = 0

n→∞
Puis on a

1 1 1 1 1 1
lim n sin = lim n − = lim 1 − =1
n→∞ n n→∞ n 6 n3 n→∞ 6 n2

Et enfin on a
p n 1
lim n2 + n − n = lim √ =
n→∞ n→∞ n2 +n+n 2

1
Ainsi lim xn = 0, 1,
n→∞ 2

I 3 par
Exercice 58 Etudier de la convergence de la suite (xn )n définie dans R

n2 + 1 1 + 3 + 5 + . . . + (2n − 1)

1
xn = e ,n ,
1 − n2 n2

Critère de convergence
Soit A une matrice carrée d’ordre N
On dit que la matrice A est convergente si on a limk→∞ Ak = 0 c’est à dire
(k)
si limk→∞ aij = 0 pour tout i, j
Pour la convergence de la suite de puissance de matrice, on a le critère :

Théorème 5.2.2 Soit A une matrice carrée. On a les conditions équivalentes

suivantes :
(1) limk→∞ Ak = 0,
(2) I N
limk→∞ Ak x = 0, ∀x ∈ K
(3) ρ(A) < 1
(4) kAk < 1 pour au moins une norme induite

117
Preuve
La démontration est laissée en exercice.

Théorème 5.2.3 Soit A une matrice carrée. Soit k.k une norme matricielle
(induite) quelconque. Alors on a
1
lim kAk k k = ρ(A)
k→∞

Preuve
La démontration est laissée en exercice.

Théorème 5.2.4 Soit A une matrice carrée. La série I + A + A2 + . . . converge

−1
vers (I − A) si et seulement si ρ(A) < 1

Preuve
La démontration est laissée en exercice.
Exemple
Etude de la convergence de la matrice A et déduire éventuellement celle de la
Xn
suite (Bn )n où Bn = Ai , pour
i=0

− 12

0
A= 5 1
4 3

A cet effet, cherchons le rayon spectral de A.

1 1 1 1
detA − λI = 0 ⇔ − − λ −λ =0 ⇔ λ∈ − ,
2 3 2 3
1
D’où ρ(A) = . D’après le théorème 5.2.2 la matrice A est convergente. Par
2
n
−1
X
conséquent Bn = Ai → B = (I − A) que nous allons évaluer.
i=0

3 2

0 0
I −A= 2
5 2 ⇒ B= 3
5 3
−4 3 4 2

Exercice 59 Etude de la convergence de A et déduire éventuellement celle de

n
X
la suite (Bn )n où Bn = Ai pour :
i=0
1.
1

0 2
A= 1
2 0
2.
1 0
A= 1 1
4 2

118
Chapitre 6

Conditionnement

6.1 Conditionnement d’une matrice

Définition 6.1.1 Soit k.k une norme matricielle. Soit A ∈ K I N ×N une matrice
régulière. On appelle conditionnement de la matrice A associé à cette norme, le
réel
cond(A) = kAkkA−1 k
Exemple
On considère la matrice :
 
2 −1 0
A =  −1 2 −1 
0 −1 2
dont l’inverse est
3 1 1
 
4 2 4
A−1 =  1
2 1 1
2

1 1 3
4 2 4
Evaluons Cond∞ (A)
Si A = (aij ) alors
3
X
|a1j | = |2| + |−1| + |0| = 3
j=1
X3
|a2j | = |−1| + |2| + |−1| = 4
j=1
X3
|a3j | = |0| + |−1| + |2| = 3
j=1

et
3
X
kAk∞ = max |aij | = 4
1≤i≤3
j=1

119
Si A−1 = (aij ) alors
3
X 3 1 1 3
|a1j | = + + =
j=1
4 2 4 2
3
X 1 1
|a2j | = + |1| + =2
j=1
2 2
3
X 1 1 3 3
|a3j | = + + =
j=1
4 2 4 2

et
3
X
kA−1 k∞ = max |aij | = 2
1≤i≤3
j=1

Ainsi
Cond∞ (A) = kAk∞ kA−1 k∞ = 4 × 2 = 8

Exercice 60 On considére la matrice

 
4 −1 7
A =  −1 4 0 
−7 0 4

Sachant que
 
0.0625 0.015625 −0.109375
A−1 =  0.015625 0.25390625 −0.02734375 
0.109375 0.02734375 0.05859375

on évalue Cond1 (A).

On rappelle que si (aij )
n
X
kAk1 = max |aij |
1≤j≤n
i=1

Proposition 6.1.1 1. cond(αA) = cond(A) pour toute matrice A et tout

scalaire α 6= 0.
2. cond(A) ≥ 1 pour toute norme matricielle induite
µmax
3. cond2 (A) = où µmax et µmin sont respectivement la plus grande
µmin
et la plus petite valeur singulière de A.
4. cond2 (A) = 1 si et seulement si A = αQ où α est un scalaire et Q une
matrice unitaire.

Preuve

120
1.
−1
cond(αA) = kαAkk(αA) k
= |α| kAk α−1 kA−1 k
= kAkkA−1 k = cond(A)
2. I = AA−1 Donc 1 = kIk = kAA−1 k ≤ kAkkA−1 k = cond(A)
3. On sait que
12
t
kAk2 = max λi (A A) = max µi = µmax
1≤i≤N

et
−1
2
−1 t
kA k2 = max λi (A A)
1≤i≤N
1 1
= max =
µi µmin
µmax
donc cond2 (A) = kAk2 kA−1 k2 =
µmin
Dans le cas où A est symétrique on a µi = |λi (A)|. Ainsi
max λi (A)
cond2 (A) =
min λi (A)

4. D’après 3) cond2 (A) = 1 si et seulement si toutes les valeurs singulières

de A sont égales entre elles.
Or d’après les théorèmes de réduction de matrice il existe deux matrices
unitaires U et V et une matrice diagonale Σ = diag(µi ) = αI telles que
A = U ΣV ∗ = αU V ∗ = αQ avec Q = U V ∗ qui est unitaire.
Remarque :
On dit qu’une matrice est bien conditionnée si son conditionnement est proche
de 1.
Les matrices unitaires sont les mieux conditionnées possibles.

6.2 Conditionnement d’un système linéaire

Dans cette section, on met en évidence l’influence que peut avoir le nombre
de conditionnement de la matrice d’un système linéaire sur la propagation de
l’erreur éventuelle sur les données du problème.

Proposition 6.2.1 Soit A une matrice inversible. Soient x et x + ∆x les solu-

tions des systèmes linéaires :

Ax = b et A(x + ∆x) = b + ∆b

on a
k∆xk k∆bk
≤ cond(A)
kxk kbk

121
Preuve
Par différence des deux équations Ax = b et A(x + ∆x) = b + ∆b on obtient
A(∆x) = ∆b soit ∆x = A−1 ∆b. Ainsi

k∆xk = kA−1 ∆bk ≤ kA−1 kk∆bk

Par ailleurs kbk = kAxk ≤ kAkkxk soit

1 kAk
≤
kxk kbk

d’où par produit on a

k∆xk kAk kAk k∆bk

≤ k∆xk ≤ kA−1 kk∆bk = cond(A)
kxk kbk kbk kbk

Proposition 6.2.2 Soit A une matrice inversible. Soient x et x + ∆x les solu-

tions des systèmes linéaires :

Ax = b et (A + ∆A)(x + ∆x) = b

on a
k∆xk k∆Ak
≤ cond(A)
kx + ∆xk kAk

Preuve
En effet on a Ax = b = (A + ∆A)(x + ∆x) donc

A∆x + ∆A(x + ∆x) = 0 c’est à dire ∆x = −A−1 ∆A(x + ∆x)

donc
k∆xk = kA−1 ∆A(x + ∆x)k ≤ kA−1 kk∆Akkx + ∆xk
d’où
k∆xk k∆Ak
≤ kA−1 kk∆Ak = cond(A)
kx + ∆xk kAk

122
Bibliographie

[1] J. D. Faires et R. Burden , Numerical Methods, Edition 2, ed. Brooks/Cole

Publishing Company (1998)
[2] J. D. Faires et R. Burden , Numerical Analysis, Edition 9, ed. Brooks/Cole
Cengage Learning (2010)
[3] W. Chenay et D. Kincay, Numerical Mathmatics and Computing, ed
Brooks/Cole Publishing company (2000)
[4] M. Crouzeix et A. L. Mignot , Analyse Numérique des équations
différentielles , Edition 2, ed. Masson (1989)
[5] P. Ciarlet , Introduction à l’analyse numérique matricielle et à l’optimisa-
tion, Edition 4, ed. Masson (1990)
[6] P. Lascaux et R. Theodor , Analyse numérique matricielle appliquée, Tome
1 et 2, ed. Masson (1986)
[7] S. L. Ross , Differential Equations, Edition 3, ed. John Wiley & Sons Inc.
(1984)
[8] P. Blanchard, R. L. Devaney et G. R. Hall , Differential Equations, ed.
Brooks/Cole Publishing Company (1997)

123

Vous aimerez peut-être aussi

Notes Dec Ours
Pas encore d'évaluation
Notes Dec Ours
95 pages
Polycope AN1 2021 2022
Pas encore d'évaluation
Polycope AN1 2021 2022
120 pages
Guide Complet sur l'Analyse Numérique
Pas encore d'évaluation
Guide Complet sur l'Analyse Numérique
168 pages
An Numerique
Pas encore d'évaluation
An Numerique
168 pages
Un Cours D'Analyse Numérique
Pas encore d'évaluation
Un Cours D'Analyse Numérique
168 pages
Introduction à l'Analyse Numérique
100% (1)
Introduction à l'Analyse Numérique
156 pages
Exercices Résolus en Analyse Numérique
Pas encore d'évaluation
Exercices Résolus en Analyse Numérique
5 pages
Lakrib AN PDF
0% (1)
Lakrib AN PDF
5 pages
Méthodes Numériques en Mathématiques
Pas encore d'évaluation
Méthodes Numériques en Mathématiques
46 pages
Cours de Méthodes Numériques
Pas encore d'évaluation
Cours de Méthodes Numériques
80 pages
Poly PDF
Pas encore d'évaluation
Poly PDF
168 pages
Calcul Numérique et Analyse d'Erreur
Pas encore d'évaluation
Calcul Numérique et Analyse d'Erreur
133 pages
Exercices Résolus en Analyse Numérique
Pas encore d'évaluation
Exercices Résolus en Analyse Numérique
6 pages
Livre Analyse Numerique Chap6
Pas encore d'évaluation
Livre Analyse Numerique Chap6
11 pages
Exercices Résolus en Analyse Numérique
Pas encore d'évaluation
Exercices Résolus en Analyse Numérique
6 pages
Cours Analnum
Pas encore d'évaluation
Cours Analnum
42 pages
Cours AnalNum Agreg FBoyer$$$$$$
Pas encore d'évaluation
Cours AnalNum Agreg FBoyer$$$$$$
129 pages
Cours d'Analyse Numérique 2019/2020
Pas encore d'évaluation
Cours d'Analyse Numérique 2019/2020
84 pages
Cours Anal Num31-10-2018
Pas encore d'évaluation
Cours Anal Num31-10-2018
37 pages
Notes de cours sur les méthodes numériques
Pas encore d'évaluation
Notes de cours sur les méthodes numériques
94 pages
Erreurs en Analyse Numérique
Pas encore d'évaluation
Erreurs en Analyse Numérique
15 pages
polyAnaNum PDF
Pas encore d'évaluation
polyAnaNum PDF
176 pages
Notes de Cours
Pas encore d'évaluation
Notes de Cours
47 pages
Arithmétique des ordinateurs et erreurs
Pas encore d'évaluation
Arithmétique des ordinateurs et erreurs
170 pages
Introduction à l'Analyse Numérique
Pas encore d'évaluation
Introduction à l'Analyse Numérique
45 pages
COURS MethodesNumeriques2015 - SOUK - AHRAS PDF
0% (1)
COURS MethodesNumeriques2015 - SOUK - AHRAS PDF
125 pages
TD Et Cours
Pas encore d'évaluation
TD Et Cours
125 pages
DownloadAnalyse Numérique
100% (1)
DownloadAnalyse Numérique
137 pages
Support Cours PHY200
Pas encore d'évaluation
Support Cours PHY200
34 pages
Travaux Pratiques D'analyse Numérique Avec Python
Pas encore d'évaluation
Travaux Pratiques D'analyse Numérique Avec Python
25 pages
Cours d'Analyse Numérique - Université Burundi
Pas encore d'évaluation
Cours d'Analyse Numérique - Université Burundi
71 pages
PolyMethNumBentobache2020 2021
Pas encore d'évaluation
PolyMethNumBentobache2020 2021
50 pages
Integration Et Derivation Numeriques
Pas encore d'évaluation
Integration Et Derivation Numeriques
26 pages
Méthodes Numériques et Systèmes Linéaires
Pas encore d'évaluation
Méthodes Numériques et Systèmes Linéaires
73 pages
Introduction à l'analyse numérique
Pas encore d'évaluation
Introduction à l'analyse numérique
63 pages
Ananum
Pas encore d'évaluation
Ananum
49 pages
Analyse Numérique
Pas encore d'évaluation
Analyse Numérique
51 pages
Erreurs et stabilité en analyse numérique
Pas encore d'évaluation
Erreurs et stabilité en analyse numérique
117 pages
Cours MN
Pas encore d'évaluation
Cours MN
47 pages
Analyse Numérique : Méthodes et Applications
Pas encore d'évaluation
Analyse Numérique : Méthodes et Applications
97 pages
Syllabus Analyse Numerique 2020
Pas encore d'évaluation
Syllabus Analyse Numerique 2020
209 pages
Analyse numérique : cours et méthodes
Pas encore d'évaluation
Analyse numérique : cours et méthodes
261 pages
Introduction à l'Analyse Numérique
Pas encore d'évaluation
Introduction à l'Analyse Numérique
43 pages
Analyse Numerique
Pas encore d'évaluation
Analyse Numerique
59 pages
Agrégation Externe en Analyse Numérique
Pas encore d'évaluation
Agrégation Externe en Analyse Numérique
115 pages
Résolution Numérique de Systèmes Linéaires
Pas encore d'évaluation
Résolution Numérique de Systèmes Linéaires
72 pages
SUPPORT DE COURS Introduction Au Calcul Scientifique
Pas encore d'évaluation
SUPPORT DE COURS Introduction Au Calcul Scientifique
52 pages
Cours Analyse Numérique
Pas encore d'évaluation
Cours Analyse Numérique
57 pages
Méthodes Numériques pour Systèmes Linéaires
Pas encore d'évaluation
Méthodes Numériques pour Systèmes Linéaires
68 pages
Cours MN
Pas encore d'évaluation
Cours MN
47 pages
An Non Lineaire 2
Pas encore d'évaluation
An Non Lineaire 2
86 pages
Equations Diff
Pas encore d'évaluation
Equations Diff
4 pages
Concours Commun Mines-Ponts (CCMP) 1989 MP-PC-TSI Mathématiques Pratique Ea
Pas encore d'évaluation
Concours Commun Mines-Ponts (CCMP) 1989 MP-PC-TSI Mathématiques Pratique Ea
2 pages
1-Equations Diff
Pas encore d'évaluation
1-Equations Diff
10 pages
BTS 2024 Mathématiques : Épreuve et Exercices
Pas encore d'évaluation
BTS 2024 Mathématiques : Épreuve et Exercices
8 pages
Équations Différentielles en Calculus 2
Pas encore d'évaluation
Équations Différentielles en Calculus 2
51 pages
Résolution d'EDO par Runge-Kutta
Pas encore d'évaluation
Résolution d'EDO par Runge-Kutta
5 pages
Exercices sur les équations diff.
Pas encore d'évaluation
Exercices sur les équations diff.
3 pages
Utilisation Des Équations Différentielles en Finance
Pas encore d'évaluation
Utilisation Des Équations Différentielles en Finance
4 pages
Optimisation et contrôle en mathématiques
Pas encore d'évaluation
Optimisation et contrôle en mathématiques
213 pages
Analyse - 2 - Les Équations Différentielles
Pas encore d'évaluation
Analyse - 2 - Les Équations Différentielles
12 pages
TD 7
Pas encore d'évaluation
TD 7
2 pages
Devoir 11 Espaces Vectoriels Applications Lineaires
Pas encore d'évaluation
Devoir 11 Espaces Vectoriels Applications Lineaires
2 pages
Sujet Bts Ol Mathematiques
Pas encore d'évaluation
Sujet Bts Ol Mathematiques
8 pages
Équations Différentielles Terminale
Pas encore d'évaluation
Équations Différentielles Terminale
11 pages
Résolution numérique des EDO en 3ème année
Pas encore d'évaluation
Résolution numérique des EDO en 3ème année
19 pages
TD3 Sce
Pas encore d'évaluation
TD3 Sce
6 pages
TP Matlab
100% (2)
TP Matlab
12 pages
Dérivabilité 1bac
Pas encore d'évaluation
Dérivabilité 1bac
7 pages
EDP Avec Laplace
Pas encore d'évaluation
EDP Avec Laplace
8 pages
Équations Différentielles : Exercices et Solutions
Pas encore d'évaluation
Équations Différentielles : Exercices et Solutions
2 pages
Cours
Pas encore d'évaluation
Cours
10 pages
Électronique de puissance et convertisseurs
Pas encore d'évaluation
Électronique de puissance et convertisseurs
160 pages
Équations Différentielles : Théorie et Pratique
Pas encore d'évaluation
Équations Différentielles : Théorie et Pratique
9 pages
Équations Différentielles Linéaires
Pas encore d'évaluation
Équations Différentielles Linéaires
24 pages
22 Équations Différentielles Ordinaires EDO
Pas encore d'évaluation
22 Équations Différentielles Ordinaires EDO
36 pages
Révision: Équations Différentielles
Pas encore d'évaluation
Révision: Équations Différentielles
3 pages
Cours MATLAB pour Ingénieurs
Pas encore d'évaluation
Cours MATLAB pour Ingénieurs
163 pages
COURS EQAUDIIF Niveau 3PDF
Pas encore d'évaluation
COURS EQAUDIIF Niveau 3PDF
39 pages
Equations Differentielles 3
Pas encore d'évaluation
Equations Differentielles 3
1 page
Cour Régulation Numérique.
Pas encore d'évaluation
Cour Régulation Numérique.
26 pages