Anlyse M3
Anlyse M3
Anlyse M3
60
50
40
30
20
10
−10
3
2.5 60
50
2 40
30
1.5 20
10
1 0
Abdesslam B OUTAYEB
Abdelaziz C HETOUANI
Mohamed D EROUICH
Table des matières
1 Rappels et définitions 6
1.1 Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2 Quelques rappels sur les matrices . . . . . . . . . . . . . . . . . . . . 6
1.3 Valeurs et vecteurs propres . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Normes vectorielles et normes matricielles . . . . . . . . . . . . . . . 11
1.4.1 Normes vectorielles . . . . . . . . . . . . . . . . . . . . . . . . 11
1.4.2 Normes matricielles . . . . . . . . . . . . . . . . . . . . . . . . 12
1.4.3 Normes compatibles . . . . . . . . . . . . . . . . . . . . . . . . 12
1.5 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.6 Complément bibliographique sur les matrices positives . . . . . . . . 21
1.7 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1
2.6 Factorisation de Householder (matrice unitaire ) . . . . . . . . . . . . 35
2.7 Conditionnement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.8 Matrices creuses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.9 Résultats sur les matrices non carrées . . . . . . . . . . . . . . . . . . 43
2.10 Résolution des systèmes à matrices non carrées . . . . . . . . . . . . 44
2.11 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.12 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2
4.2.1 Schéma de Horner . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.2.2 Méthode de Laguerre . . . . . . . . . . . . . . . . . . . . . . . 92
4.2.3 Méthode de Bairstow . . . . . . . . . . . . . . . . . . . . . . . . 94
4.2.4 Méthode de Maehly . . . . . . . . . . . . . . . . . . . . . . . . 95
4.2.5 Localisation des racines . . . . . . . . . . . . . . . . . . . . . . 95
4.2.6 Nombre de racines réelles d’un polynôme . . . . . . . . . . . 98
4.2.7 Cas des racines isolées . . . . . . . . . . . . . . . . . . . . . . . 99
4.3 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
4.4 Complément bibliographique . . . . . . . . . . . . . . . . . . . . . . . 107
4.5 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
3
7.5.2 Consistance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
7.5.3 Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
7.5.4 Méthode d’Euler . . . . . . . . . . . . . . . . . . . . . . . . . . 162
7.5.5 Méthodes de Taylor dans le cas scalaire . . . . . . . . . . . . . 164
7.5.6 Méthodes de Runge-Kutta (R.K) dans le cas scalaire . . . . . . 165
7.5.7 Méthodes de Runge-Kutta explicites . . . . . . . . . . . . . . . 165
7.6 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
7.7 Complément bibliographique . . . . . . . . . . . . . . . . . . . . . . . 182
7.8 Exercices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184
Annexe 190
Bibliographie 199
4
Liste des figures
5
Chapitre 1
Rappels et définitions
1.1 Notations
A matrice carrée d’ordre n.
ai j élément de la ieme ligne et la jeme colonne.
A−1 inverse de A.
A> transposée de A.
det( A) déterminant de A.
trace( A) trace de A.
ρ( A) rayon spectral de A.
I matrice idendité d’ordre n.
0 matrice nulle d’ordre n.
x vecteur colonne d’éléments xi , i = 1, 2, · · · , n.
x> vecteur ligne d’éléments x j , j = 1, 2, · · · , n.
k Ak norme de A.
k xk norme de x.
h, i produit scalaire.
h x, yi = ∑ xi yi = y H x.
i
6
• A est dite symétrique si A est réelle et A> = A.
• Si A est une matrice réelle avec ai j ≤ 0 pour tout i 6= j, A est dite M-matrice
si A est non singulière et A−1 ≥ 0.
• Si A est une matrice réelle avec ai j ≤ 0 pour tout i 6= j, A est dite matrice de
Stieljes si A est symétrique et définie positive.
7
• A est dite à diagonale strictement dominante en lignes si elle vérifie la pro-
priété
n
∑ |ai j | < |aii | ∀i = 1, · · · , n.
j=1
j 6 =i
pour tout i = 1, · · · , n.
On a en particulier
Définition 1.3.2. Un vecteur x 6= 0 est dit vecteur propre de A associé à une valeur
propre λ si Ax = λ x.
Définition 1.3.4. Deux matrices A et B sont dites équivalentes s’il existe deux ma-
trices P et Q telles que A = PBQ.
A et B sont dites semblables s’il existe une matrice inversible P telle que A = PBP−1 .
8
Définition 1.3.5. Une matrice A est diagonalisable si elle est semblable à une ma-
trice diagonale.
Définition 1.3.6. Une matrice A est dite définie positive si elle vérifie
Remarque 1.3.1.
Pour toute matrice A, A∗ A est une matrice hermitienne définie positive. Il s’ensuit
que les valeurs propres de A∗ A sont réelles positives. Si ces valeurs propres sont
√ √ √
notées µ1 ≥ µ2 ≥ · · · ≥ µn ≥ 0 les racines carrées positives µ1 , µ2 , · · · , µn
sont appelées valeurs singulières de A.
Preuve:
Soit (e1 , e2 , · · · , en ) la base canonique de Cn .
On procède par récurrence sur n. Pour n = 1 le théorème est trivial.
On suppose que le résultat est vrai pour les matrices d’ordre n − 1 et soient A
une matrice carrée d’ordre n et λ1 une valeur propre de A associé à un vecteur
propre x1 6= 0 avec k x1 k2 = 1. On peut trouver n − 1 vecteurs x2 , · · · , xn tels que
( x1 , x2 , · · · , xn ) forment une base de Cn et la matrice X = [ x1 , · · · , xn ] soit unitaire,
X H X = I. Comme
X H AXe1 = X H Ax1 = λ1 X H x1 = λ1 x1 ,
9
existe une matrice unitaire U1 d’ordre n − 1 telle que
λ2 ∗ · · · ∗
.. ..
λ . .
3
U1H A1 U1 = .
..
. ∗
λn
La matrice à !
1 0
U=X ,
0 U1
est une matrice unitaire d’ordre n qui satisfait
à ! à !
1 0 1 0
U H AU = X H AX
0 U1H 0 U1
à !à !à !
1 0 λ1 a 1 0
=
0 U1H 0 A1 0 U1
λ1 ∗ · · · ∗
.. ..
λ . .
2
= .. .
. ∗
λn
Ce qui termine la preuve.
Si la matrice A est hermitienne, le théorème de Schur assure que
Corollaire 1.3.1.
Pour toute matrice hermitienne A, il existe une matrice unitaire U = [ x1 , · · · , xn ] telle
que:
λ1
λ2
U AU = U AU =
−1 H
..
,
.
λn
les valeurs propres de A sont réelles et A est diagonalisable.
Théorème 1.3.2. Si les valeurs propres λi , i = 1, · · · , n d’ une matrice hermitienne sont
rangées dans un ordre décroissant
λ1 ≥ λ2 ≥ · · · ≥ λn ,
10
Preuve:
Si U H AU = D = diag(λ1 , · · · , λn ), U unitaire, alors pour tout x 6= 0 on a
x H Ax ( x H U )U H AU (U H x) y H Dy
∑ λi | ηi |2
i
= = = ≤ λ1 ,
xH x ( x H U )(U H x) yH y ∑ | ηi |2
i
h Ax, xi
µR = , pour tout x 6= 0.
h x, xi
Remarque 1.3.2. Si A est une matrice hermitienne dont les valeurs propres sont
rangées par ordre croissant µ1 ≤ µ2 ≤ · · · ≤ µn , alors d’après le théorème 1.3.2 on
a
µ1 ≤ µ R ≤ µn .
Preuve:
n
λ étant valeur propre associée à un vecteur propre x, on a ∑ ai j x j − λ xi = − aii xi
j=1
j 6 =i
Définition 1.4.1. Une norme définie sur un espace vectoriel E est une application
notée k.k de E dans R+ telle que
i) k xk = 0 si et seulement si x = 0.
11
ii) kλ xk = |λ |k xk ∀λ ∈ R.
iii) k x + yk ≤ k xk + k yk ∀ x ∈ E ∀ y ∈ E.
Exemple 1.4.1.
à !1/2 à !1/ p
2 p
k x k 1 = ∑ | xi | , k x k 2 = ∑ | xi | , k xk p = ∑ | xi | , k xk∞ = max | xi |.
i i i i
Définition 1.4.2. Une norme matricielle est une norme qui vérifie, en plus de i),ii)
et iii) la condition suivante
Les normes vectorielles k.k1 , k.k∞ et k.k2 induisent, respectivement, les normes
p
matricielles suivantes k Ak1 = max ∑ | ai j |, k Ak∞ = max ∑ | ai j |, k Ak2 = ρ( A> A).
j i i j
Remarque 1.4.1. Il existe des normes matricielles qui ne sont subordonnées à au-
cune norme vectorielle.
µ ¶1/2
Exemple 1.4.2. Norme de Schur k Ak = ∑ ∑ ai2j .
i j
Définition 1.4.4. Une norme vectorielle k.k et une norme matricielle k.k∗ sont com-
patibles si pour toute matrice A et tout vecteur x; k Axk ≤ k Ak∗ k xk.
Théorème 1.4.1. Soit k.k une norme matricielle et B une matrice complexe vérifiant
k Bk < 1. Alors la matrice I ± B est inversible et on a:
i) ( I − B)−1 = I + B + B2 + · · ·
12
° ° 1
ii) °( I ± B)−1 ° ≤ .
1 − k Bk
Inversement, si la série de i) converge alors ρ( B) < 1.
Preuve:
Si I ± B n’était pas inversible alors elle admettrait λ = 0 pour valeur propre et par
suite on aurait k Bk ≥ 1. Soit G = ( I ± B)−1 il s’ensuit que ( I ± B) G = I c.à.d
1
k G k ≤ 1 + k Bkk G k soit k G k ≤ .
1 − k Bk
En écrivant S p = I + B + B2 + · · · + B p , on obtient
S p − BS p = ( I − B) S p = I − B p+1 .
Remarque 1.4.3. Sous MATLAB les fonctions matricielles les plus courantes sont:
det(A), eig(A), poly(A), inv(A), rank(A), norm(A), norm(A,1), norm(A,2), norm(A,inf),
norm(A,’fro’), cond(A,k).
Exemple 1.4.3. µ ¶
1
Considérons la matrice de Hilbert d’ordre 6, H =
i+ j−1 i, j=1,··· ,6
rang( H ) = 6.
det( H ) = 5.3673e−18 .
k H k1 = 2.4500 si k = 1,
norm( H, k) = k H k2 = 1.6189 si k = 2,
k H k∞ = 2.4500 si k = ∞,
° −1 °
° ° 7
k H k1 ° H °1 = 2.907e si k = 1,
cond( H, k) = k H k2 ° H −1 °2 = 1.4951e7 si k = 2,
° °
k H k∞ ° H −1 °∞ = 2.907e7 si k = ∞,
ρ( H ) = 1.6189.
Pour les matrices positives on a les résultats suivants
13
Théorème 1.4.2 (Perron-Frobenius ).
Si A est positive, alors A admet une valeur propre λ0 ≥ 0 égale au rayon spectral de A, à
λ0 correspond un vecteur propre positif.
Si m = 1 ° °
° n(t)
°
°
>° |λ2 | t
° λt − uv ° ≤ c
1 λ1
λ1
est dit rapport d’amortissement.
|λ2 |
Remarque 1.4.4. Toute matrice stochastique admet 1 comme valeur propre associée
au vecteur propre I dont toutes les composantes sont égales à 1.
Théorème 1.4.6. Soit A une matrice associée à une chaı̂ne de Markov. Alors il existe un
n(k)
vecteur p > 0 avec k pk = 1, lim Ak = A et lim = p, si n(k + 1) = An(k).
k−→∞ k−→∞ k n ( k )k
14
1.5 Applications
Application 1.5.1 (Modèle de Markov). [41]
Considérons une stratégie de gestion des terres en supposant quatre types:
T1=Terres agricoles non irriguées, T2=Terres agricoles irriguées, T3=Terres de con-
struction , T4=Terres non exploitables. Supposons que nous adoptions une stratégie
à long terme qui consiste à contrôler chaque année le passage d’un type de terre
à un autre pour arriver à un certain équilibre. Soit M une matrice de transition
d’un type de terre à un autre, si n(t) désigne l’état du système à l’instant t alors
¡ ¢
n(t + 1) = Mn(t) où M = pi j i, j=1,··· ,4 , pi j est la probabilité de passage du type
de terre i au type j au bout d’une année.
A titre d’exemple, si cette stratégie est appliquée durant une quarantaine d’années
avec
0.75 0.1 0.15 0
0.1 0.75 0.15 0
M= 0
,
0.2 0.7 0.1
0 0.05 0.25 0.7
on arrive à l’équilibre suivant:
0.1468 0.3714 0.36 0.12
0.1468 0.3714 0.36 0.12
M 37
=
0.1468 0.3714 0.36 0.12
0.1468 0.3714 0.36 0.12
Interprétation: à l’équilibre le domaine des terres non irriguées sera formé de 59%
de son étendue initiale, les domaines des terres irriguées et de construction auront
augmenté de presque 50% chacun (1.49 et 1.44). Enfin, les terres non exploitées
auront diminué de la moitié (0.52%). Ces chiffres ayant été obtenus en prenant
pour chaque colonne, les sommes sur les quatre lignes.
• La première classe d’âge est formée par les nouveaux nés, sa taille est donnée
par:
m
n1 (t + 1) = ∑ bi ni ( t ) .
i =1
15
• Chaque groupe d’âge i + 1 est formé par les individus du groupe d’âge i qui
survivent
ni+1 (t + 1) = pi ni (t) ; i ≥ 1.
n1 (t + 1) = b2 n2 (t) + b3 n3 (t)
ni +1 ( t + 1 ) = pi ni ( t )
où
0 b2 b3 0 0 1000
p1 0 0 0 0 900
L=
0 p2 0 0 0 et n(0) =
800
0 0 p3 0 0 700
0 0 0 p4 0 600
avec b1 = 1.5 , b3 = 1, p1 = 0.98, p2 = 0.96, p3 = 0.93, p4 = 0.9, comme
ni ( k + 1 )
lim = λ1 , donc d’après le tableau 1.1 on a λ1 = 1.4548. D’autre
k−→+∞ ni ( k )
n(k + 1)
part lim = v1 , donc d’après le tableau 1.2 on a
k−→+∞ k n ( k )k
v1 = (0.387, 0.261, 0.172, 0.1102, 0.068)> .
16
période classe 1 classe 2 classe 3 classe 4 classe 5
1 2.15 1.0889 1.08 1.0629 10.5
2 1.0856 2.15 1.0889 1.08 1.029
3 1.7572 1.0856 2.15 1.0889 1.08
4 1.3297 1.7572 1.0856 2.15 1.088
··· ··· ··· ··· ··· ···
··· ··· ··· ··· ··· ···
16 1.4550 1.455 1.4548 1.4548 1.4557
17 1.4548 1.4548 1.455 1.455 1.4547
18 1.4549 1.4549 1.4548 1.4548 1.4548
19 1.4548 1.4548 1.4549 1.4549 1.455
20 1.4549 1.4549 1.4548 1.4548 1.4548
17
dépend de l’équilibre entre les trois paramètres pré-cités et le problème réside à
juste titre dans le déséquilibre qui conduit à long terme aux complications diabé-
tiques( cécité, amputations, insuffisance rénale,etc ...). Les études récentes ont montré
que l’éducation diabétique et la bonne maı̂trise des paramètres permet une réduction
significative des complications. Par le modèle donné dans le présent exemple, les
auteurs entendent illustrer l’importance du contrôle du passage entre la classe des
diabétiques sans complications, Dk (t) et celle des diabétiques avec complications,
Ck (t) en tenant en compte de la structure d’âge. Pour les statistiques concernant
l’incidence, la prévalence, le coût du diabète et d’autres données, nous renvoyons
aux références( [22, 30, 31, 145]).
On suppose que les évenements tels que naissance, mortalité et guérison, sont
enregistrés tous les 15 ans. On néglige l’émigration et l’immigration et on suppose
que le nombre de males est égal au nombre de femelles. On suppose qu’il n’y a pas
de complication congénitale. On obtient alors le modèle compartimental suivant:
@ ¡
@ ¡
pk (t)(1 − sk ) ¡@ qk (t)(1 − ak )
qk (t)¡
ak k (t)sk
p@
? ¡ ª R
@ ?
Ik+1 (-
t) Dk + 1 ( t ) Ik0 +1 (t)
Ck+1 (t) ¾
pour k = 0, 1, · · · , m − 1
avec
ε k ( ti )
Ik+1 (ti+1 ) = D ( t i ) = α k + 1 ( t i ) Dk ( t i ) ,
γk0 (ti ) k
ε0k (ti )
Ik0 +1 (ti+1 ) = C (ti ) = αk0 +1 (ti )Ck (ti ) et I0 (ti ) = α0 (ti )n0 (ti ).
γk00 (ti ) k
On suppose que bk (ti ) = 0, b0k (ti ) = 0 et b00k (ti ) = 0 pour k < A0 et k > A1 .
18
En Additionnant Ck+1 (ti+1 ) et Dk+1 (ti+1 ) On obtient le modèle simplifié avec nk+1 :
µ ¶ µ ¶
m
θ ( ti ) 1 θ 0 (ti ) 00 00
n0 (ti+1 ) = I0 (ti ) + ∑
− 1 b k ( ti ) n k ( ti ) + bk (ti ) Dk (ti ) + bk (ti )Ck (ti )
k=0
2 γ k ( ti ) 2
r 0 ( ti +1 ) = 0
µ ¶
¡ ¢
nk +1 ( ti +1 ) = pk (ti ) + qk (ti ) − pk (ti ) rk (ti ) nk (ti ) + αk+1 (ti ) Dk (ti ) + αk0 +1 (ti )Ck (ti )
n(ti+1 ) = Ln(ti )
où
f0 f1 ··· fm
ν0 0 ··· 0
L=
.. .. .. ..
. . . .
0 ··· νm−1 0
µ ¶ µ ¶
θ 1 θ0 0 00 0 ∗ θ 1 θ0
avec f 0 = α0 + − 1 b 0 + ( b 0 + ( b 0 − b 0 )r 0 ), f k = − 1 bk + (b0k +
2 γ0 2 2 γk 2
(b00k − b0k )r∗k ).
µ ¶>
ν k = αk +1 + pk + (αk0 +1 + qk − αk+1 − pk )r∗k , n ( ti ) = n 0 ( ti ) , · · · , n m ( ti ) .
19
pk (t)µk (1 − sk ) qk (t)µk0 (1 − ak )
On obtient
n(t + 1) = Ln(t)
où
f0 f1 f2 ··· fm
ν0 ν10 0 ··· 0
L=
.. .. .. .. ..
. . . . .
0 0 · · · 0 νm−1 νm 0
µ ¶
θ 1 θ0 ¡ 0 ¢
avec f 0 = µ0 + α0 + − 1 b0 + b0 + (b000 − b00 )r∗0 .
µ ¶ 2 γ0 2
θ 1 θ0 ¡ 0 ¢
fk = − 1 bk + bk + (b00k − b0k )r∗k nk (t) pour k = 1, 2...m − 1,
2 γk 2
νk = pkωk + αk+1 + (αk0 +1 + qkω0k − αk+1 − pkωk )r∗k pour k = 0, 1 · · · , m − 1,
νk0 = pk µk + (qk µk0 − pk µk )r∗k+1 pour k = 1, · · · , m.
20
1.6 Complément bibliographique sur les matrices positives
Les matrices positives n’apparaissent pas toujours dans les livres d’analyse numé-
rique. Cependant, leur théorie et leur utilisation en écologie, statistique et numérique
de phénomènes réels revêt une grande importance comme en atteste le nombre de
publications et de livres dévoués à ce sujet. Il est difficile de donner une liste ex-
haustive mais on peut retracer briévement le ”chemin” des matrices positives à
travers les auteurs suivants:
Frobenius (1908) et Perron (1912) ont marqué la théorie des matrices positives
par leurs théorèmes bien connus concernant les valeurs et vecteurs propres des
matrices strictement positives, positives et irréductibles ou simplement positives.
L’importance de leurs résultats peut-être également mesurée par le nombre d’auteurs
qui ont proposé des démonstrations différentes des mêmes résultats, des générali-
sations ou des applications dans différents domaines: Bernadelli (1941), Collatz
(1942), Leslie (1945, 1948), Brauer(1957), Householder (1958), Gantmakher (1937,
1959), Karlin (1959), Mewborn (1960), Marcus & Minc (1964), Caswell (1980), Cia-
rlet (1984), Horn & Johnson (1985), Minc (1988), Logofet (1993), Axelsson (1994),
Berman & Plemmons (1994), Varga (2000).
Signalons enfin qu’un certain nombre de résultats ont été établis sur le lien entre
matrices positives et théorie des graphes.
21
1.7 Exercices
Exercice 1.7.1. Soit A une matrice carrée réelle, montrer que
i) k Ak1 = max ∑ | ai j |.
i
Exercice 1.7.2. Montrer que toute matrice de Leslie L admet une seule valeur propre
positive λ1 associée au vecteur propre v1 donnée par
à !>
p1 p1 p2 · · · pm−1
v1 = 1, , · · · , .
λ1 λ1m−1
3. Soient D et P deux matrices non singulières avec P telle que PAP−1 = U (tri-
° °
angulaire supérieure), en posant B = D −1 P, montrer que k Ak∗ ≤ ° D −1 UD °∞
et en déduire que, étant donné une matrice carrée A et un réel ε > 0, il existe
une matrice diagonale D = diag(di ) avec d > 0 et telle que: k Ak∗ ≤ ρ( A) + ε
4. Conclure.
° ° 1
montrer que k Axk∞ ≥ δ k xk∞ , en déduire que A−1 existe et que ° A−1 °∞ ≤ .
δ
22
7. On considère la matrice A donnée par
3 1 0 0
1 2 1 0
A= 0 1
avec ε > 0.
2 1 +ε
0 0 1 +ε 3
23
Chapitre 2
Ax = b (2.0.1)
où A est une matrice carrée supposée inversible, b un vecteur second membre et x
¡ ¢
le vecteur des inconnues, A = ai j i, j=1,··· ,n , b = (b1 , · · · , bn )> , x = ( x1 , · · · , xn )> .
Théoriquement, le fait que A soit inversible entraı̂ne que le système (2.0.1) admet
une solution unique x = A−1 b.
Mais cette écriture suppose que l’on dispose de la matrice A−1 , or l’obtention de
¡ ¢
A−1 est équivalente à la résolution de n systèmes linéaires, A. A−1 j = e j =
(0, · · · , 1, 0, · · · , 0)> en plus de la multiplication x = A−1 b.
Une autre méthode consisterait à obtenir les xi à l’aide des formules de Cramer
det( Ai )
xi = où det( Ai ) désigne le déterminant de la matrice obtenue en remplaçant
det( A)
la ieme colonne de A par le vecteur b.
Le calcul de chaque déterminant nécessite n.n! multiplications et (n! − 1) additions.
Soit au total : (n + 1)!n multiplications, (n + 1)(n! − 1) additions et n divisions.
A titre d’exemple, on a besoin de 4319 opérations si n = 5 et environ 41000 opérations
pour n = 10. Comme les problèmes d’analyse numérique donnent lieu à des ma-
trices de grandes tailles (n assez grand), la méthode de Cramer et les méthodes
similaires s’avèrent inutilisables.
24
2.1 Résolution d’un système par les méthodes de descente
ou de remontée
Ux = b (2.1.1)
ou
Lx = b (2.1.2)
u11 x1 + u12 x2 + · · · + u1n xn = b1
0 + u22 x2 + · · · + u2n xn = b2
Ux = b ⇔ .. ..
0 . .
unn xn = bn
En supposant que les ukk sont non nuls, on obtient les xi de façon évidente en
commençant
à ! et en remontant. On obtient ainsi xn = bn /unn puis
par le bas
n
xi = bi − ∑ ui j xi /uii , pour i = n − 1 à 1.
j =i + 1
L’algorithme de résolution est le suivant :
xn = bn /unn
Pour i = n − 1 à 1
xi = bi
Pour j = i + 1 à n
xi = xi − ui j ∗ x j
fin j
xi = xi /uii
fin i
n(n − 1) n(n − 1)
Le nombre d’opérations nécessaire est: multiplications addi-
2 2
2
tions et n divisions. Soit au total n opérations.
Remarque 2.1.1. Le cas d’un système avec matrice triangulaire inférieure se traite
de façon similaire en obtenant x1 d’abord puis en descendant.
25
2.2 Matrices élémentaires
k
2.2.1 Matrices élémentaires de Gauss
1
Soient
les matrices ..
1 0 . 0
k
−m21 1 0 0 1
M1 = .. , Mk =
.. .
. −mk+1k . .
. 0 ..
.
.. .. ..
−mn1 1 . . .
0 −mk+1k 1
en posant ek = (0, · · · , 1, 0, · · · , 0)> et mk = (0, · · · , 0, −mk+1k , · · · , −mnk )> ,
on obtient Mk = I − mk e> k et on vérifie facilement que Mk est inversible et que
−1 >
Mk = I + mk ek .
Pk = In − 2ωkω>
k avec ωk = ( 0, · · · , 0, ωk+1k , · · · , ωnk )
>
et ω>
k ωk = 1.
On vérifie aussi qu’on a Pk = Pk−1 = Pk> , Pk est donc une matrice orthogonale
symétrique. Elle peut aussi s’écrire sous la forme explicite et en blocs:
26
à !
Ik
Pk =
Pbk
Remarque 2.2.1. Ik,l A échange les lignes k et l de A alors que AIk,l échange les
−1 >.
colonnes k et l de A . On a encore Ik,l = Ik,l = Ik,l
1 2 3 0 0 1
Exemple 2.2.1. Soit A la matrice donnée par: A = 4 5 6 et I13 = 0 1 0
7 8 9 1 0 0
7 8 9 3 2 1
alors I13 A = 4 5 6 et AI13 = 6 5 4 .
1 2 3 9 8 7
Définition 2.2.1. Une matrice de permutation est un produit de matrices élémentaires
de permutation.
- Ii (d), une matrice In dont la ieme ligne a été multipliée par un scalaire d;
- Iil (d), une matrice In dont la ieme ligne est multipliée par: ei j + del j , où
¡ ¢
Ii j = ei j i, j=1,··· ,n ;
27
Les matrices Ii j , Ii (d), Iil (d) s’appellent matrices élémentaires de Perlis.
0 0 1 1 0 0 1 d 0
Exemple 2.2.2. I13 = 0 1 0 , I2 (d) = 0 d 0 , I12 (d) = 0 1 0 .
1 0 0 0 0 1 0 0 1
Ces matrices sont régulières et leurs inverses s’écrivent:
Ii−j 1 = Ii j
Ii−1 (d) = Ii (1/d)
Iil−1 (d) = Iil (−d)
Remarque 2.2.2. On peut aussi définir des matrices élémentaires pour les opérations
sur les colonnes, on les note Pi j , Pj (d), Pjl (d).
1 2 3
Exemple 2.2.3. A = 4 5 6 .
7 8 9
Ici la matrice élémentaire M1 de Gauss est donnée par
1 0 0
M1 = −4 1 0 .
−7 0 1
Par ailleurs
1 2 3
A0 = I21 (−4) A = 0 −3 −6
7 8 9
1 2 3
00
et A = I31 (−7) I21 (−4) A = 0 −3 −6
0 −6 −12
1 2 3
(2)
et on a A = M1 A = I31 (−7) I21 (−4) A = 0 −3 −6 .
0 −6 −12
28
2.2.6 Matrices élémentaires de Givens (ou de rotation)
−1
Remarque 2.2.3. On vérifie facilement que R>
k,l = Rk,l .
29
(2)
(2) ai2
On suppose que a22 6= 0 et on recommence avec mi2 = (2)
, i = 3, · · · , n.
a22
A l’étape k, le système se présente sous la forme :
(1) k
(1) (1)
a11 x1 + · · · · · · · · · · · · + a1n xn = b1
.. ..
.
0
.
.. (k) (k) (k)
k ( Sk ) . akk xk + · · · + akn xn = bk
..
.. ..
. . .
(k) (k) (k)
0 ank xk + · · · + ann xn = bn
(k)
(k) aik
En supposant akk 6= 0 , on pose mik = (k)
, i = k + 1, · · · , n puis on remplace
akk
la ligne Li par la ligne Li0 = Li − mik ∗ Lk .
On aboutit alors au système
(1final ( Sn ) de la forme :
) (1) (1)
a11 x1 + · · · · · · · · · · · · + a1n xn = b1
.. ..
.
0
.
.. (k) (k) (k)
( Sn ) ⇐⇒ Ux = c ⇐⇒ . 0 + akk xk + · · · + akn xn = bk
.
.. .. ..
. .
(n) (n)
0··· 0+······ + ann xn = bn
³ ´
(1) (2) (n) >
où c = b1 , b2 , · · · , bn .
30
l’applique avec ou sans pivot.
et (
10−50 x1 + x2 =0
( S2 )
(−1 − 1050 ) x2 = 10−50
qui donne pour solution approchée x1 ' 1 et x2 ' 0.
10−50
Pour lequel m21 = = 10−50 et qui conduit à la solution approchée:
1
x2 ' 1 et x1 = x2 .
A travers cet exemple simple, on voit donc le problème que peut poser un pivot
trop petit. Pour éviter de diviser par des pivots trop petits pouvant conduire à des
solutions absurdes, on peut adopter automatiquement la stratégie du pivot partiel
de la manière suivante : ¯ ¯
(k) (k) ¯ (k) ¯
A chaque étape k : choisir akk tel que : akk = max ¯ aik ¯.
i ≥k
Matriciellement, cette opération revient à multiplier la matrice A(k) par une matrice
de permutation Ikl avant d’appliquer l’élimination de Gauss. La méthode de Gauss
avec pivot partiel s’écrit donc :
où les Mi sont des matrices élémentaires de Gauss et les Iki des matrices de per-
mutation pour i ≥ k. Si à une étape k on n’a pas besoin de pivoter, l’écriture reste
valable avec Iki = I où I désigne la matrice identité.
Théorème 2.3.1. Soit A une matrice carrée, inversible ou non. Il existe (au moins) une
matrice inversible M tellle que la matrice MA soit triangulaire supérieure.
Preuve:
Si A est inversible, le résultat est déjà prouvé en appliquant la méthode de Gauss
31
sans pivot et en posant M = Mn−1 · · · M1 . Si A n’est pas inversible cela signifie
(k)
qu’à une certaine étape k on trouve aik = 0 pour tout i ≥ k. Mais dans ce cas, il
suffit de passer à l’étape suivante.
Matriciellement, cela reviendrait à prendre Iik = Mk = I.
A(k) par deux matrices de permutation P et Q, l’une à droite pour permuter les
lignes et l’autre à gauche pour permuter les colonnes .
C’est une variante qui ressemble à la méthode de Gauss sauf qu’elle aboutit
directement à une matrice diagonale. Au lieu des matrices Mk élémentaires on
considère les matrices
k
1 −m1k
..
0 ... .
.
0 . . −m
k−1k 0
fk = 0 k
M 1
.. .
. −mk+1k . .
.. .. ..
. . .
0 −mk+1k 1
2.4 Factorisation LU
Si on suppose que la méthode de Gauss sans pivot a été appliquée à toutes les
étapes du procédé on aboutit à
32
Remarque 2.4.1. Si on utilise des permutations, alors les matrices Mk Iik ne sont
plus triangulaires inférieures. On démontre dans ce cas qu’on aboutit à la forme
PA = LU.
Preuve:
Si a11 6= 0, la matrice a11 (d’ordre 1) est inversible, donc on peut choisir la matrice
de permutation égale à l’identité et appliquer la méthode de Gauss sans pivot à la
première étape. Supposons qu’on ait pu choisir toutes les matrices de permutation
égales à l’identité jusqu’à l’étape k, il s’ensuit que
i =1
A(k) = Mk−1 Mk−2 · · · M1 A = ∏ Mi A.
i =k −1
Avec
(k)
a11
..
0 .
(k) (k)
Ak =
akk ··· akn
.. .. ..
0 . . .
(k) (k)
ank ··· ann
(k) (k)
1 a11 · · · a1k
. ¶³
.. .. ..
∗∗ . .
B
(k) µ´ k
(k) (k)
= ∗∗ ∗ 1 ak1 · · · akk · · · akn
.. .. .. .. .. .. .. .. ..
. . . . . . . . .
(k) (k) (k)
∗∗ ∗ ··· ··· 1 an1 · · · ank · · · ann
33
Unicité
Supposons qu’il existe L1 , L2 , U1 et U2 telles que A = L1 U1 = L2 U2 , comme L2 et
U1 sont inversibles alors L− 1
2 L1 = U2 U1 .
−1
Ce qui impose L− 1 −1
2 L1 = U2 U1 = I et donc L1 = L2 et U1 = U2 .
Lk−1 l = v (2.5.1)
Lk−1 L>
k−1 = Ak −1 (2.5.2)
> 2
l l + lkk = akk (2.5.3)
34
ii) L’équation (2.5.3) permet d’obtenir la dernière inconnue du problème, à savoir
p
lkk = akk − l > l et on peut choisir lkk > 0.
ω>
0 ω0 = 1. (2.6.1)
P0 a = ke1 , (2.6.2)
¡ ¢1/2
Soit k = ± a> a , les équations (2.6.1) et (2.6.2) donnent :
P0 a = a − 2ω0ω0 a = ke1 et parsuite 2ω0ω>
> >
0 a = − ke1 + a = v, si on pose α = 2ω0 a.
On obtient αω0 = v, et comme on cherche ω0 tel que ω> 2 >
0 ω0 = 1, il vient : α = v v.
2 vv>
Par suite P0 = I − 2 vv> = I − 2 > .
α v v
Remarques 2.6.1. i) Le choix de k se fait au signe près , on peut choisir le signe +.
ii) Le même procédé peut être appliqué pour obtenir une matrice
Pk = Ik − 2ωkω> >
k avec ωk = ( 0, · · · , 0, ωk+1k , · · · , ωnk ) .
35
iii) La factorisation de Householder permet d’écrire :
Pn−2 Pn−3 · · · P1 P0 A = U,
2.7 Conditionnement
Exemple 2.7.1. [dû à R.S. Wilson]
Soit à résoudre le système linéaire Ax = b avec
10 7 8 7 32
7 5 6 5
A= , b = 23 .
8 6 10 9
33
7 5 9 10 31
La solution exacte est donnée par x = (32, 23, 33, 31)> , si on perturbe le second
membre d’un δ b, quel est l’effet de cette perturbation sur la solution?
Soit b + δ b = (32.1, 22.9, 33.1, 30.9)> , alors il en résulte une solution
xe = x + δ x = (9.2, −12.6, 4.5, −1.1).
1
Soit une erreur de l’ordre de sur les données et un rapport d’amplification
200
de l’erreur relative de l’ordre 2000 ce qui montre que le système n’est pas stable,
il reste vulnérable à toute perturbation, on caractérise la stabilité intrinsèque d’un
tel système indépendamment de la méthode de résolution en disant qu’il est mal
conditionné. De même si on perturbe les éléments ai j de la matrice A de δ A, on
aboutit à une solution approchée xe qui est très différente de la solution exacte.
Définition 2.7.1. Soit k.k une norme subordonnée et A une matrice inversible . On
appelle conditionnement de A relativement à la norme k.k, le nombre k Akk A−1 k
noté C ( A) ou cond( A).
Propriétés du conditionnement
i) C ( A) ≥ 1.
¡ ¢
ii) C ( A) = C A−1 .
iii) C (α A) = |α |C ( A), α 6= 0.
° ° µn ( A)
iv) C2 ( A) = k Ak2 ° A−1 °2 = , µi valeur singulière de A.
µ1 ( A)
36
v) C ( QA) = C ( A) pour toute matrice orthogonale Q.
Théorème 2.7.1. Soit A une matrice inversible, x et x + δ x les solutions respectives de
Ax = b et A xe = b + δ b où xe = x + δ x, alors on a
µ ¶ µ ¶
1 kδ bk kδ xk kδ bk
≤ ≤ C ( A) . (2.7.1)
C ( A) kbk k xk kbk
De même la solution obtenue après perturbation de A par δ A vérifie :
kδ xk C ( A) kδ Ak
≤
.
k xk kδ Ak k Ak
1 − C ( A)
k Ak
Si on perturbe en même temps A et b alors on obtient :
µ ¶
kδ xk C ( A ) kδ Ak kδ bk
≤
.
k xk kδ Ak k Ak kbk
1 − C ( A)
k Ak
Preuve:
On a Ax = b et Ax + Aδ x = b + δ b d’où Aδ x = δ b ou encore δ x = A−1 δ b, comme
kδ xk kδ bk
kbk ≤ k Akk xk, on déduit que ≤ C ( A) .
k xk kbk
Les autres inégalités s’obtiennent de façons similaires.
Remarques 2.7.1. i) L’équation (2.7.1) donne une estimation de l’erreur relative
kδ bk
de la solution en fonction de l’erreur relative connue .
kbk
ii) Tous les calculs sont effectués sur un ordinateur, des erreurs d’arrondi sont
accumulées à chaque étape de calcul. Si ε désigne la précision numérique
relative (dépend de la machine), l’erreur relative de la solution explose si
C ( A) × ε ' 1.
µ ¶n
1
iii) La matrice de Hilbert d’ordre n, H = présente un exemple
i + j − 1 i, j=1
classique de matrices mal conditionnées.
37
- tridiagonales: les éléments non nuls sont sur la diagonale et de part et d’autre
de celle-ci sur les deux lignes adjacentes. On a ai j = 0 pour |i − j| > 1,
Matrices tridiagonales
Soit le système
d1 e1
x1 b1
c2 d2 e2 .. ..
. .
.. .. ..
. . . .. ..
. = . (2.8.1)
.. .. ..
. . . .. ..
. .
cn−1 dn−1 en−1
xn bn
cn dn
Il y’a un vaste choix de bibliothèques disponibles pour calculer les solutions qui
prennent un temps de calcul proportionnel à n.
De manière générale, on peut obtenir un algorithme proportionnel à n pour une
matrice à bande. Le préfacteur de l’algorithme est proportionnel à m.
Formule de Sherman-Morison
Supposons que l’on ait une matrice A dont on a facilement calculé l’inverse (cas
d’une matrice triangulaire ). Si on fait un petit changement dans A en modifiant
par exemple un ou quelques éléments de l’ensemble de la matrice, peut on calculer
facilement l’inverse de cette nouvelle matrice?
En effet soit B une matrice telle que
B = A + uv> (2.8.2)
38
La formule de Sherman-Morison donne l’inverse de B.
³ ´−1
B−1 = I + A−1 uv> A−1
³ ´
= I − A−1 uv> + A−1 uv> A−1 uv> + · · · A−1
³ ´
= A−1 − A−1 uv> A−1 1 − λ + λ 2 + · · ·
A−1 uv> A−1
= A−1 −
1+λ
¡ ¢>
où λ = v> A−1 u. Posons z = A−1 u et w = A−1 v on a λ = v> z et
zw>
B−1 = A−1 − .
1+λ
Remarque 2.8.1. Dans le cas où les vecteurs u et v sont des matrices U et V d’ordre
respectif n × k et k × n et si de plus la matrice I + VA−1 U est inversible, alors
et que
ci ei−1 6= 0, i = 2, · · · , n (2.8.4)
Sous ces conditions la matrice A est irréductible à diagonale dominante donc A est
inversible. Wendroff a montré qu’il existe une factorisation de la matrice A sous la
forme suivante
A=b LRb (2.8.5)
b est une matrice triangulaire supérieure avec éléments diagonaux égaux à 1 et
où R
b
L est une matrice triangulaire inférieure. Si A est tridiagonale, les matrices b b
R et L
α1 1 γ1
.. ..
c2 α2 . .
sont données par b
L= .. ..
et R
b= ,
. . . .. γ
n
cn αn 1
où
e1
i ) α1 = d1 , γ1 =
d1
ii ) αi = di − ci γi−1 , i = 2, · · · , n (2.8.6)
ei
iii ) γi = , i = 2, · · · , n − 1
αi
39
n
Puisque A est inversible et det( A) = ∏ αi , alors tous les αi sont non nuls et le
i =1
schéma récursif (2.8.6) est bien défini.
Remarque 2.8.2. On note que la factorisation (2.8.5) est obtenue par application de
la méthode d’élimination de Gauss usuelle à A> .
La résolution du système lnéaire est équivalente à
b b = v,
Lv = b, Rx
i ) |γi | < 1
(2.8.7)
ii ) 0 < |di | − |ci | < |di | + |ci |, i = 2, · · · , n
Preuve:
ei e |e |
On a γi = = i , la condition i) de (2.8.3) donne |γ1 | < 1 < 1.
αi di |d1 |
Supposons que γ j < 1 pour j = 1, · · · , i − 1, en remplaçant (2.8.6 (iii)) dans (2.8.6
(ii)) on obtient
ei
γi =
di − ci γi −1
et
| ei |
|γi | <
|di | − |ci ||γi−1 |
| ei |
<
| di | − | ci |
par hypothèse. Ensuite en considérant (2.8.3 (ii)) il s’ensuit que |γi | < 1.
Pour montrer (2.8.7 (ii)) il suffit de considérer (2.8.6 (ii)) et (2.8.3 (ii)).
40
L’algorithme présenté précédemment peut être simplifié en éliminant les αi ,
c’est le cas où la factorisation b LR b n’est pas nécessaire. L’algorithme se présente
comme suit
v1 = b1 /d1
t = d1
Pour i = 2 à n
γi −1 = ei −1 / t
t = di − ci γi −1
vi = (bi − ci vi−1 )/t
fin
xn = vn
Pour i = 1 à n-1
x n −i = v n −i − γ n −i x n −i + 1
fin
Si A est définie positive, elle possède une factorisation de la forme L> DL qui est
obtenue facilement du schéma (2.8.6). Si
1
..
l .
1
L= (2.8.8)
.. ..
. .
ln−1 1
41
Pour i = 1 à n-1
xn−i = (vn−i − en−i xn−i+1 )/δn−i
fin
Si la factorisation n’est pas nécessaire on obtient l’algorithme simplifié suivant
v1 = b1
δ1 = d1
Pour i = 1 à n-1
t = ei /δi
δi+1 = di+1 − tei
vi+1 = bi+1 − tvi
fin
xn = vn /δn
Pour i = 1 à n-1
xn−i = (vn−i − en−i xn−i+1 )/δn−i
fin
Exemple 2.8.1. Soit le système linéaire Ax = b, où A est une matrice carrée d’ordre
n de la forme
a1 · · · · · · · · · · · · an
1 λ 1
.. .. ..
. . .
A=
.
.. .. ..
. . .
1 λ 1
b1 · · · · · · · · · · · · bn
Ce type de matrices est issu de la discrértisation de certains problèmes paraboliques
voir [32, 130].
Considérons la matrice
α 1
1 α 1
.. .. ..
A0 = . . . .
1 α 1
1 α
√ √
λ− λ2 − 4 λ + λ2 − 4
avec α = et β = . Il s’ensuit que
2 2
A0 = LU
42
avec
α 1 β
.. .. ..
1 . . .
.. .. .. ..
L=
. . , U =
. . .
.. .. ..
. . . β
1 α 1
La solution y de A0 y = ω, est donnée par
yn = vn
ym = vm − β ym+1 , m = 0, · · · , n − 1,
vn = βω0
vm = β(ωm − ωm−1 ), m = 1, · · · , n.
Soit z = x − y, on obtient
à !>
n n
Az = ω0 − ∑ ak yk , · · · , ωn − ∑ bk yk .
k=0 k=0
Comme
zm−1 + λ zm + zm+1 = 0, m = 1, · · · , n − 1.
Il vient
zm = c0 γ n−m + c1 γ m , m = 0, · · · , n,
√
−λ − λ2 − 4
où γ = et c0 , c1 sont des constantes à déterminer en résolvant le
2
système
n n n
c0 ∑ ak γ n−k + c1 ∑ ak γ k = ω0 − ∑ ak yk
k =0 k=0 k=0
n n n
c0 ∑ b k γ n−k + c 1 ∑ b k γ k = ωn − ∑ bk yk
k =0 k =0 k=0
43
matrice rectangulaire (m, n) donnée par :
n
µ1
..
.
µr
0 m>n
Σ=
..
.
n
0
0
où les µi sont les valeurs singulières de A.
Corollaire 2.9.1. • Le rang de A est égal au nombre de valeurs singulières non nulles.
Preuve:
e1 · · · en e1 · · · em
1 ··· 0 1 ··· 0
. .. . ..
Soient In = .. . et
Im = .. .
0 1 0 1
on a par définition
n n
V= ∑ vi ei∗ , V∗ = ∑ ei vi∗ , Σei = µiεi pour i = 1 · · · n,
i =1 i =1
n r r
ΣV ∗ = ∑ µiεi vi∗ , UΣV ∗ = ∑ µi ui vi∗ et A∗ A = ∑ µi2 ui vi∗ . .
i =1 i =1 i =1
44
Cas où m ≥ n et rang A = n
n
ri = ∑ ai j x j − bi i = 1, · · · , m.
j=1
( Vectoriellement r = Ax − b).
La méthode des moindres carrées consiste à minimiser krk22 .
Posons
I ( x) = r> r = ( Ax − b)> ( Ax − b).
∂I ( x)
= 0.
∂xi
D’autre part
∂I ( x)
= ei> A> Ax + xA> Aei − 2ei> A> b
∂xi
= 2ei> ( A> Ax − A> b).
46
8.1 0 0 0 0 0 0 0 0
0 0.64 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
w=
0
,
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0.15 −0.98 0 0 0 0 0 0 0
0.34 0.05 0.93 0 0 0 0 0 0
0.34 0.05 −0.13 0.92 0 0 0 0 0
0.34 0.05 −0.13 −0.15 0.91 0 0 0 0
V=
0.34 0.05 −0.13 −0.15 −0.18 −0.89 0 0 0
0.34 0.05 −0.13 −0.15 −0.18 0.22 −0.5 −0.5 −0.5
0.34 0.05 −0.13 −0.15 −0.18 0.22 0.83 −0.16 −0.16
0.34 0.05 −0.13 −0.15 −0.18 0.22 −0.16 0.83 −0.16
0.34 0.05 −0.13 −0.15 −0.18 0.22 −0.16 −0.16 0.83
n
2. AA+ = ∑ ui ui∗ matrice de projection orthogonale sur ImA.
i =1
n
3. A+ A = ∑ vi vi∗ matrice de projection orthogonale sur ImA∗ .
i =1
i) Si A est une matrice triangulaire , le système est résolu par simple sub-
stitution.
47
ii) Si A est une matrice symétrique ou hermitienne, définie positive , la
résolution est effectuée par la méthode de Choleski.
iii) Si A est une matrice carrée, mais n’entrant pas dans les deux cas précé-
dents, une factorisation LU est réalisée en utilisant la méthode d’élimina-
tion de Gauss avec stratégie de pivot partiel.
iv) Si A n’est pas une matrice carrée, la méthode QR est utilisée.
4. La matrice A peut être creuse, elle comporte une forte proportion de co-
efficients nuls (de nombreux problèmes issus de la physique conduisent à
l’analyse de systèmes linéaires à matrices creuses ), l’intérêt de telles matrices
résulte non seulement de la réduction de la place mémoire (on ne stocke pas
les zéros) mais aussi de la réduction des nombres d’opérations. Dans le cas
de ces matrices des algorithmes particuliers sont mis en œuvre.
2.11 Conclusion
Avant d’entamer la résolution des systèmes linéaires de grandes tailles, il est
impératif de commencer par une analyse des propriétés de la matrice afin de détermi-
ner la méthode la plus adaptée afin d’obtenir une solution avec une précision cor-
recte et pour un temps de calcul qui sera minimal. Les différentes méthodes dans
ce chapitre ne sont qu’une introduction à ce très vaste sujet.
48
2.12 Exercices
1. Montrer que la matrice I − α ab> admet une matrice inverse de la même forme
en imposant une condition sur α . En déduire l’expression de ( A + ab> )−1 en
fonction de A−1 et des vecteurs a et b.
3. Montrer que:
( I + m1 e> > > > >
1 ) . ( I + m 2 e2 ) · · · ( I + m n−1 en−1 ) = I + m1 e1 + · · · + mn−1 en−1
2. Déduire le rang de A
49
2 1 0 0
1 5 3 4
Exercice 2.12.4. Soit A =
0 3 −20 15
0 4 15 −5
50
Chapitre 3
3.1 Introduction
Pour résoudre le système
Ax = b, (3.1.1)
• Méthode de Jacobi: M = D, N = L + U.
51
• Méthode de Gauss-Seidel: M = D − L, N = U.
1
• Méthode de relaxation : A = A(ω) = M(ω) − N (ω), avec M(ω) = D − L,
ω
1−ω
N (ω) = D − U où ω est un scalaire.
ω
La deuxième sera consacrée aux matrices positives.
Théorème 3.2.1. Soit A une matrice carrée. Les assertions suivantes sont équivalentes :
i) lim Ak = 0.
k−→+∞
iii) ρ( A) < 1.
52
l’exercice 1.5.3 assure que k Ak ≤ ρ( A) + ε et par suite k Ak < 1.
iv)=⇒ i)
Elle est évidente car si k Ak < 1 alors lim k Akk = 0 et par ailleurs k Ak k ≤ k Akk
k−→+∞
d’où lim k Ak k = 0 et lim Ak = 0
k−→+∞ k−→+∞
Théorème 3.2.2. Soit A une matrice carrée et k.k une norme quelconque,
° °1/k
° °
alors lim ° Ak ° = ρ( A).
k−→+∞
Preuve:
¡ ¢
i) On a ρ( A) ≤ k Ak et comme ρ Ak = (ρ( A))k
¡ ¢ ° ° ° °1/k
il s’ensuit que ρ Ak = (ρ( A))k ≤ ° Ak ° et par suite ρ( A) ≤ ° Ak ° donc
° °1/k
° °
ρ( A) ≤ lim ° Ak ° .
k−→+∞
° °1/k
° °
ii) Pour montrer que lim ° Ak ° ≤ ρ( A), introduisons pour tout ε > 0 la ma-
k−→+∞
1 ρ( A)
trice Aε = A, pour cette matrice on a ρ( Aε ) = < 1 et d’après
ρ( A) + ε ρ( A) + ε
le théorème précédent on obtient : k Aε k < 1 pour au moins une norme matricielle
induite, donc lim Aεk = 0 et lim k Aε kk = 0.
k−→+∞ k−→+∞ ° °
Donc ∀ ε > 0, ∃ N (ε), tel que pour tout k ≥ N (ε) on ait: ° Aεk ° < 1
° ° ° °1/k
ou encore ° Ak ° ≤ (ρ( A) + ε)k soit encore ° Ak ° ≤ (ρ( A) + ε) pour tout ε > 0
° °1/k
° k°
finalement lim ° A ° ≤ ρ( A).
k−→+∞
53
Définition 3.2.4. On appelle taux moyen de convergence sur k itérations le nombre
° °
e (k, T ) = − log ° T k °1/k et taux asymptotique de convergence le nombre
R
R( T ) joue le rôle de vitesse de convergence, plus R( T ) est grand plus rapide est la
convergence.
x(k+1) = D −1 ( L + U ) x(k) + D −1 b.
Théorème 3.3.1. Si A est une matrice carrée à diagonale strictement dominante en lignes
alors la méthode de Jacobi converge.
Preuve:
On a
n
∑ |ai j | < |aii | i = 1, · · · , n (3.3.1)
j=1
j 6 =i
54
ai j
d’autre part on a : ti j = − pour i 6= j et tii = 0 d’où
aii ½ ¾
1
= max ∑ |ti j | = max
| aii | ∑
k TJ k∞ | ai j | et d’après (3.3.1) on a k TJ k∞ < 1.
i j i j 6 =i
Preuve:
On procède d’une manière analogue à celle de la preuve du théorème 3.3.1 pour
montrer que k TJ k∞ ≤ 1.
Donc
| TJ |e ≤ e, | TJ |e 6= e, e = (1, 1, · · · , 1)> . (3.3.2)
| TJ |n e ≤ e,
| TJ |2 e ≤ | TJ |e < e,
et par suite
| TJ |i+1 e ≤ | TJ |i e ≤ · · · < e,
55
Montrons que le nombre de composantes non nulles τi de t(i) croı̂t avec i.
Si ce n’était pas le cas, (3.3.3) impliquerait qu’il existe i ≥ 1 tel que τi = τi+1 .
t(i) peut s’écrire sous la forme
à !
a
, a > 0, a ∈ R p .
0
Il s’ensuit que
à !
b
= t(i+1) = e − | TJ |i+1 e ≥ | TJ |e − | TJ |i+1 e
0
à !à !
(i ) | T11 | | T12 | a
= | TJ |t =
| T21 | | T22 | 0
Puisque a > 0, ceci n’est possible que si T21 = 0 ie TJ est réductible. Ceci contredit
les hypothèses du théorème.
Donc 0 < τ1 < τ2 < · · · , et t(n) > 0. La preuve est ainsi achevée.
Pour remédier au problème du stockage et dans l’espoir d’améliorer les résultats
en accélérant la convergence, on cherche une méthode qui utilise les composantes
de x(k+1) au fur et à mesure qu’elles sont calculées. C’est ce que réalise la méthode
de Gauss-Seidel.
ou encore
x(k+1) = ( D − L)−1 Ux(k) + ( D − L)−1 b (3.3.5)
56
en supposant que D − L est inversible .
Les équations (3.3.4) et (3.3.5) peuvent aussi être présentées sous les formes :
Théorème 3.3.3. Si A est une matrice carrée à diagonale strictement dominante en lignes
alors la méthode de Gauss-Seidel converge .
Preuve:
k Txk∞
Posons T = ( D − L)−1 U et montrons que k T k∞ < 1 où k T k∞ = max .
x6=0 k x k∞
Soit y = Tx = ( D − L)−1 Ux on a alors ( D − L) y = Ux ou encore Dy = Ly + Ux et
y = D −1 Ly + D −1 Ux. Considérons l’indice i0 tel que
Il vient :
i0 −1 n
yi0 = ∑ ( D −1 L )i0 j y j + ∑ ( D −1 U )i0 j x j .
j=1 j =i 0 + 1
Par suite
i0 −1 ¯¯a
¯
¯ n
¯ ¯
¯ ai0 j ¯
| yi0 | = k y k ∞ ≤ ∑ ¯ i0 j ¯ k yk∞ + ¯
∑ ¯ ai i ¯¯ kxk∞ .
¯a ¯
j=1 i0 i0 j =i 0 + 1 0 0
57
En regroupant les termes
à ¯!
i0 −1 ¯¯ ai0 j ¯ k y k ∞ n
¯ ¯
¯ ai0 j ¯
1− ∑ ¯ ¯ ∑ ¯¯ ai i ¯¯
¯ a ¯ k xk∞ ≤
j=1 i0 i0 j =i 0 + 1 0 0
i0 −1 ¯¯a
¯
¯
Par hypothèse, le terme 1 − ∑ ¯ i0 j ¯ est strictement positif d’où on en tire :
¯a ¯
j=1 i0 i0
à ¯ ¯! à ¯ !−1
n ¯ ai0 j ¯ i0 −1 ¯¯ a ¯
k Txk∞ k yk∞ i j
k xk∞
=
k xk∞
≤ ∑ ¯¯ ai i ¯¯ 1 − ∑ ¯¯ ai 0i ¯¯ ,
j =i0 + 1 0 0 j=1 0 0
finalement
k Txk∞
max < 1.
x6=0 k xk∞
Remarques 3.3.2. 1. Un résultat de convergence similaire a lieu si A est à diag-
onale dominante en colonnes.
58
• Si ω < 1, on parle de sous-relaxation.
Ici la condition de convergence k Tω k < 1 dépendra du paramètre ω et par
conséquent, on est amené à chercher tous les ω pour lesquels il y a convergence et
ensuite choisir la valeur optimale ω0 de telle sorte que la vitesse de convergence
soit la meilleure possible.
Preuve:
A = A∗ =⇒ ( M − N )∗ = M∗ − N ∗ = M − N.
59
Par suite R est hermitienne définie positive.
En remarquant que A = R + T ∗ AT, on obtient
A = R + T ∗ ( R + T ∗ AT ) T = R + T ∗ RT + T ∗2 ( R + T ∗ AT ) T 2
k−1
= ∑ T ∗i RTi + T ∗k RT k .
i =0
lim T k = lim T ∗k = 0.
k−→∞ k−→∞
Donc
∞ ∞
A= ∑ T ∗i RTi = R + ∑ T ∗i RTi .
i =0 i =0
Preuve:
D’après le théorème précédent, si A est hermitienne définie positive et M∗ + N
définie positive alors la méthode converge. Il suffit donc que M∗ + N soit définie
positive.
2−ω 2−ω
Or M∗ + N = D et par suite M∗ + N est définie positive si > 0 c.à.d
ω ω
si 0 < ω < 2.
Preuve:
µ ¶−1 µ ¶
1 1−ω
Tω = D−L D+U
ω ω
60
Si les valeurs propres de Tω sont notées λi (ω) on a :
µ ¶
1−ω
n det D+U
ω
det Tω = ∏ λi (ω) = µ ¶ = (1 − ω)n .
i =1
1
det D−L
ω
1
D’où ρ( Tω ) ≥ [|1 − ω|n ] n= |1 − ω|.
Pour que la méthode converge, il est nécessaire d’avoir ρ( Tω ) < 1 et par conséquent
|1 − ω| < 1 d’où ω ∈]0, 2[.
Remarque 3.3.4. Le résultat reste valable si A est tridiagonale par blocs (voir Ciar-
let(1984) ou Sibony(1986)).
Théorème 3.4.1. Soit A une matrice tridiagonale. Alors les méthodes de Jacobi et de
Gauss-Seidel convergent ou divergent simultanément; lorsqu’elles convergent, la méthode
de Gauss-Seidel est plus rapide que celle de Jacobi, plus précisément, on a
ρ( TGS ) = (ρ( TJ ))2 où TGS et TJ sont les matrices de Gauss-Seidel et Jacobi (respective-
ment).
Preuve:
TGS = ( D − L)−1 U et TJ = D −1 ( L + U )
On est donc amené à chercher les valeurs propres de ces deux matrices.
1) Soit λ une valeur propre de TJ . Alors λ est racine du polynôme caractéristique
PJ (λ )
³ ´
PJ (λ ) = det D −1 ( L + U ) − λ I
³ ´
= det − D−1 det(λ D − ( L + U ))
= K1 det(λ D − ( L + U ))
2) Soit α une valeur propre de TGS . Alors α est racine du polynôme caractéristique
PGS (α ).
³ ´
−1
PGS (α ) = det ( D − L) (U ) − α I
³ ´
= det −( D − L)−1 det(α D − α L − U )
= K2 det(α D − α L − U ).
61
Comme il s’agit de matrices tridiagonales, on a le résultat suivant (voir exercice
(3.7.2)) ³ ´
PGS (α ) = K2 det α D − αµ L − µ −1 U pour tout µ 6= 0.
Théorème 3.4.2. Soit A une matrice tridiagonale, telle que toutes les valeurs propres de
la matrice de Jacobi associée soient réelles. Alors les méthodes de Jacobi et de relaxation,
pour ω ∈]0, 2[, convergent où divergent simultanément. Lorsqu’elles convergent, on peut
déterminer une valeur optimale ω0 du paramètre ω telle que
2
ω0 = q , ρ( Tω0 ) = inf ρ( Tω ) = ω0 − 1.
ω∈]0,2[
1+ 1 − ρ( TJ )2
0.95
0.9
0.85
0.8
ρ(Tω)
0.75
0.7
0.65
0.6
0.55
ω0 − 1
0.5
0 0.2 0.4 0.6 0.8 1
ω
1.2 1.4 ω1.60 1.8 2
62
Remarques 3.4.1. 1. ρ( TGS ) = ρ( T1 ) est obtenue comme cas particulier avec
ω = 1 on sait d’après le théorème précédent que ρ( TGS ) = (ρ( TJ ))2 .
PJ (α ) = K1 det(α D − L − U )
µ ¶
1−ω−λ
et Pω (λ ) = Kω det D + λL + U
ω
µ ¶
n λ +ω−1
= (−1) Kω det D − λL − U
ω
à !
2
eω det λ + ω − 1
Pω (λ 2 ) = K D − λ2 L − U
ω
à !
2
eω det λ + ω − 1
= K D − λ L − λU , λ 6= 0
ω
Soit encore :
(ω − 1)
λ+
eω λ n det
Pω (λ 2 ) = K λ
D − L − U pour tout λ 6= 0
ω
Il s’ensuit qu’on peut déterminer une relation liant les valeurs propres non nulles
de Tω et celles de TJ .
λ2 + ω − 1
Soit λ 2 6= 0 une valeur propre de Tω donc ± est une valeur propre de TJ .
λω
Inversement, si α est valeur propre de TJ , on sait que −α est aussi valeur propre de
TJ et ceci implique que λ12 et λ22 sont valeurs propres de Tω où λ1 et λ2 sont données
63
par :
1³ 2 2 ´ αω q
=λ12 α ω − 2ω + 2 − (α 2ω2 + 4 − 4ω)
2 2
1³ 2 2 ´ αω q
λ22 = α ω − 2ω + 2 + (α 2ω2 + 4 − 4ω)
2 2
λ1 et λ2 proviennent de la résolution de l’équation du second degré en λ, à savoir
λ2 + ω − 1
α=
λω
ou encore
λ 2 − αωλ + ω − 1 = 0 (∗)
on a
∆(α ) = α 2ω2 − 4ω + 4 = α 2ω2 − 4(ω − 1) = α 2ω2 + 4(1 − ω)
2. Si |α | < 1
On est amené à étudier
n ³¯ ¯ ¯ ¯´o
¯ 2 ¯ ¯ 2 ¯
ρ( Tω ) = max max ¯λ1 (α , ω)¯ , ¯λ2 (α , ω)¯
α ∈ Sp( TJ )
64
pour cela, considérons la fonction f définie par :
f : R+ ×]0, 2[−→ R+
n¯ ¯ ¯ ¯o
¯ ¯ ¯ ¯
(α , ω) −→ max ¯λ12 ¯ , ¯λ22 ¯
65
p p
en utilisant 2λ2 = αω + ∆(α ) donc ∆(α ) = 2λ2 − αω d’où
à !
∂f α ωα 2 − 2
= + 2λ2
∂ω 2 2(2λ − αω)
à !
2αλ − ωα 2 + ωα 2 − 2
= 2λ
2(2λ − αω)
µ ¶
αλ − 1
= 2λ2
2λ − αω
µ ¶
λ2α − 1
= 2λ2
2λ2 − αω
∂f
Donc < 0 car αλ2 < λ2 < 1 pour 0 < α < 1 en effet
∂ω
q
α 2ω2 − 2ω + 2 αω
λ22 = + (α 2ω2 − 4ω + 4)
2 2
q
ω2 − 2ω + 2 ω
< + (ω2 − 4ω + 2)
2 2
ω2 − 2ω + 2 ω
< + (2 − ω)
2 2
< 1
Par ailleurs
∂f
1- (0) = 0 donc ∆(α )(ω = 0) = 4
∂ω q
√
∆(α )(ω = 0) 4
et λ2 (ω = 0) = = = 1par conséquent
2 2
2
∂f ωλ − 2
(ω = 0) = 2λ2 λ2 + q
∂ω ∆(α )(ω = 0)
2
= 2λ2 (0) λ2 (0) − q
∆(α )(ω = 0)
2
= 2(1 − )
2
= 0
∂f 1
2- (ω0 ) = ∞ car ω0 est racine de ∆(α ) d’où tend vers l’infini
∂ω ∆(α )
quand ω tend vers ω0 .
66
1
0.95
0.9
0.85
0.8
ρ(Tω)
0.75
0.7
0.65
0.6
0.55
ω0 − 1
0.5
0 0.2 0.4 0.6 0.8 1
ω
1.2 1.4 ω1.60 1.8 2
2
α −→ p est croissante en fonction de α
1 + 1 − α2
0.9
0.8
α3
0.7
0.6 α2
ρ(Tω)
0.5 α1
0.4 α0=0
0.3
0.2
0.1
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
ω
67
3.5 Méthodes semi-itératives
Dans la première partie de ce chapitre, nous avons considéré les méthodes itératives
du type x(k+1) = Tx(k) + C, impliquant deux itérés successifs x(k) et x(k+1) .
A présent, nous considérons des méthodes semi-itératives permettant d’exprimer
y(k) comme combinaison algébrique des k itérés précédents x(0) , x(1) , · · · , x(k−1) .
Nous obtenons l’expression:
k
y(k) = ∑ θ j (k ) x( j) , k ≥ 0,
j=0
k
avec ∑ θ j (k) = 1, k ≥ 0.
j=0
Sie( j) =x( j) − x est l’erreur de la jeme itération et ε( j) = y( j) − x est l’erreur de la
méthode semi-itérative alors on a:
k
ε(k) = ∑ θ j ( k ) e( j) ,
j=0
ou encore
ε(k) = Qk ( T ) e(0) , (3.5.1)
k
où Qk ( x) est le polynôme Qk ( x) = ∑ θ j (k) x j .
j=0
1
Remarques 3.5.1. 1. Si on prend θ j (k) = pour j = 0, 1, · · · , k, on retrouve
k+1
y(k) comme moyenne arithmétique des x(k) .
Théorème 3.5.1. Si les coefficients des polunômes Qk ( x) sont réels et si la matrice T est
hermitienne et ses valeurs propres vérifient −1 ≤ a ≤ λ1 ≤ · · · ≤ λn ≤ b ≤ 1, alors
68
En revenant à l’équation (3.5.1), nous sommes donc amenés au problème de
minimisation suivant:
min | Qk ( x)|.
Qk (1)=1
dont la solution est donnée par les polynômes de Chebyshev (voir Varga(2000),
Golub(1989))
(
cos(k cos−1 ( x)), −1 ≤ x ≤ 1, k ≥ 0
Ck ( x) =
cosh(k cosh−1 ( x)), x ≥ 1, k ≥ 0
Théorème 3.6.1. Soit B une matrice carrée; si B ≥ 0 alors les assertions suivantes sont
équivalentes:
i) β > ρ( B).
(β I − B)−1 ≥ 0
Preuve:
Supposons que β > ρ( B).
1
En posant M = (β I − B) = β( I − Bβ ) où Bβ = B, il est évident que ρ( Bβ ) < 1 et
β
il découle du théorème 1.4.1 que B est inversible et (β I − B)−1 ≥ 0.
Réciproquement, en supposant ii) et en utilisant le théorème 1.4.1 avec B ≥ 0, si
v ≥ 0 est un vecteur propre de B associé à la valeur propre ρ( B) alors v est aussi
vecteur propre de (β I − B) associé à la valerur propre (β − ρ( B)).
De l’inversibilité de (β I − B) il s’ensuit que (β − ρ( B)) 6= 0 et par suite
1
(β I − B)−1 v = v.
(β − ρ( B))
Enfin, v étant un vecteur propre ≥ 0 (non identiquement nul) et (β I − B)−1 ≥ 0
entraı̂ne que (β − ρ( B)) > 0.
Théorème 3.6.2. Soit A une matrice réelle avec ai j ≤ 0 pour tout i 6= j et D = ( aii ), on
a équivalence entre:
69
i) A est inversible et A−1 ≥ 0
et
Preuve:
Supposons que A est inversible et A−1 = (αi j ) avec αi j ≥ 0 ∀ i, j = 1, n , en
explicitant A−1 A = I, il vient
(α I − A)−1 > 0
Théorème 3.6.4. Soit A une matrice réelle avec ai j ≤ 0 pour tout i 6= j et D = ( aii ), on
a équivalence entre:
Théorème 3.6.5. Soit A une matrice réelle irreductible à diagonale dominante avec ai j ≤
0 pour tout i 6= j et aii > 0 ∀i = 1, · · · , n alors A−1 > 0.
Théorème 3.6.6. Soit A une matrice réelle symétrique inversible et irreductible avec
ai j ≤ 0 pour tout i 6= j, alors A−1 > 0 si et seulement si A est définie positive.
70
3.6.1 Décomposition régulière des matrices
ρ ( A−1 N )
ρ ( M −1 N ) = <1
1 + ρ ( A−1 N )
et par conséquent, si A est inversible avec A−1 ≥ 0 alors la méthode x(k+1) = M−1 Nx(k) +
c est convergente pour n’importe quel choix initial.
Preuve:
Supposons A monotone.
En posant B = A−1 N, on peut aisement voir que
• M−1 N = ( I + B)−1 B,
λ µ
• si M−1 Nv = λ v alors Bv = µ v avec µ = et λ = ,
1−λ 1+µ
• M−1 N ≥ 0 et B ≥ 0,
µ
• la fonction µ −→ est strictement croissante pour µ ≥ 0 d’où ρ( M−1 N ) =
1+µ
ρ ( A−1 N )
< 1 et par conséquent la méthode converge.
1 + ρ ( A−1 N )
Réciproquement, si ρ( M−1 N ) < 1 alors( I − M−1 N ) est inversible et on a succes-
sivement
M−1 A = I − M−1 N et A−1 = ( I − M−1 N )−1 M−1 ,
Preuve:
On a A−1 N1 ≤ A−1 N2 et par suite ρ( M1−1 N1 ) ≤ ρ( M2−1 N2 ) ≤ 1.
71
3.7 Comparaison des méthodes classiques dans le cas des
matrices positives
Soit la décomposition A = D − E − F où D = ( aii )i=1,··· ,n et E etF sont respective-
ment triangulaire inférieure et supérieure. En supposant D et ( D − E) inversibles
et en posant L = D −1 E et U = D −1 F, on retrouve les matrices de Jacobi TJ = L + U
et de Gauss-Seidel TGS = ( I − L)−1 U qu’on peut comparer dans certains cas parti-
culiers.
1. ρ( TJ ) = ρ( TGS ) = 0,
3. ρ( TJ ) = ρ( TGS ) = 1,
R( TJ ) < R( TGS ).
Remarques 3.7.1. 1. Le théorème de Stein & Rosenberg affirme donc que les
matrices de Jacobi et de Gauss-Seidel convergent ou divergent simultanément.
Lorsqu’elles convergent, la matrice de Gauss-Seidel converge asymptotique-
ment plus rapidement.
M1 = I; N1 = L + U,
M2 = I − L; N2 = U.
Le théorème 3.6.8 permet d’obtenir directement le point 2 du théorème de
Stein & Rosenberg.
Mω = ( I − ω L)−1 (ωU + (1 − ω) I )
72
est convergente pour 0 < ω ≤ 1. De plus, si 0 < ω1 < ω2 ≤ 1, alors
73
3.9 Exercices
Exercice 3.9.1. Soit B une approximation de l’inverse A−1 d’une matrice carrée A
et x̂ une solution approchée du système Ax = b.
En posant R = I − BA et r = b − A x̂, montrer que si k Rk < 1 alors on a:
° ° k Bk
i) ° A−1 ° ≤
1 − k Rk
° ° k Bk.k Rk
ii) ° A−1 − B° ≤
1 − k Rk
k B k . kr k
iii) k x − x̂k ≤
1 − k Rk
Exercice 3.9.2. Soit θ un scalaire non nul et A(θ ) la matrice tridiagonale suivante:
a1 b1θ −1
c2θ a 2 b 2θ −1
c3θ .
A(θ ) =
.
. bn−1θ −1
cnθ an
k Ax − Ayk ≤ α k x − yk ∀ x ∈ E, ∀ y ∈ E;
avec 0 ≤ α < 1.
1. Soit x0 ∈ E donné, montrer que l’itération xn+1 = Axn définit une suite de
Cauchy dans E. On notera x∗ la limite de cette suite .
74
a) Montrer que l’équation x = Mx + b( x) admet une solution unique x∗
b) Montrer que si on utilise l’itération xn+1 = Mxn + b( xn ) alors on obtient
k x∗ − xn+1 k∞ ≤ k k xn+1 − xn k∞ ,
Exercice 3.9.5. Soit A une matrice symétrique définie positive et T la matrice définie
par: T = 2D −1 − D −1 AD −1 . On suppose que 2D − A est définie positive.
Exercice 3.9.6. Soit A une matrice carrée d’ordre n à coefficients dans Cn . Pour
chercher la solution du système Ax = b on considère le schéma itératif suivant:
(
x(n+1) = ( I − α A ) x(n) + α b
(∗)
x(0) ∈ Cn donné
avec α > 0.
II- On suppose que A est hermitienne définie positive dont les valeurs propres
sont rangées par ordre décroissant: 0 ≤ λn ≤ λn−1 ≤ · · · ≤ λ1 .
75
1. Quelle condition doit vérifier α pour que la méthode converge?
2. Déterminer la valeur optimale de α .
2. Donner
un algorithme de décomposition
A = LU oùAest donnée par:
a1 d2 1 u1 d2
c1 a2 d3 l1 1 u2 d3
c2 . l2 . .
A=
. . .
. d . . d
n n
cn−1 an ln−1 1 un
en posant: cn+1 = d1 = l1 = 0 et en supposant que | ai | > |ci+1 | + |di |
Montrer que: |li | < 1 et |ui | > |ci+1 |
a b b
Exercice 3.9.8. Soit A la matrice donnée par A = b a b avec a > b > .0
b b a
2. Donner une condition sur a et b pour que la méthode de Jacobi soit conver-
gente.
3. On suppose que a > 2b, montrer que A est inversible que le conditionnement
α a + βb
de A pour la norme k k∞ vérifie C ( A) ≤ où α et β sont des constantes
α a − βb
à déterminer.
Exercice 3.9.9. Soit A une matrice carrée hermitienne et définie positive, pour résoudre
le systéme Ax = b, on pose A = D − E − F, L = D −1 E et U = D −1 F et on considère
1
le procédé itératif B(ω) x(k+1) = ( B(ω) − A) x(k) + b avec B(ω) = D ( I − ω L).
ω
76
1. Ecrire l’itération sous la forme x(k+1) = T (ω) x(k) + c (∗) et montrer que
T (ω) = ( I − ω L)−1 ((1 − ω) I + ωU ).
2. Montrer que la matrice B(ω) + B H (ω) − A est définie positive pour 0 <
ω < 2.
3. Montrer que si λ est valeur propre de Q(ω) = A−1 (2B(ω) − A) alors Reλ > 0.
4. Vérifier que Q(ω) + I est inversible et que ( Q(ω) − I )( Q(ω) + I )−1 = T (ω).
Exercice 3.9.1. Soit A = ( ai j ) une matrice carrée inversible dont les éléments di-
agonaux sont non nuls. A est écrite sous la forme A = D − L − U où D est une
matrice diagonale et L ( respectivement U) est triangulaire inférieure (respective-
ment supérieure). Pour résoudre le système Ax = b, on utilise la méthode itétative
suivante:
i −1 i −1
(k+1) (k) (k +1) (k)
aii xi = aii xi − θ ∑ ai j x j − (1 − θ )ω ∑ ai j x j
j=1 j=1
i −1 n
(k) (k)
+ (1 − ω)θ ∑ ai j x j − ω ∑ ai j x j + ω bi
j=1 j =i
1. (a) Montrer que la méthode proposée peut s’écrire sous la forme matricielle:
2. On prend θ = 1
(a) Trouver une relation entre les valeurs propres de M(θ , ω) et celles de la
matrice de Gauss-Seidel associée à A.
77
(b) Peut-on comparer la vitesse de convergence de la méthode (∗) à celle de
Gauss-Seidel ?
3. On prend U = Lt et θ = 0
(a) Trouver une relation entre les valeurs propres de M(θ , ω) et celles de la
matrice de Jacobi associée à A.
(b) En supposant que la méthode de Jacobi converge, montrer que la méthode
2
(∗) avec θ = 0 converge pour 0 < ω < où µ1 est la plus petite
1 − µ1
valeur propre de la matrice de Jacobi associée à A.
4. En supposant que A est une matrice tridiagonale, trouver une relation entre
les valeurs propres M(0, ω) et celles de M(1, ω).
( I + E + F )( x + δ x) = b
Exercice 3.9.12. Soit A = ( ai j ) une matrice carrée inversible dont les éléments di-
agonaux sont non nuls. A est écrite sous la forme A = D − L − U où D est une
matrice diagonale et L (respectivement U) est triangulaire inférieure (respective-
ment supérieure). Pour résoudreÃle système Ax ! = b, on utilise la méthode itérative
n i −1 ³ ´
(k+1) (k) (k) (k) (k+1)
suivante: aii xi = aii xi + ω bi − ∑ ai j x j + r ∑ ai j x j − x j , où r et
j=1 j=1
ω sont des réels fixés (ω non nul ) et k = 0, 1, · · ·
78
1. Montrer que la méthode proposée peut s’écrire sous la forme matricielle:
x(k+1) = M(r, ω) x(k) + c avec: M(r, ω) = ( D − rL)−1 ( aD + bL + eU )
où a , b sont des réels qu’on exprimera en fonction de r et/ou de ω.
2. Vérifier que cette méthode permet d’obtenir les méthodes de Jacobi, Gauss-
Seidel et de relaxation pour des choix appropriés de r et ω.
3. Montrer que les valeurs propres de M(r, ω) sont les racines de l’équation:
det(α D − β L − ωU ) avec α = λ + ω − 1 et β = (λ − 1)r + ω.
4. En supposant que A est une matrice tridiagonale, montrer que les valeurs
propres µ de M(0, 1) sont liées aux valeurs propres λ de la matrice générale
M(r, ω) par la relation (λ + ω − 1)2 = ωµ 2 ((λ − 1)r + ω).
79
¯ ¯ Ã√ √ !2
¡ ¯ r − x ¯2 b − a
On pourra montrer que min max ¯¯ ¯ = √
¯ √ et qu’il
r≥0 0< a≤ x≤b r + x b+ a
√ ¢
est atteint pour r = ab
Exercice 3.9.14. Soit A une matrice carrée autoadjointe. On suppose que A est
symétrique.
2. Théorème d’Ostrowski.
En supposant A définie positive et en la décomposant sous la forme
A = D − L − L>
80
b. On suppose que M + M∗ − A est définie positive.
Montrer que si x est un vecteur et y = Bx alors:
81
Chapitre 4
Solutions numériques de
l’équation f ( x) = 0
| ei +1 |
lim =c
i −→∞ | ei | p
c (constante).
(pour p = 1 on imposera la condition c < 1), c est appeleé constante de l’erreur
asymptotique.
p = 1 : convergence linéaire.
p = 2 : convergence quadratique.
p = 3 : convergence cubique.
82
4.1.2 Méthodes du point fixe ou d’approximations successives
i) Φ( I ) ⊂ I.
ii) Il existe une constante 0 < L < 1 telle que |Φ( x) − Φ( y)| ≤ L| x − y| pour tous
x, y de I.
Alors le procédé itératif xn+1 = Φ( xn ) converge vers une solution x∗ qui est la solution
unique du problème du point fixe xn+1 = Φ( xn ). Cette convergence a lieu pour n’importe
quel point de départ x0 ∈ I. De plus on a une majoration de l’erreur :
Ln
| xn+1 − xn | ≤ Ln | x1 − x0 | , | en | ≤ | x1 − x0 |.
1−L
Preuve:
On montre que la suite ( xn ) est de cauchy :
83
Soit encore µ ¶
n 1 − Lp
| xn+ p − xn | ≤ L | x1 − x0 | (4.1.1)
1−L
pour n assez grand Ln tend vers zéro et la suite xn est de Cauchy .
Soit x∗ une limite de ( xn ) alors par la continuité de Φ on a x∗ = Φ( x∗ ).
Unicité :
Supposons que ( xn ) admet deux limites x∗ et y∗ alors : x∗ = Φ( x∗ ) et y∗ = Φ( y∗ )
et par suite | x∗ − y∗ | = |Φ( x∗ ) − Φ( y∗ )| ≤ L| x∗ − y∗ | qui contredit l’hypothèse
L < 1, si on suppose que x∗ 6= y∗ .
ii) On vérifie facilement que |Φ01 ( x)| > 1 et |Φ20 ( x)| < 1.
Théorème 4.1.2. Si l’itération xn+1 = Φ( xn ) produit une suite qui converge linéairement
vers x∗ et si Φ est de classe C 1 alors
| ei +1 |
C = lim = |Φ0 ( x∗ )|.
i −→∞ | ei |
Preuve:
Il suffit d’appliquer le théorème des accroissements finis
|ei+1 | = | xi+1 − x∗ | = |Φ( xi ) − Φ( x∗ )| = |( xi − x∗ )Φ0 (ξ )| avec x∗ < ξ < xi , d’où
| ei +1 |
lim = lim |Φ0 (ξ )| = |Φ0 ( x∗ )|.
i −→∞ | ei | i −→∞
Remarque 4.1.3. Le théorème 4.1.1 donne une convergence globale c.à.d pour tout
choix x0 ∈ I. On peut avoir une condition analogue à celle du théorème 4.1.1 mais
qui n’est pas vérifiée sur I tout entier .
84
Corollaire 4.1.1 (Convergence locale ).
Si Φ est de classe C 1 dans un ouvert contenant x∗ et si |Φ0 ( x∗ )| < 1, alors :
il existe ε > 0 tel que l’itération xn+1 = Φ( xn ) converge vers x∗ si x0 satisfait :
| x0 − x∗ | ≤ ε. En d’autres termes, on obtient convergence pour des choix de x0 assez proche
de la solution x∗ .
Preuve:
Puisque Φ0 est continue dans un voisinage de x∗ et que |Φ0 ( x∗ )| < 1 alors pour
tout K tel que |Φ0 ( x)| ≤ K < 1, il existe ε > 0 tel que |Φ0 ( x)| ≤ K pour tout
x ∈ [ x∗ − ε, x∗ + ε].
Soit K fixé et ε(K ) donné alors I = [ x∗ − ε, x∗ + ε] répond aux conditions du
théorème 4.1.1
85
Le meilleur choix α̂ pour α qui optimiserait la convergence du processus itératif
(4.1.3) serait tel que xn+1 = x∗ . Or ce point est a priori inconnu, donc α̂ l’est
également.
Toutefois, on peut obtenir un estimé de α̂ de la manière suivante:
La distance entre Φ( xn ) et x∗ est égale à Φ( xn ) − x∗ = (α − 1)∆xn .
Si θ désigne l’angle de l’intersection entre la droite passant par ( xn , Φ( xn )) et
parallèle à l’axe des abscisses et la droite passant par ( xn , Φ( xn )) et ( x∗ , Φ( x∗ ))
alors
(α − 1)∆xn α−1
tan(θ ) = = (4.1.4)
α ∆xn α
Φ( x∗ ) − Φ( xn )
d’autre part tan(θ ) = .
x∗ − xn
Le théorème des accroissements finis donne:
1
α= .
1 − Φ0 (e)
Φ ( xn ) − Φ ( xn−1 ) Φ( xn ) − xn
Φ0 (e) = = .
xn − xn−1 xn − xn−1
Elle consiste à prendre h( x) = f 0 ( x). Le procédé itératif de Newton est alors donné
par :
f ( xn )
xn+1 = xn − 0 = Φ( xn )
f ( xn )
on voit aisément que f ( x∗ ) = 0 à condition que f 0 ( x∗ ) 6= 0 ce qui nous assure au
moins une convergence locale (|Φ0 ( x∗ )| < 1).
86
f ( xn ) − f ( x∗ ) 1 00
∗ )2 f (η) .
En effet xn+1 − x∗ = xn − x∗ − = ( x n − x
f 0 ( xn ) 2 f 0 ( xn )
00 00 ∗
| en+1 | 1 | f (η)| 1 | f ( x )|
donc = 0 tend vers lorsque n tend vers l’infini.
|en | 2 2 | f ( xn )| 2 | f 0 ( x∗ )|
Exemple 4.1.4. Soit à calculer la racine carrée d’un réel positif A par le schéma
itératif · ¸
1 A
xn+1 = xn +
2 xn
i) f ( a) f (b) < 0,
0
ii) f ( x) 6= 0 ∀ x ∈ [ a, b] . (Monotonie)
| f ( a)| | f (b)|
iv) 0 < b−a , 0 < b − a.
| f ( a)| | f (b)|
Alors la méthode de Newton converge vers l’unique solution x∗ de f ( x) = 0 dans [ a, b]
pour n’importe quel choix de x0 ∈ [ a, b] .
87
Preuve:
Considérons le cas f ( a) < 0 , f (b) > 0 , f 0 ( x) > 0 et f 00 ( x) ≥ 0 ( les autres cas se
traitent de façon similaire).
D’après i ) et ii ) , il existe une solution unique x∗ ∈ [ a, b] qui vérifie:
00
f ( x∗ ) − f ( xn ) 1 ∗ 2 f (ξ )
xn+1 − x∗ = xn − x∗ + = ( x n − x ) avec xn ≤ ξ ≤ x∗
f 0 ( xn ) 2 f 0 ( xn )
00
∗ f (ξ )
Par conséquent, le signe de xn+1 − x est celui de 0 .
f ( xn )
0 00
Si f > 0 et f ≥ 0 alors xn+1 ≥ x∗ ∀n, (xn ) est donc minorée par x∗ à partir de
x1 .
Par ailleurs on a deux possibilités de choix de x0 :
f ( x0 )
- ou bien x0 > x∗ et alors x1 = x0 − < x0 (puisque f ( x0 ) > f ( x∗ ) = 0)
f 0 ( x0 )
et xn > x∗ avec f ( xn ) > f ( x∗ ) = 0 ∀n ∈ N.
Dans ce cas, la suite (xn ) est décroissante minorée donc elle converge.
f ( x0 )
- ou bien x0 < x∗ et alors f ( x0 ) < 0 d’où x1 = x0 − 0 > x0 ,
f ( x0 )
f ( x1 )
mais on aussi x1 > x∗ et f ( x1 ) > 0 et par suite x2 = x1 − < x1
f 0 ( x1 )
Donc mis à part le point x0 , on aura encore une suite décroissante minorée
donc convergente.
88
Cette méthode est particulièrement interessante dans le cas où f 0 ne varie pas
trop.
Entre la méthode de Newton qui nécessite le calcul de f 0 ( xn ) à chaque étape
0
et la méthode de Newton modifiée qui utilise seulement la valeur f ( x0 ), on peut
aussi utiliser une méthode qui estime f 0 ( xϕ(n) ) pour des valeurs intermédiaires de
n.
En général, on obtient des méthodes de Newton modifiées (encore dites discrètes)
en considérant le schéma:
f ( xn )
xn+1 = xn − ,
an
où an est une approximation de f 0 ( xn ).
f ( xn ) − f ( xn−1 )
En prenant an = on obtient la méthode de la sécante (voir exercice
xn − xn−1
4.4.9).
a+b
Soit c= si f ( a) f (c) < 0 on prend b=c
2
si f (b) f (c) < 0 on prend a = c,
a−b
et on reitére le procédé, après n itérations l’intervalle [ a, b] est réduit à .
2n
89
4.1.8 Méthode ∆2 d’Aitken (accélération de la convergence)
Soit xn une suite générée par un procédé itératif, on pose ∆xn = xn+1 − xn
(∆ opérateur de différence), la méthode d’Aitken consiste à génèrer une suite yn à
partir de la suite xn par la formule:
(∆xn )2 ( xn+1 − xn )2
yn = xn − = x n − .
∆2 xn xn+2 − 2xn+1 + xn
On a alors le résultat suivant:
Théorème 4.1.4. Supposons qu’il existe un réel k tel que |k| < 1 pour lequel la suite xn
vérifie:
xn 6= x∗ et xn+1 − x∗ = (k + δn )( xn − x∗ ) avec lim δn = 0.
n−→∞
alors la suite yn est bien définie pour n assez grand et on a
yn − y∗
lim = 0.
n−→∞ xn − x∗
Autrement dit, la suite yn converge plus rapidement que la suite xn vers la limite commune
x∗ .
| en+1 |
Remarque 4.1.7. On a lim = lim (k + δn ) = k, donc il y’a convergence
n−→∞ |en | n−→∞
linéaire avec k < 1, encore dite géometrique.
Preuve:
Par conséquent, pour n assez grand, ∆2 xn est différent de 0 car k < 1, µn tend vers
0 et par suite, yn est bien définie.
Par ailleurs on a
90
d’où
(∆xn )2 ((k − 1)2 + σn2 )e2n
yn − x∗ = xn − x∗ − = e n −
∆2 xn ((k − 1)2 + µn2 )en
ou encore µ ¶
∗ ((k − 1)2 + σn2 )
yn − x = 1 − en
((k − 1)2 + µn2 )
i.e
yn − x∗ ((k − 1)2 + σn2 )
= 1 −
xn − x∗ ((k − 1)2 + µn2 )
yn − x∗
et par passage à la limite on obtient le résultat annoncé lim = 0.
n−→∞ xn − x∗
Posons yn = ϕ( xn )
zn = ϕ( yn ) = ϕoϕ( xn ).
( yn − xn )2
xn+1 = xn − ,
zn − 2yn + xn
ceci conduit à un procédé itératif avec la nouvelle fonction Ψ, soit Ψ( xn ) = xn+1
avec Ψ définie par
xϕoϕ( x) − (ϕ( x))2
Ψ( x) = .
ϕoϕ( x) − 2ϕ( x) + x
On obtient alors les deux résultats:
Théorème 4.1.5.
i) Si Ψ( x∗ ) = x∗ alors ϕ( x∗ ) = x∗ ,
Théorème 4.1.6.
Soit ϕ une fonction de classe C p+1 définissant un procédé itératif xn+1 = ϕ( xn ) dont la
convergence est d’ordre p. Alors la fonction Ψ donne un procédé itératif dont la convergence
est d’ordre 2p − 1 si p > 1 et d’ordre au moins égal à 2 si p = 1 et ϕ0 ( x∗ ) 6= 1.
91
4.2 Cas des polynômes: Solutions numériques des équations
algébriques
On s’intéresse à l’équation P( x) = 0, avec
P( x) = a0 xn + a1 xn−1 + · · · + an , ai ∈ R et a0 6= 0. (4.2.1)
P( x) = (· · · ((( a0 x + a1 ) x + a2 ) x + a3 ) · · · ) x + an ) (4.2.3)
P( x) = ( x − r1 ) Q( x) (4.2.4)
92
où Q( x) est un polynôme de degré n − 1. Ainsi théoriquement, une fois obtenue
une première racine, on peut commencer la recherche d’une autre racine par un
polynôme de degré strictement inférieur. Successivement, on poursuit cette procéd-
ure jusqu’à l’obtention de l’ensemble des n racines du polynôme P( x). Rappelons
que les polynômes à coefficients complexes se factorisent en un produit de monômes
de degré 1. Cette propriété exprime le fait que les polynômes à coefficients com-
plexes ont l’ensemble de leurs racines dans le plan complexe.
n
P( x) = ∏ ( x − ri ) (4.2.5)
i =1
n
ln(| P( x)|) = ∑ ln(|x − ri |) (4.2.6)
i =1
n
d ln(| P( x)|) 1
dx
= ∑ x − ri
i =1
P0 ( x)
= (4.2.7)
P( x)
= G (4.2.8)
n
d2 ln(| P( x)|)) 1
−
dx2
= ∑ ( x − ri )2
i =1
µ ¶2
P0 ( x) P00 ( x)
= −
P( x) P( x)
= H (4.2.9)
x − r1 = a (4.2.10)
x − ri = b, i = 2, · · · , n (4.2.11)
93
en insérant les équations (4.2.10) et (4.2.11) dans les équations (4.2.7) et (4.2.9), on
en déduit respectivement les relations suivantes
1 n−1
+ = G
a b
1 n−1
+ 2 = H (4.2.12)
a2 b
Après éliminition de b, la valeur de a est
n
a= p (4.2.13)
G± (n − 1)(nH − G2 )
Le signe placé devant la racine du dénominateur est choisi tel que le dénominateur
soit le plus grand possible, x − a devient alors la nouvelle valeur de départ et on
itére le processus.
P ( x ) = xn + a1 xn−1 + · · · + an−1 x + an
P( x) = ( x2 + bx + c) Q( x) + rx + s
94
4.2.4 Méthode de Maehly
P ( x ) = xn + a1 xn−1 + · · · + an−1 x + an
P( xk )
xk+1 = xk −
P0 ( xk )
P( x)
Soit r1 la première racine trouvée et P1 ( x) = , on a
x − r1
P0 ( x) P( x)
P10 ( x) = −
x − r1 ( x − r1 )2
P( xk )
xk+1 = xk −
P( xk )
P0 ( xk ) −
xk − r1
P( x)
Pj ( x) =
( x − r1 ) · · · ( x − r j )
alors
j
P0 ( x) P( x) 1
Pj0 ( x)
( x − r1 ) · · · ( x − r j ) ( x − r1 ) · · · ( x − r j ) i∑
= −
=1 x − ri
P( xk )
xk +1 = xk −
j
1
P0 ( xk ) − P( xk ) ∑
i =1
x k − ri
A
|α | < R avec R = 1 + (4.2.14)
| a0 |
95
Preuve :
Soit α une racine de P.
Si |α | ≤ 1 alors α < R.
Si |α | > 1 la formule (4.2.1) entraı̂ne :
| P(α )| = | a0α n + · · · + an |
≥ | a0 ||α |n − | a1 ||α |n−1 − · · · − | an |
n
≥ | a0 ||α |n − A ∑ |α |k
k=1
µ ¶
n |α |n − 1
≥ | a0 ||α | − A
|α | − 1
· ¸
A
≥ | a0 | − |α |n
(|α | − 1)
A
Il s’ensuit que | P(α )| > 0 si | a0 | − ≥ 0, c’est à dire | P(α )| > 0 si |α | ≥ R.
(|α | − 1)
Donc si |α | ≥ R on obtient P(α ) 6= 0 ce qui contredit notre hypothèse. Par
conséquent, toute racine α de (4.2.1) vérifie (4.2.14).
r < |α | < R.
96
Théorème 4.2.2 (Théorème de Lagrange).
Si a0 > 0 et ak est le premier des coefficients négatifs du polynôme P( x)r
, on peut prendre
B
pour limite supérieure des racines réelles positives le nombre R = 1 + k où B est la
a0
plus grande des valeurs absolues des coefficients négatifs de P( x).
Preuve :
Elle est analogue à celle du théorème 4.2.1. En effet, soit α une racine de P, sup-
posée réelle positive.
Si α < 1 on a bien α < R.
Si α > 1, en écrivant l’équation (4.2.1) dans laquelle on remplace tous les coeffi-
cients positifs a1 , · · · , ak−1 par zéro et tous les coefficients ak , ak+1 , · · · , an par − B,
on obtient la majoration :
α n−k +1
P(α ) ≥ a0α n − B(α n−k + · · · + α 0 ) = a0α n − B .
α−1
α n−k +1 h i
P(α ) > a0 (α − 1)k − B
α−1
r
B
et par suite, pour α ≥ R = 1 + k , on a P(α ) > 0. Donc toute racine positive α
a0
de (4.2.1) vérifie α < R.
Preuve:
En fait, on peut construire une suite croissante 0 < c1 ≤ c2 ≤ · · · ≤ cn qui vérifie
( x − c )n (n)
P( x) = P(c) + ( x − c) P0 (c) + · · · + P (c)
n!
et par suite, pour toute racine positive α on a α ≤ c , sinon on aboutit à une contra-
diction avec P(α ) > 0.
97
Théorème 4.2.4 (Théorème de Cauchy ).
Soient
Alors Q1 et Q2 admettent chacun une seule racine positive, notons les R et r, alors toute
racine α de P( x) vérifie r < |α | < R.
Corollaire 4.2.4. Si ∆N = 1 alors une et une seule racine réelle dans l’intervalle [ a, b].
98
Preuve :
On applique le théorème 4.2.5 à P(k) sur [ a, b] = [0, ∞] en remarquant que
an−k = 0, k = 0, 1, · · · , n. Les an−k ont même signe que P(k) (0) d’où les change-
ments
³ de signe dans ( a0 , ´ a1 , · · · , an ) Sont les mêmes que ceux de la suite
0 (
P(0), P (0), · · · , P (0) ils sont égaux à N (0) . Par ailleurs, P(k) (+∞) ont toutes
k )
Corollaire 4.2.6. Le nombre de racines réelles négatives est égal au nombre de changements
de signes −2k.
Exemple 4.2.2. P( x) = x4 − x3 − x2 + x − 1
Ici N = N = N = 3. Le nombre de racines positives est égal à 3 − 2k soit 3 ou 1.
Posons Q( x) = P(− x) = x4 + x3 − x2 − x − 1.
Le nombre de racines négatives est égal à 1 − 2k d’où une seule possibilité : 1.
En résumé P( x) admet soit 3 racines positives et une négative soit une positive, une
négative et 2 racines complexes conjuguées.
99
b- Relation de Newton
p p p
Si S p désigne la somme S p = x1 + x2 + · · · + xn = − p on a
a0 S1 + a1 = 0
a0 S2 + a1 S1 + 2a2 = 0
. ..
..
.
a0 Sn + a1 Sn−1 + · · · + an−1 S1 + nan =0
c- Méthode de Graeffe :
On suppose ici que les racines sont nettement séparées c’est à dire :
Définition 4.2.2. soit f 0 , f 1 , · · · , f n une suite de fonctions définies sur [ a, b], sup-
posées continues (avec f 0 de classe C 1 ). On dit que cette suite de fonctions est une
suite de Sturm si
100
1. f 0 n’admet que des racines simples dans [ a, b].
Exemple 4.2.4. Soit P un polynôme de degré n qui n’admet que des racines simples
et réelles α1 , α2 , · · · , αn , on obtient une suite de Sturm de la façon suivante :
f 0 ( x) = P( x)
f 1 ( x) = P0 ( x)
..
.
f n−2 ( x ) = f n−1 ( x ) qn−1 ( x ) − f n ( x ) avec d◦ f n = 0
• f n = constante.
Preuve :
Pour tout x ∈ [ a, b], posons N ( x) = C ( a) − C ( x), N ( x) est un entier qui vérifie :
101
signe en α )
α −ε α α +ε α −ε α α +ε α −ε α α +ε α −ε α α +ε
f i −1 − − − + + + − − − + + +
fi − 0 − − 0 − + 0 + + 0 +
f i +1 + + + − − − + + + − − −
chang. de signes 1 1 1 1 1 1 1 1 1 1 1 1
α −ε α α +ε α −ε α α +ε
f0 − 0 + + 0 −
f1 + + + − − −
chang. de signes 1 . 0 1 . 0
Dans les deux cas, on remarque que la suite ( f 0 (α − ε), f 1 (α − ε), · · · , f n (α − ε))
présente un changement de signe de plus que la suite ( f 0 (α + ε), f 1 (α + ε), · · · , f n (α + ε))
dans Iε . D’où N (α + ε) = N (α − ε) + 1 et N (b) = C ( a) − C (b) = n.
La fonction N ( x) se présente comme une fonction en escalier.
102
4.3 Applications
Application 4.3.1. Dynamique élémentaire de population
Considérons une population de taille P(t) à un certain temps t et supposons
que nous nous interessions à son évolution dans le temps. Quatre paramètres peu-
vent contribuer au changement de P(t) :
• B: naissances( birth),
• D: mortalités( death),
• E: émigration,
• I: immigration.
Pn+1 − Pn = R (4.3.1)
où R = B + I − D − E.
On voit facilement que la population croı̂t ou décroı̂t avec le temps selon que R est
positif ou négatif.
Pn+1 = (1 + r) Pn (4.3.2)
103
une croissance géométrique de la population + une croissance arithmétique des moyens de
subsistance = grande misère humaine.
Remarque 4.3.1. La théorie de Malthus a été beaucoup critiquée. Elle est notam-
ment mise en cause par les économies des pays industrialisés mais aussi et surtout
par la répartition inégale et injuste des richesses dans d’autres pays.
Pn+1 = g( Pn )
on voit aisement que l’équation logistique admet deux points d’équilibre (les points
fixes de g) qui sont:
l’équilibre trivial P∗ = 0 et l’équilibre non trivial atteint pour la capacité d’environn-
ement P∗∗ = K.
En étudiant la dérivée de g aux points d’équilibre, on déduit que si r > 0 alors le
point d’équilibre trivial est instable (g0 (0) > 1) tandis que pour le point d’équilbre
non trivial ( g0 (K ) = r − 1), la zone de stabilité est donnée par 0 < r < 2.
L’équation logistique constitue un exemple pédagogique interessant. La varia-
tion de r exhibe différentes situations allant de l’équilibre stable (0 < r < 1) voir
figure 4.1, au chaos (r = 3) voir figure 4.2, en passant par des oscillations cycliques
(2 < r < 3) voir figure 4.2.
104
Convergence
30
25
20
Population
15
10
0
0 5 10 15 20 25 30 35 40 45 50
Temps
Cycle Chaos
30 30
25 25
20 20
Population
Population
15 15
10 10
5 5
0 0
0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50
Temps Temps
105
Application 4.3.4. Equation logistique et captivité (harvest)(Schaefer, 1954)
Si on considère une population de poissons vivant dans un environnement aux
ressources limitées et assujettie à une captivité (harvest) on obtient l’équation
r 2
Pn+1 = (1 + r) Pn − P − qEPn (4.3.4)
K n
L’équation (4.3.4) admet deux points d’équilibre qui sont l’équilibre trivial P∗ = 0,
qE
et l’équilibre non trivial P∗∗ = K (1 − ).
r
L’équilibre trivial correspond à l’épuisement du stock de poissons et l’équilibre
non trivial permet d’élaborer des stratégies optimales en fonction du stock Pn , de
l’effort de pêche E et de la constante q mesurant la captivité.
106
4.4 Complément bibliographique
• Population de rats et suites de Fibonacci (Fibonacci, 1202)
Fibonacci a proposé un modèle de croissance d’une population de rats en
partant d’un couple de rats (une femelle et un male) qui atteignent maturité
pour donner un nouveau couple, qui, après maturité, donne un autre couple
qui s’ajoute au couple que donne le premier couple s’il survit, etc... Les suites
de Fibonacci vérifient
B0 = 1, B1 = 1,
Bn+1 = Bn + Bn−1 , n = 1, 2, · · ·
Bn = hn + λ1 b1 Bn−1 + λ2 b2 Bn−2 + · · · + λn bn B0 .
Cette équation peut être vue comme une généralisation des suites de Fibonacci
(en prenant hn = 0, λ1 = λ2 = 1 et λ j = 0 pour j ≥ 3).
107
a- Equation de Ricker(1954): Bn+1 = exp [r(1 − Bn /K )] Bn .
R0
b- Equation de Beverton-Holt(1957): Bn+1 = Bn .
1 + [( R0 − 1)/K ] Bn
c- Equation généralisée de Beverton-Holt (1957):
Rβ0
Bn+1 = Bn .
{1 + [( R0 − 1)/K ] Bn }β
108
4.5 Exercices
Exercice 4.5.1. Etudier la convergence de l’itération définie par
xs + α sx
xn+1 = Φ ( xn ) avec Φ( x) = .
sxs−1 + α
Exercice 4.5.2. Soit f une fonction réelle admettant un minimum m en un point x∗
de l’intervalle [ a, b].
On cherche à encadrer ce minimum par une méthode de dichotomie. La méthode
du nombre d’or consiste à prendre L1 = b − a puis à localiser x∗ dans un intervalle
de longueur Ln en calculant Li par la relation de récurrence: Li = Li−2 − Li−1 et en
lui imposant que le rapport de réduction soit constant c’est à dire
Li L
= i−1 = constante.
Li −1 Li −2
a) Trouver le nombre d’or α .
( xn+1 − xn )2
zn = xn +
xn+2 − 2xn+1 + xn
zn − x∗
1. Chercher lim et comparer la vitesse de convergence des suites ( xn )
n−→∞ xn − x∗
et ( zn ).
Exercice 4.5.4. Soit f une fonction réelle de classe C 2 sur un intervalle [ a, b] et telle
que f 00 ( x) = m ∀ x ∈ [ a, b] .
109
1. Donner l’ordre de convergence de la méthode de fausse position avec
xn − x0
x0 fixé (x0 = a), xn+1 = xn − f ( xn ), quelle est la constante
f ( xn ) − f ( x0 )
d’erreur asymptotique.
Pn+1 = Pn {1 + r(1 − Pn /K )}
Pn+1 − Pn = f ( Pn ) − g( Pn )
avec
r 2 β P2
f ( Pn ) = rPn − Pn , g( Pn ) = 2 n 2 .
K α + Pn
En traçant les graphes de f et g en fonction de Pn , montrer que le système admet
un point non trivial d’équilibre stable ou 3 points d’équilibre dont 2 stables et un
instable.
dP(t) r
= rP(t) − P2 (t) − E
dt K
1
a) Montrer que si E < rK; il ya deux possibilités d’équilibre, l’un stable, l’autre
4
instable.
1
b) Montrer que si E > rK alors la population disparaitra en un temps fini.
4
110
c) Que deviennent les résultats a) et b) si E est remplacé par EP.
1. Montrer que cette équation admet une racine dans l’intervalle [1, 2]
2. Montrer que: l’itération xn+1 = x3n − x2n − 1 ne converge pour aucune valeur
1 1
initiale x0 ∈ [1, 2] mais que la méthode xn+1 = 1 + + 2 converge pour
· ¸ x n x n
7
toute valeur initiale x0 ∈ ,2 .
4
Exercice 4.5.9. On suppose que l’équation x = g( x) admet une seule racine α et que
g est de classe C p dans l’intervalle I = { x tels que: | x − α | ≤ r avec r > 0} . Pour
chercher α , on utilise le procédé itératif : xn+1 = g( xn ).
111
3. Etant donné deux réels a et b tels que f ( a) f (b) < 0, de quelle façon peut-on
adapter la méthode de la sécante pour garantir la convergence vers x∗ ∈ [ a, b].
f ( x) = x + λ (Φ( x) − x).
xn+1 = f ( xn ) .
3. Soit à résoudre f ( z) = 0, z ∈ C.
En écrivant z = x + iy et f ( z) = u( x, y) + iv( x, y) et sachant que:
∂u ∂v ∂u ∂v
= et =− .
∂x ∂y ∂y ∂x
P ( λ ) = λ n + a1 λ n−1 + · · · + an ,
112
..
.
¯ ¯
¯ a1 a3 a5 ··· ··· ··· ¯
¯ ¯
¯ ¯
¯ 1 a2 a4 ··· ··· ··· ¯
¯ ¯
¯ 0 a1 a3 ··· ··· ··· ¯
¯ ¯
Dk = ¯ ¯ > 0 k = 1, 2, · · · , n
¯ 0 1 a2 ··· ··· ··· ¯
¯ ¯
¯ ··· ··· ··· ··· ··· ··· ¯
¯ ¯
¯ ¯
¯ 0 0 ··· ··· ··· ak ¯
α1 ≥ α2 ≥ · · · ≥ αn .
4. Conclure.
113
2. On suppose que P est un polynôme de degré 1 vérifiant la condition C1 , quel
procédé d’itération obtient-on en evaluant P( y) et en prenant xn+1 = P(0)?
P( x) = 2x5 − 100x2 + 2x − 1.
En posant P0 = 1 , P1 = α1 − λ et Pk = det( Ak − λ Ik ).
2. Montrer que les zéros de Pn+1 sont distincts et strictement séparés par les
zéros de Pn .
114
3. En déduire que ( Pn , Pn−1 , · · · , P0 ) est une suite de Sturm.
1 2 0 0
2 3 −1 0
4. Application: soit A la matrice donnée par 0 −1 4
1
0 0 1 1
115
Chapitre 5
Méthodes numériques de
résolution des systèmes non
linéaires
i) Φ : M −→ M continue.
Preuve :
Analogue à celle du cas scalaire en remplaçant |.| par ||.||.
116
Théorème 5.1.2. On suppose ¯ que Φ ¯possède un point fixe α et que les composantes Φi
¯ ∂Φi ( x) ¯
sont de classe C 1 et vérifient : ¯¯ ¯ ≤ λ avec λ < 1.
∂xi ¯ n ³ ´
Alors pour tout x(0) ∈ V = { x/|| x − α || ≤ d}, la suite x(k) définie par
³ ´
x(k) = Φ x(k−1) vérifie x(k) ∈ V et x(k) converge vers α .
Preuve :
n
∂Φi (ξ (t))
Soient x et y ∈ V, Φi ( x) − Φi ( y) = ∑ ∂xi
(x j − y j )
j=1
où ξ (t) = (1 − t) y + tx, t ∈]0, 1[, on a ξ (t) − α = (1 − t)( y − α ) + t( x − α ) et
kξ (t) − α k ≤ (1 − t)d + td = d donc ξ (t) ∈ V.
Par suite
n ¯
¯ ¯
¯ ∂Φi (ξ (t)) ¯¯
|Φi ( x) − Φi ( y)| ≤ ∑ ¯ ¯ k x − yk∞
j=1
∂xi
λ
≤ n k x − yk∞ = λ k x − yk∞
n
117
³ ´ ³ ´ ³ ´
et F x(0) = (6, 2)> , on résoud donc : JF x (0 ) C (0) = − F x ( 0)
ce qui donne : C(0) = (−3/2, −1/2)> , x(1) = x(0) + C (0) = (−1/2, 5/2)>
Posons Φ( x) = x − JF−1 ( x) F ( x)
Si x ∈ B( x∗ , d) donc
° °
° °
kΦ( x) − x∗ k = ° x − JF−1 ( x) F ( x) − x∗ °
° °
° °
= ° x − x∗ − JF−1 ( x) [ F ( x) − F ( x∗ )]°
° °
° °
= ° JF−1 ( x) [ F ( x∗ ) − F ( x) − JF ( x)( x∗ − x)]°
° ° °Z 1 °
° −1 ° ° °
≤ ° JF ( x)° ° ( JF (ξ (t)) − JF ( x))( x − x)dt°
° ∗
°
0
° °Z 1
° °
≤ ° JF−1 ( x)° k JF (ξ (t)) − JF ( x)k k x − x∗ k dt
0
Z 1
≤ αγ k x − x∗ k kξ (t) − xk dt
0
Z 1
∗
≤ αγ k x − x k tdt k x − x∗ k
0
αγ
≤ k x − x∗ k2
2
αγ
≤ d k x − x∗ k
2
d’où kΦ( x) − x∗ k < k x − x∗ k ≤ d si la condition iii ) est vérifiée et comme on a
° ° αγ ° °2 ³ ´
° (k +1) ° ° (k) °
°x − x∗ ° ≤ ° x − x∗ ° alors la suite x(k) converge vers x∗ et la conver-
2
gence est quadratique.
118
Remarques 5.1.1. 1. D’après le théorème 5.1.3, la méthode de Newton converge
avec ordre de convergence quadratique mais à condition que le choix de
départ soit assez proche de x∗ , x(0) ∈ B( x∗ , d). Il s’agit donc d’une conver-
gence locale. Comme dans le cas scalaire, on peut obtenir une convergence
globale mais au prix d’hypothèses supplémentaires sur F.
on obtient ainsi une méthode de Newton modifiée qui est évidemment moins
précise que la précédente mais moins coûteuse.
B( x0 , r) = { x/k x − x0 k < r} ⊂ C0 ,
119
αβγ
h= < 1,
2
α
r= ,
1−h
si F ( x) vérfie les conditions suivantes
° °
i) JF ( x) est inversible et ° JF−1 ( x)° ≤ β ∀ x ∈ C0
ii) k JF ( x) − JF ( y)k ≤ γ k x − yk ∀ x, y ∈ C0 .
° °
iii) ° JF−1 ( x0 ) F ( x0 )° ≤ α .
ii) k JF ( x) − JF ( y)k ≤ γ k x − yk ∀ x, y ∈ C0 .
° °
iii) ° JF−1 ( x0 ) F ( x0 )° ≤ α ,
h = αβγ ,
√
1 ± 1 − 2h
r12 = α.
h
1
Si h ≤ et B( x0 , r12 ) ⊂ C0 alors la méthode de Newton converge vers la solution unique
2
x∗ .
120
Ces deux problèmes sont liés de la manière suivante
A- Méthodes de descente :
Nous cherchons des solutions du problème 2 en utilisant des méthodes itératives
donnant une suite de points x(1) , x(2) , · · · à partir d’un point initial x(0) et qui con-
verge vers un³minimum ∗
´ ³local´x . On voudrait
³ ´obtenir ³ une convergence
´ monotone
c’est à dire f x ( 0 ) > f x ( 1 ) > ··· > f x ( i ) > f x ( i + 1 ) > · · · > f ( x∗ ).
Pour cela, on a besoin de deux paramètres pour passer d’un itéré x(i) à l’itéré suiv-
ant x(i+1) .
Nous
³ cherchons ´ une³direction
´ de manière à aller en descendant, en d’autres termes:
f x(i) + ερ(i) < f x(i) pour tout ε assez petit.
En appliquant la formule de Taylor à f ( x(i) + ερ(i) ) on obtient :
³ ´ ³ ´ ³ ´>
f x(i) + ερ(i) = f x(i) + ε ∇ f ( x(i) ) ρ(i) + O(ε2 )
121
³ ´>
et par suite nous exigeons que ∇ f ( x(i) ) ρ(i) < 0 (∗)
En posant g( x) = ∇ f ( x) pour alléger les notations, l’inégalité (∗) s’écrit :
³ ´>
g (i ) ρ(i) < 0,
c’est donc la condition qui permet de choisir à chaque étape i la direction ρ(i) .
Une fois la direction ρ(i) choisie, nous sommes ramenés à un problème de min-
imisation à une dimension de f dans la direction ρ(i) . Le problème se pose donc
comme suit : ³ ´ ³ ´
Trouver α ∗ tel que f x(i) + α ∗ ρ(i) < f x(i) + αρ(i) pour tout α dans un voisi-
nage de α ∗ .
Ce procédé s’appelle une recherche en ligne. On peut utiliser des méthodes numéri-
ques pour encadrer α ∗ . ³ ´
Si on utilise une recherche en ligne exacte de telle sorte que f x(i) + α ∗ ρ(i) Soit
un minimum local de f dans la direction ρ(i) , alors on doit avoir
∂f
= 0 pour α = α ∗ ,
∂α
µ ¶>
∂f
ou encore ρi = 0 au point x = x(i) + α ∗ ρ(i) .
∂x
³ ´>
Ce qui donne la condition d’orthogonalité suivante : g(i+1) ρi = 0. En d’autres
termes, le gradient de f au point x(i+1) doit être orthogonal à la direction de recherche
si on utilise la recherche en ligne exacte.
Choix des directions:
Exemple 5.2.1. f ( x1 , x2 ) = 2 2
³ 2x1´− 6x1 x2 + 5x2 − 2x1 + 2x2 .
Prenons x(0) = (0, 0)> , f x(0) = 0,
µ ¶³ ´
(0) (1) ∂f (1)
1. On garde x2 = 0 et en cherche x1 tel que x1 , 0 = 0 on ob-
∂x1
tient
122
∂ f ³ (1) ´ (1) 1
x1 , 0 = 4x1 − 6x2 − 2 d’où x1 = et x(1) = (1/2, 0)> , f ( x(2) ) =
∂x1 2
−1
2
µ ¶
(1) 1 (2) ∂ f ³ (1) (2) ´
2. On garde x1 = et en cherche x2 tel que x1 , x2 =0
2 ∂x2
∂ f ³ (1) (2) ´ (2) 1
Ce qui donne x1 , x2 = −6x1 + 10x2 + 2 d’où x2 = d’où
∂x2 10
³ ´ −11
x(2) = (1/2, 1/10)> et f x(2) =
20
µ ¶
(2) 1 (3) ∂ f ³ (3) (2) ´
3. On garde x2 = et on cherche x1 tel que x1 , x2 = 0···
10 ∂x1
b- Méthode de la meilleure descente (de la plus grande pente).
³ ´ ³ ´ ³ ´>
En écrivant f x(i) + ρ = f x(i) + g(i) ρ + · · ·
³ ´>
On choisit ρ de telle sorte à minimiser g(i) ρ avec une condition de min-
imisation ||ρ|| = C pour une certaine norme ||.||. On obtient alors (avec la
norme euclidienne)
³ ´> ° ° ³ ´> ° °
° ° ° °
− g(i) ρ < ° g(i) ° ||ρ||2 pour tout ρ, d’où g(i) ρ > − ° g(i) ° kρk2 pour
2
C ³ ´> ° ° ³2 ´>
( ) ( ) ° °
tout ρ. Soit ρe = − ° i
° g alors g i e = −C ° g ° , donc g(i) ρ >
ρ ( i )
° (i ) ° 2
°g °
° ° 2
³ ´>
° °
−C ° g(i) ° = g(i) ρe ∀ ρ.
2
Dans la méthode classique on prend ρ(i) = − g(i) .
123
³ ´>
En effet ∇ f = JF> F et on vérifie aussi que g(k) ρ(k) < 0 puisque
³ ´> ¡ ¢> ¡ −1 ¢
g ( k ) ρ(k) = Jk> Fk − Jk Fk = − Fk> Jk Jk−1 Fk = − Fk> Fk .
B- Méthodes des directions conjuguées
Ces méthodes sont aussi basées sur les algorithmes de descente :
x(k+1) = x(k) + αk ρ(k) mais elles choisissent des directions conjugueés.
On sait que pour les méthodes de descente avec direction en ligne exacte on
obtient : ³ ´>
g(k+1) ρ(k) = 0 k = 0, 1, · · ·
posons s(k) = x(k+1) − x(k) et t(k) = g(k+1) − g(k) et supposons que f est une
fonction quadratique f ( x) = 1/2x> Gx + b> x + c de telle sorte que
g( x) = ∇ f ( x) = Gx + b et ∇2 f ( x) = G Hessien . Il s’ensuit alors que :
³ ´
t(k) = G s(k) = αk Gρ(k)
en particulier, on obtient :
³ ´>
g(n+1) ρ( j) = 0 j = 1, 2, · · · , N
( j) n
³ ρ sont
et comme les ´ linéairement indépendantes et donc engendrent R , on con-
clue que g x(n+1) = 0, nous venons donc de prouver le théorème suivant
Théorème 5.2.1. Si f est une fonction quadratique avec Hessien définie positive alors les
directions de recherche sont conjuguées par rapport à G.
124
L’algorithme est donc le suivant:
(i + 1 ) = x (i ) + α ρ (i )
x i
ρ (i ) = − g (i ) + β i ρ (i − 1 )
(0)
ρ = − g(0)
On démontre alors que pour les fonctions quadratiques les αi et les βi sont données
³ ´> ³ ´>
g (i ) g (i ) g (i ) g (i )
par αi = ¡ ¢> , où G est le Hessien et βi = ¡ ¢> .
ρ(i) Gρ(i) g (i − 1 ) g (i − 1 )
Remarques 5.2.2. 1. La méthode du gradient conjugué appliquée à la résolution
du système Ax = b se présente sous une forme similaire .
• on choisit x(0) ∈ Rn
• on prend ρ(0) = r(0) = b − Ax(0)
• pour k = 1, 2, · · ·
– si ρ(k) = 0 stop, x(k) est la solution de Ax = b
³ ´>
r(k) r(k)
– sinon calculer ak = ¡ ¢> , x(k+1) = x(k) + ak ρ (k) ,
ρ ( k ) Aρ ( k )
³ ´>
r(k +1) r(k +1)
r(k+1) = r(k) − ak Aρ(k) , bk = ³ ´> et
r(k) r(k)
ρ (k +1) = r(k +1) + bk ρ (k) .
125
5.3 Applications
Une poupulation ne peut vivre isolée dans un environnement, elle entre en contact
avec d’autres populations d’une manière ou d’une autre, et la nature de ces inter-
actions contribue à son évolution. Un ensemble de deux ou plusieurs espèces est
dit communauté. Une classification traditionnelle d’un système de populations est
basée sur la nature des interactions interespèces et le comportement intraespèce.
Une grande et importante classe de modèles écologiques consacrés à la modéli-
sation de la dynamique des communautés de populations, notamment avec inter-
action, peut être obtenue à partir du modèle général de Lotka-Volterra
n
dNi
= Ni (ri − ∑ ai j N j ) i = 1, · · · , n
dt j=1
126
+ + Coopération ou mutualisme ou symbiose
+ - Proie-Prédateur ou ressource-consommateur ou hôte-parasite
+ 0 Commensalisme
- - compétition
- 0 Amensalisme
0 0 Neutralisme
interessé trouvera dans la littérature les études complètes des comportements des
différents systèmes avec les outils mathématiques et écologiques nécessaires.
-Cas continu : (
x0 = r1 x − b1 xy
y0 = −r2 x + b2 xy
x: densité de populations des proies.
y: densité de populations des prédateurs.
r1 : taux de croissance des proies.
r2 : taux de mortalité des prédateurs.
b1 et b2 : taux de prédation.
∆x x(t + ∆t) − x(t)
L’approximation de = avec ∆t = 1 donne:
∆t ∆t
(
xn+1 = ( 1 + r1 − b1 yn ) xn
yn+1 = ( 1 − r2 + b2 xn ) yn
En l’absence des prédateurs les proies croı̂ent infiniment mais en réalité cette crois-
sance est soumise à la capacité d’accueil de l’environnement, on introduit alors un
terme −cx2n qui caractérise la compétition intraespèce au sein de la communauté
des proies, les nouvelles équations sont :
(
xn+1 = (1 + r1 − cxn − b1 yn ) xn
yn+1 = b2 xn yn
127
Exemple 5.3.1. ([41])
0.45
0.3
0.4
0.35
0.25
Population
Population
0.3
0.2
0.25
0.2
0.15
0.15
0.1 0.1
0 10 20 30 40 50 60 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Temps
Figure 5.1: Convergence Temps
PREDATION PREDATION
0.5 0.3
Proie
Predateur
0.45 0.28
0.26
0.4
0.24
0.35
Population
Population
0.22
0.3
0.2
0.25
0.18
0.2
0.16
0.15
0.14
0.1 0.12
0 10 20 30 40 50 60 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Temps Temps
128
PREDATION
0.35
0.3
0.25
0.2
population 0.15
0.1
0.05
Proie
Predateur
0
0 10 20 30 40 50 60
Temps
0.55
0.5
0.5
0.45
0.4
Population
Population
0.4
0.3
0.35
0.2 0.3
0.25
0.1
0.2
0
5 10 15 20 25 30 35 40 45 50 55 60 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65
Temps
Figure 5.4: 2 cycles Temps
PREDATION PREDATION
0.8 0.8
Proie
Predateur
0.7
0.7
0.6
0.6
0.5
0.5
Population
Population
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0 0.1
5 10 15 20 25 30 35 40 45 50 55 60 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Temps Temps
Figure 5.5: 4 cycles
129
PREDATION PREDATION
0.8 0.7
Proie
Predateur
0.7
0.6
0.6
0.5
0.5
0.4
Population
Population
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0 0
5 10 15 20 25 30 35 40 45 50 55 60 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Temps Temps
L’interaction entre les deux espèces a un effet négatif sur les deux densités; le
modèle est régi par le système d’équations suivantes:
(
xn+1 = (1 + r1 − a11 xn − a12 yn ) xn
yn+1 = (1 + r2 − a22 yn − a21 xn ) yn
r1 = 0.9 , a11 = 0.5 , a12 = 0.4 , a21 = 0.3 , r2 = 0.8 , a22 = 0.6 =⇒ coexistence des
deux espèces voir figure(5.8);
130
COMPETITION COMPETITION
0.6 0.8
0.7
0.5
0.6
0.4
0.5
0.3
Population
Population
0.4
0.2
0.3
0.1
0.2
0
0.1
Population1 Population1
Population2 Population2
−0.1 0
0 10 20 30 40 50 60 5 10 15 20 25 30 35 40 45 50 55 60
Temps Temps
COMPETITION
1.6
Population1
Population2
1.4
1.2
1
population
0.8
0.6
0.4
0.2
0 10 20 30 40 50 60
Temps
131
5.3.5 Modèle à coopération
COOPERATION
0.8
Population1
Population2
0.7
0.6
0.5
Population
0.4
0.3
0.2
0.1
0
5 10 15 20 25 30 35 40 45 50 55 60
Temps
N 0 = rN − cNP
P0 = bNP − mP.
Et depuis lors, ce modèle avec ses variantes et ramifications n’arrête pas de sus-
citer l’intérêt des chercheurs. A titre indicatif, nous citons les travaux de Scuda &
Ziegler(1978) et Oliveira & Conolly(1982) qui donnent des traductions et des com-
mentaires des principaux écrits de Volterra.
132
Pour les modèles à compétition, la formulation classique due à Lotka(1932)
et Volterra(1926) suppose une compétition par interference dans laquelle chaque
espèce empêche la croissance de l’autre. Le modèle se présente comme suit:
N10 = r1 N1 (1 − N1 /K )
N20 = r2 N2 (1 − N2 /K )
May & Leonard(1975) ont considéré un modèle de compétition non transitive en-
tre trois espèces, d’autres auteurs ont étudié la possibilité de coexistence de deux
prédateurs en compétition par rapport à une troisième espèce. Le lecteur intéressé
pourra consulter la review des recherches récentes sur ce sujet donnée par Grover(1997).
Contrairement aux cas compétition et prédateur-proie, dans le modèle de Mu-
tualisme la présence mutuelle de deux espèces favorise la croissance de chacune
d’elles. Une telle symbiose a été notamment décrite en 1874 par le naturaliste
Thomas Belt qui a observé un mutulalisme formidable entre les fourmis et les ar-
bres acacias. Une formulation du modèle est donnée par le système suivant:
N1 − α12 N2
N10 = r1 N1 (1 − )
K1
N2 − α21 N1
N20 = r2 N2 (1 − )
K2
Pour l’étude théorique des systèmes écologiques en général, le lecteur pourra con-
sulter les auteurs suivants: Hirsch(1982)[88, 89, 90, 91, 92]: Systèmes d’équations
differentielles qui sont coopératif et compétitif, Smilatova & Sujan (1991): Systèmes
préda-teur-proie ( continu, discret, aléatoire, chaos , formes de stabilité) Logofet(1993)
:Différentes formes de stability passant par l’étude matricielle et les graphes, Kot(2000):
Une revue globale des mathémtiques pour l’écologie ( interaction, dispersion, struc-
ture d’âge) et références incluses.
133
5.5 Exercices
Exercice 5.5.1. Supposons que f est une fonction quadratique montrer pour une
reche-rche en ligne exacte on obtient
³ ´> ³ ´>
g(k) g(k) g(k) ρ(k)
αk = ³ ´> ; αk = − ³ ´> .
g(k) Gg(k) ρ(k) Gρ(k)
Exercice 5.5.2. Démontrer que le système suivant admet une solution unique
2 sin( x)
ex + + ye y = 0
x 2
3 y cos( y)
e − =0
y 2
2. Trouver une matrice D telle que xk −→ α (on pourra utliser la norme k.k∞ ).
xk+1 = xk − J −1 ( xk ) F ( xk )
134
(a) En posant sk = xk+1 − xk , écrire le système en sk qui permet d’eviter
d’inverser J.
(b) Rappeler l’expression de l’inverse P−1 d’une matrice de la forme
P = I − α uv> .
(c) Pour éviter le calcul de J −1 à chaque itération, on considère une matrice
Bk qui approche J et on utilise l’itération: xk+1 = xk − Bk−1 F ( xk ).
B−1 uv> Bk−1
Montrer que si Bk+1 = Bk + uv> alors Bk−+11 = Bk−1 − k .
1 + v> Bk−1 u
(d) On considère l’algorithme de Broyden qui consiste à prendre:
sk
u = yk − Bk sk et v = >
sk sk
f ( x) = 1/2x> Gx + b> x
xk+1 = xk + αk pk ; αk ∈ R et pk ∈ Rn
tel que
pi> Gp j = 0 pour tout i 6= j; i, j = 1, · · · , n
gi>+1 pi = 0 pour i = 1, · · · , n
135
Si les deux populations sont présentes on obtient le système suivant:
Xn+1 = g( Xn )
donné par:
xn+1 = ( 1 + r1 − b1 yn ) xn
yn+1 = ( 1 − r2 + b2 xn ) yn
2- Etudier la convergence.
3- Donner les conditions qui favorisent la première espèce par rapport à la deux-
ième.
136
Exercice 5.5.9. ( Modèle Maynard Smith)
Pn+1 = bNn Pn
Nn cPn
1- En posant xn = , yn = et a = bK on obtient le système
K bK
xn+1 = (1 + r) xn − rx2n − axn yn
yn+1 = axn yn
Problème 5.5.1. 1. Soit A une matrice carrée d’ordre n symétrique, définie pos-
itive. Deux vecteurs a et v sont dits A-conjugués si h Av, ai = 0.
(a) Montrer que si les vecteurs propres v0 , · · · , vn−1 sont A-conjugués, ils
forment une base de Rn .
(b) On définit les deux suites de matrices suivantes:
k−1 vi vi>
Ck = ∑ h Avi , vi i
i =0
Dk = I − Ck A.
Montrer que
Ck Av j = v j
Dk v j = 0
Dk> Av j = 0
Que valent Dn et Cn .
(c) Supposons que v0 , · · · , vk−1 soient connues.
Si Dk = 0 que peut-on conclure.
Sinon, soit v ∈ Rn tel que Dk v 6= 0, montrer que vk = Dk v est A-conjugué
par rapport à v0 , · · · , vk−1 .
(d) Ecire un algorithme qui à partir de v0 , donné constuit la suite v0 , · · · , vn−1 .
On peut considérer un vecteur de la forme Dk ei où ei est le ieme vecteur
de la base canonique.
En déduire un algorithme de calcul de A−1 .
137
Chapitre 6
6.1 Introduction
Comme on l’a vu dans les chapitres précédents, l’analyse spectrale ( calcul des
valeurs et vecteurs propres , rayon spectral, · · · ) joue un rôle important quant aux
décisions concernant la convergence , la stabilité etc · · · Par ailleurs, de nombreux
problèmes (statistique, modélisation,· · · ) se ramènent à l’étude spectrale des ma-
trices ou d’opérateurs en général .
Il est donc nécessaire de savoir calculer numériquement les valeurs et les vecteurs
propres dans des situations concrètes impliquant des grands systèmes où les métho-
des théoriques sont pratiquement inutilisables.
La stratégie générale pour déterminer les valeurs propres et les vecteurs propres
consiste à construire une suite de transformations de similarité jusqu’à l’obtention
d’une matrice diagonale, ou plus simplement jusqu’à l’obtention d’une matrice
triangulaire à partir de laquelle il est possible de déterminer assez facilement les
vecteurs propres .
pour réaliser ces transformations, deux grandes classes de méthodes
sont disponibles : les méthodes directes et les méthodes itératives.
Les premières consistent à effectuer une suite de transformations similaires et
ne s’appliquant que sur des matrices de taille modeste, car le temps de calcul croı̂t
comme le cube de la dimension linéaire de la matrice.
Pour les matrices de grande taille et généralement creuse, les méthodes itératives
sont plus adaptées. Dans la mesure ou la plupart du temps, le spectre complet
d’une très grande matrice n’est pas recherché mais seulement une partie, les métho-
138
des itératives peuvent réaliser aisément cette tâche.
Nous allons voir dans la suite de ce chapitre quelques algotithmes de base, tout
en ayant à l’esprit qu’il existe une vaste littérature sur ce sujet, et pour un problème
particulier il est nécessaire de commencer par analyser précisément le type de ma-
trice dont on souhaite obtenir le spectre, afin de choisir la méthode la plus adaptée
pour résoudre ce problème.
Soit Pn (λ ) = det( A − λ I )
et P(λ ) = (−1)n Pn (λ )
= λ n − a1 λ n−1 − · · · − an−1 λ − an
n
= λn − ∑ ak λ n−k
k=1
n
On a : An x(0) = ∑ ak A n−k x(0)
k=1
Alors, en posant x(1) = Ax(0) et x(k) = Ax(k−1) = Ak x(0) ,
on obtient : x(n) = a1 x(n−1) + · · · + an x(0) , ou encore
(n) (n−1) (n−2) (1) (0)
x1 x1 x1 · · · x1 x1 a1
. ..
. .. ..
. . ··· ··· ··· . .
(n) (n−1) (n−2) (1)
(0)
x = x · · · ai (6.2.1)
i i x i x i x i
. . .
..
.. .. · · · · · · · · · ..
.
(n) (n−1) (n−2) (1) (0)
xn xn xn ··· xn xn an
139
Soit en condensant
Ba = An x(0) .
P = S−1 AS.
1 ani
où mn−1n−1 = et mn−1i = − pour i 6= n − 1.
ann−1 ann−1
140
On verifie que la matrice inverse Mn−−11 a la même forme que Mn−1 et qu’elle
est donnée par
1 0 ··· ··· 0
0 1 ··· ··· 0
.. .. .. ..
..
Mn−−11 = . . . . . .
..
an−11 an−12 · · · . an−1n
0 ··· ··· 0 1
La méthode de Danilevski permet aussi de calculer les vecteurs propres d’une ma-
trice A si ses valeurs propres sont connues.
En effet si λ est une valeur propre de A, elle est aussi valeur propre de la matrice
de Frobenius semblable à A Soit donc v = (v1 , v2 , · · · , vn )> un vecteur propre de P
associé à la valeur propre λ, alors en écrivant Pv = λ v et en identifiant composante
à composante, on obtient le système
p1 v1 + p2 v2 + · · · + pn vn = λ v1
v1 = λ v2
v2 = λ v3
.. .. ..
. . .
vn−1 = λ vn .
141
6.3 Méthodes itératives
Soient λ1 , λ2 , · · · , λn les valeurs propres de A. On suppose que |λ1 | > |λ2 | >
n
· · · > |λn |. Tout vecteur x(0) peut s’écrire sur la base des vecteurs propres x(0) = ∑ αi vi .
i =1
³ ´ Ax(k−1)
On considère alors les deux suites x(k) et mk : x(k) = où mk est l’élément
° mk°
° °
de Ax(k−1) le plus grand en valeur absolue,mk = ° Ax(k−1) ° de telle sorte que 1
∞
est la plus grande composante de x(k) .
³ ´
Théorème 6.3.1. Dans ces conditions, les suites mk et x(k) convergent vers λ1 et v1 ,
respectivement, où v1 est le vecteur propre associé à λ1 .
Preuve : " µ ¶k #
n n n
λi
soit x(0) = ∑ αi vi , Ak x(0) = ∑ αi λik vi = λ1k α1 v1 + ∑ αi
λ1
vi
i =1 i =1 i =1
d’où
Ax(k−1)
x(k) =
mk
1 1
= A2 xk−2
mk mk−1
Ak x(0)
=
k
∏ mi
i =1
k
1 k (0)
=
Ck
A x avec Ck = ∏ mi
i =1
" µ ¶k #
λk n
λi
Par suite x(k) = 1 α 1 v 1 + ∑ αi vi donc lorsque k tend vers l’infini x(k)
Ck i =1
λ1
λ1k α1 λ1k
tend vers α1 v1 . Par identification, on obtient lim = 1 et donc lim x(k) = v1 ,
Ck k−→∞ Ck k−→∞
(k−1)
Ax λ1 v1
comme lim Ax(k) = Av1 = λ1 v1 on a lim x(k) = lim = d’où
k−→∞ k−→∞ k−→∞ mk lim mk
k−→∞
lim mk = λ1 .
k−→∞
4 1 1 1 6
(0) ( 0 )
Exemple 6.3.1. A = 1 0 −1 x = 1 , Ax = 0 m1 = 6 donc
1 −1 4 1 4
142
1 14/3 1
x(1) = 0 , et Ax(1) = 1/3 , m2 = 14/3, x(2) = 1/14 , · · ·
2/3 11/3 11/14
λ1
• Si les valeurs propres ne sont pas bien séparées, i.e ' 1, la convergence
λ2
sera très lente.
Méthode LR
Soit A une matrice carrée d’ordre n, la méthode LR est basée sur le procédé suiv-
ant: A1 = A = LR où L est triangulaire inférieure avec lii = 1 et U triangulaire
supérieure, Ai = Li Ri , Ai+1 = Ri Li . On suppose que toutes les décompositions
Li Ri existent. Alors on a:
3. Ai = Ti Ui avec Ui = Ri Ri+1 · · · R1 .
Preuve:
143
Théorème 6.3.3. On suppose que les valeurs propres de A vérifient: |λ1 | > |λ2 | >
· · · > |λn | > 0 et que la matrice de passage P ( P−1 AP = D ) ainsi que P−1 admettent les
décomposition P = L P R P et P−1 = Q = LQ RQ , (( L P )ii = 1 et ( LQ )ii = 1). Alors les
suites de matrices ( Ai ), ( Li ) et ( Ri ) convergent et on a :
λ1 ∗ ∗
..
lim Ai = lim Ri =
0 . ∗ et lim Li = I.
i −→∞ i −→∞ i −→∞
0 0 λn
Preuve:
¡ ¢
Soient A = PDQ et Ai = PDi Q = L P R P Di LQ RQ qui s’écrit A = L P R P Di LQ D −i Di RQ .
µ ¶i
i − i i
λj
Comme D LQ D = (l jk ) = l jk avec |λ j | < |λk | pour j ≥ k, il s’ensuit que
λk
lim l ijk = 0 pour j 6= k; On peut donc écrire Di LQ D −i = I + Ei avec lim Ei = 0.
i −→∞ i −→∞
Maintenant on remplace Di LQ D −i par I + Ei dans l’expression de Ai , on obtient:
Ai = L P R P ( I + Ei ) D i R Q
³ ´
= L P I + R P Ei R − P
i
R P Di RQ
= L P ( I + Fi ) R P Di RQ en posant Fi = R P Ei R−
P
1
λ1 ∗ ∗ · · ·
..
lim Ai = lim Li Ri = R P DR−
P
1
=
. ∗···
i i
λn
144
Méthode QR
Ai = Qi Ri
ii) Q1 · · · Qi Ri · · · R1 = Ai1 ,
iii) la décomposition QR de Ai1 est unique si R est triangulaire supérieure avec diagonale
strictement positive.
Théorème 6.3.4.
lim Ai = T
i→∞
où T est triangulaire supérieure contenant en diagonale les valeurs propres λ j rangées par
ordre décroissant en module.
P−1 A1 P = D = diag(λ j )
145
Méthode de Jacobi
Théorème 6.3.5. Soit θ un nombre réel et p et q tels que 1 ≤ p < q ≤ n deux entiers
aux quels on associe une matrice élémentaire Ω
1. Si A est une matrice symétrique alors la matrice B = Ω> AΩ est également symétrique
et vérifie ∑ bi2j = ∑ ai2j soit encore k Bk2E = k Ak2E .
ij ij
π π
2. Si a pq 6= 0, il existe une seule valeur de θ dans I =] − , 0[∪]0, [ telle que
4 4
b pq = 0, c’est la solution dans I de l’équation:
aqq − a pp
cotg(2θ ) = .
2a pq
n n
θ ainsi choisi, on obtient ∑ bii2 = ∑ aii2 + 2a2pq .
i =1 i =1
Preuve:
° °2 ¡ ¢ ¡ ¢
1. k Bk2E = °Ω> AΩ° E = k Ak2E = tr A> A = tr A2 .
146
Méthode de Jacobi classique
Théorème 6.3.6. La suite ( Ak+1 ) des matrices obtenues par la méthode de Jacobi classique
converge vers diag(λσ (i) ) où σ est une permutation convenable.
Preuve: ³ ´
(k) 2
Posons εk = ∑ ai j = k Ek k2E d’après le théorème 6.3.5 on a: εk+1 = εk −
i 6= j
³ ´2 ³ ´ ³ ´ ³ ´
(k) (k+1) 2 (k+1) 2 (k) 2
2 a pq ∑ = ε +
ai j
k+1 ∑ ii a = ε k + ∑ ii . a
i6= j i =1 i =1
¯ ¯ ¯ ¯
¯ (k) ¯ ¯ (k) ¯
Comme on a supposé que ¯ a pq ¯ = max ¯ ai j ¯ on a :
i, j
³ ´2 ³ ´2 ³ ´ µ ¶
(k) (k) 2 (k) 2 2
ε k = ∑ ai j ≤ ∑ a pq = (n − n) a pq d’où εk+1 ≤ 1 − ε
i6= j i6= j
n(n − 1) k
i.e εk+1 ≤ αεk avec 0 < α < 1 ce qui entraı̂ne que lim εk = 0 et donc lim k Ek k = 0.
k−→∞ k−→∞
Lemme 6.3.3. Soit E un espace vectoriel de dimension finie et ( xk ) une suite bornée de
E, admettant un nombre fini de valeurs d’adhérence et telle que: lim k xk+1 − xk k = 0.
k−→∞
Alors la suite ( xk ) converge vers l’une des valeurs d’adhérence.
Lemme 6.3.4. La suite ( Dk ) vérifie les hypothèses du lemme 6.3.3, par conséquent ( Dk )
converge vers une des valeurs d’adhérence, qui est de la forme diag(λσ (i) ).
Preuve: ³ ´
(k+1) (k)
Puisque Dk+1 − Dk = diag aii − aii alors on a
(
(k+1) (k) 0 si i 6= p i 6= q
aii − aii = (k)
± tan(θk ) a pq si i=p ou i = q.
π
avec |θk | < .
¯ ¯ 4 ³ ´
¯ (k) ¯ (k+1) (k)
¯ a pq ¯ ≤ k Ek k d’où lim aii − aii = 0 parsuite k Dk k E ≤ k Ak kE = k Ak k E est
k
bornée. Dk n’a qu’un nombre fini des valeurs d’adhérence, qui sont nécessairement
de la forme diag(λσ (i) ).
λ1 · · · λn supposées rangées dans un certain ordre .
147
en effet:
si ( Dk0 ) est une suite extraite de ( Dk ) et telle que lim
0
Dk0 = D, alors
k
lim
0
Ak0 = lim
0
( Dk0 + Bk0 ) = D et
k k
Donc A et D ont mêmes valeurs propres. Comme D est diagonale, elle contient les
λi en diagonale.
Finalement les lemmes 6.3.2, 6.3.3 et 6.3.4 entraı̂nent que
λ1
..
lim Ak = lim Dk = .
k k
λn
148
6.4 Applications
Notions de persistance et d’extinction
Définition 6.4.1. Une population sera dite persistante si lim n(t) > 0 pour toute
t−→+∞
valeur initiale n(0) et elle tendra vers l’extinction si lim n(t) = 0 pour un certain
t−→+∞
n ( 0 ).
Théorème 6.4.1. Soit L une matrice à coefficients constants non négatifs et réguliere avec
valeur propre dominante λ
149
Exemple 6.4.1. (Modèle ”Island”)
1 − (n − 1)d1 d2 ··· d2
d2 1 − (n − 1)d1 · · · d2
D=
.. .. .. ..
. . . .
d2 ··· ··· 1 − (n − 1)d1
λ ( Nc , p) = 1
on obtient µ ¶
1−d
Nc = E 1 + .
d1 − d2
avec
µ ¶1/k
1
d= .
λ1
Nc est le nombre maximal des régions qui assure la persistance puisque λ (n, p)
est décroissante.
150
La valeur critique Nc vérifie l’équation
µ µ ¶¶k
πN
λ1 1 − 2d1 − 2d2 cos = 1.
N+1
1 − b ≤ 2(d1 + d2 )
avec µ ¶1/k
1
b= .
λ1
µ ¶
b
arccos 1 − 2d1 − 2d2
Nc = E
µ ¶ + 1.
b
π − arccos 1 − 2d1 −
2d2
Nc est la valeur minimale qui assure la persistance.
151
6.5 Exercices
Exercice 6.5.1. 1. Rappelez la méthode de Householder permettant d’obtenir
une matrice réelle orthogonale .P telle que:
Px = ke1
Qi Ri = Ai − pi I;
Ai+1 = Ri Qi + pi I, i = 1, 2, · · · .
où Qi est une matrice réelle orthogonale, Ri une matrice triangulaire supérieure
et pi un paramètre de translation ( shift).
Exercice 6.5.2. Soit A une matrice symétrique ayant les valeurs propres (λi )in=1 as-
sociées aux vecteurs propres (ui )in=1 .
On pose ρ(µ , x) = Ax − µ x ∀ x ∈ Rn (ρ = 0 si x est vecteur propre associé à µ )
1. Montrer que le spectre de A contient au moins une valeur propre λi telle que
kρ(µ , x)k
Montrer qu’il existe une constante α telle que k x − α ui k ≤ .
d − kρ(µ , x)k
Exercice 6.5.3. Soit A une matrice dont les vecteurs propres v1 , · · · , vn sont associés
aux valeurs propres λ1 , · · · , λn vérifiant
152
1. On suppose que la valeur propre dominante λ1 et le vecteur propre associé v1
aient été déterminés et soit
v1 v>
1
A1 = A − λ1
v> v
1 1
λ1
w1 = v 1 , wi = vi − v1 x> vi i 6= 1.
λi
153
Chapitre 7
154
Alors pour tout α ∈ Rm , il existe une solution unique y( x) du problème (7.1.3), où y est
continue differentiable pour tout couple ( x, y) ∈ D.
0 0 0
En posant: Y1 = y, Y2 = Y1 (= y ), · · · , Yq = Yq−1 (= y(q−1) ) on obtient
0
³ ´>
Y = F ( x, Y) avec Y = Y> 1 , Y > , · · · , Y>
2 q ∈ Rqm
³ ´>
F = Y>
2 , Y >
3 , · · · , Y >
q , ϕ >
∈ Rqm
Exemple 7.1.1.
y(iv) = f ( x, y), f : R × R −→ R
0 0 0 00 0 000 0
Y1 = y, Y2 = Y1 = y , Y3 = Y2 = y , Y4 = Y3 = y , Y4 = y(iv)
f ( x, y) = A( x) y + ψ( x) (7.2.1)
155
Si les valeurs propres de A sont distinctes, la solution du système homogène est de
la forme
k
y( x) = ∑ C j exp(λ j x)V j + ϕ(x) (7.2.4)
j=1
où λ j est valeur propre de A associée au vecteur propre V j et les C j sont des con-
stantes arbitraires et ϕ( x) est une solution particulière de l’équation (7.2.2).
Exemple 7.2.1.
0
y = Ay + ψ( x) , y(0) = α
à !
1 2
avec A = , ψ( x) = (2x − 1, x + 1)> , α = (0, 0)> , les valeurs propres de
2 1
A sont λ1 = 3 et λ2 = −1.
Les vecteurs propres de A sont V1 = (1, 1)> et V2 = (1, −1)> . Ã !
ax + b
En cherchant une solution particulière de la forme: ϕ( x) = , on
cx + d
obtient
à ! à ! à !
1 1 −5/3
y( x) = C1 exp(3x) + C2 exp(− x) + .
1 −1 − x + 4/3
dX (t)
= F ( X (t)) (∗∗)
dt
156
On parle aussi de point stationnaire, état stationnaire et solution stationnaire
qui désignent la même notion.
2. Un point d’équilibre est dit asymptotiquement stable s’il est stable et en plus
∃ R0 tel que pour tout X (t0 ) ∈ B( X̄, R0 ) on a lim X (t) = X̄;
t−→+∞
3. Un point d’équilibre est dit marginalement stable s’ il est stable et non asymp-
totiquement stable;
Théorème 7.3.1. Si F est linéaire, une condition nécessaire et suffisante pour que le
système (∗∗) admette 0 comme point d’équilibre asymptotiquement stable est que
Re(λ ) < 0, pour toute valeur propre λ de F.
Si au moins une des valeurs propres de F vérifie Re(λ ) > 0 alors le point d’équilibre 0 est
instable .
Définition 7.3.3. Si le système non linéaire (∗∗) admet un point d’équilibre X̄, on
appelle matrice de linéarisation du système la matrice
∂ f1 ∂ f1
∂X1
···
∂f ∂Xn
2 ∂ f2
···
A= ∂X
.1
∂Xn ,
.. .. ..
. .
∂f ∂f
n n
···
∂X1 ∂Xn X̄
dY (t)
le système = AY (t) est dit système linearisé obtenu à partir du système
dt
(∗∗).
Théorème 7.3.2. Si le système non linéaire (∗∗) admet l’origine comme point fixe unique
alors dans un voisinage de l’origine, les comportements du système non linéaire et du
système linearisé sont équivalents à condition que le système n’admette pas de point centre
(valeurs propres imaginaires pures).
157
Remarque 7.3.1. Le théorème peut s’appliquer aux points d’équilibre qui sont dis-
tincts de l’origine, il suffit d’introduire des coordonnées locales.
i) V ( X̄ ) = 0 et V ( X ) > 0 si X 6= X̄,
ii) V̇ ( X ) ≤ 0 ∀ X ∈ U \ { X̄ }.
Alors X̄ est stable. Si de plus
iii) V̇ ( X ) < 0 ∀ X ∈ U \ { X̄ }
Pour les preuves des théorèmes 7.3.1, 7.3.2 et 7.3.3, voir ([10, 88, 103]).
La solution générale du système (7.4.1) s’obtient de façon similaire à celle qui donne
la solution du système d’équations differentielles linéaires (7.2.1). Elle est donnée
par la somme d’une solution (Yn ) du système homogène (7.4.2) et d’une solution
particulière ϕn du système (7.4.1) (yn = Yn + ϕn ).
158
Si r1 est une racine de π (r) de multiplicité m, la solution générale du système
m k
(7.4.1) est donnée sous la forme: yn = ∑ n j−1θ j rn1 + ∑ θ j rnj + ϕn .
j=1 j=m+1
1
Exemple 7.4.1. yn+2 − yn+1 + yn = 1, K = 2 est une solution particulière
2
1
π (r ) = r 2 − r + , ∆ = 1 − 2 = i 2 .
2
Les racines sont complexes conjugués r1 = (1 + i )/2 et r2 = (1 − i )/2. La solution
générale est donnée par yn = θ1 (1 + i )n /2n + θ2 (1 − i )n /2n + 2.
xi = a + ih, i = 0, 1, · · · , N
b−a
où h = est le pas de discrétisation (ici le pas est supposé constant mais on
N
peut envisager des pas hi variables).
La solution exacte au point xi est notée y( xi ), la solution approchée est notée yi
( y( xi ) ' yi ), une méthode numérique est un système d’équations aux differences
impliquant un certain nombre d’approximations successives yn , yn+1 , · · · , yn+k où
k désigne le nombre de pas de la méthode.
Si k = 1 on parle de méthode à un pas.
Si k > 1 la méthode est dite à pas multiples ou multi-pas.
Exemple 7.5.1.
yn+1 − yn = h f n
h
yn+2 + yn+1 − 2yn = ( f n+2 + 8 f n+1 + 3 f n )
4
h
yn+2 − yn+1 = ( 3 f n+1 − 2 f n )
3
h
yn+1 − yn = ( k1 + k2 )
4
159
avec µ ¶
1 1
k1 = f n = f ( xn , yn ) et k2 = f xn + h, yn + hk1 + hk2 .
2 2
Ces exemples peuvent être donnés sous la forme générale
k
∑ α j yn+ j = φ f ( yn+k , · · · yn , xn ; h ) . (7.5.2)
j=0
7.5.1 Convergence
Considérons une méthode numérique de type (7.5.2) avec des valeurs initiales ap-
propriées:
k
∑ α j yn+ j = φ f ( yn+k , · · · yn , xn ; h ) (7.5.3)
j=0
yκ = γκ (h) , κ = 0, 1, · · · , k − 1
Définition 7.5.1. La méthode (7.5.3) est dite convergente si pour tout problème de
condition initiale vérifiant les hypothèses du théorème 7.1.1 ( existence et unicité)
on a
max k y( xn ) − yn k = 0 quand h −→ 0
0≤n≤ N
k
Rn+k ( h ) = ∑ α j y(xn+ j ) − hφ f ( y(xn+k ), y(xn+k−1 ), · · · , y(xn ), xn , h) (7.5.4)
j=0
1
Remarque 7.5.1. Parfois on utilise aussi τn (h) = Rn+k (h) comme définition d’erreur
h
de troncature locale, mais dans le cadre de ce chapitre c’est la première définition
qui est adoptée.
¡ ¢
Définition 7.5.3. La méthode (7.5.3) est dite d’ordre p si Rn+k = O h p+1 .
7.5.2 Consistance
La méthode numérique est dite consistante si son ordre est au moins 1 ou encore,
pour tout p.c.i vérifiant les hypothèses du théorème 7.1.1 d’existence et d’unicité
on a
1
lim Rn+k (h) = 0, x = a + nh.
h−→0 h
160
k
i) ∑ α j = 0.
j=0
k
ii) (φ f ( y( xn ), y( xn ) · · · y( xn ), xn , 0))/ ∑ jα j = f (xn , y(xn )).
j=0
7.5.3 Stabilité
z0 = f ( x, z), z( a) = α , x ∈ [ a, b]
z0 = f ( x, z) + δ ( x), z( a) = α + δ , x ∈ [ a, b] .
161
Définition 7.5.7. On dit que la méthode numérique satisfait les conditions aux
racines si tous les zéros du 1er polynôme catactérisitique sont de module inférieur
ou égal à 1 et ceux ayant un module égal à 1 sont des zéros simples.
Théorème 7.5.1. Une condition nécessaire et suffisante pour que la méthode soit zéro-
stable et que la méthode vérifie la condition aux racines.
Théorème 7.5.2. Une condition suffisante et nécessaire pour que la méthode (7.5.2) soit
convergente est qu’elle soit zéro- stable et consistante.
Pour les preuves des théorèmes 7.5.1 et 7.5.2, voir ([48, 103]).
La plus simple et la plus connue des méthodes d’approximation des solutions des
e.d.o est la méthode d’Euler donnée par:
yn+1 − yn = h f n , y0 = α (7.5.5)
En considérant
1 2
y( xn + h) − y( xn ) − h f ( xn , y( xn )) = h y”(θn ) avec xn ≤ θn ≤ xn+1 (7.5.6)
2
On voit que la méthode d’Euler est une méthode explicite d’ordre 1 (donc consis-
tante).
0 ≤ (1 + x)m ≤ exp(mx)
t
z0 ≥ − et zn+1 ≤ (1 + s) zn + t ∀n = 1, · · · , k
s
Alors on a µ ¶
t t
zn+1 ≤ exp((n + 1)s) + z0 − .
s s
Théorème 7.5.3. Soit f : R × R −→ R une fonction continue et vérifiant la condition
de Lipschitz pour y sur D = {( x, y); a ≤ x ≤ b, −∞ < y < ∞} avec a et b finis. On
suppose qu’il existe une constante M telle que | y”( x)| ≤ M pour tout x ∈ [ a, b] . Soit
0
y( x) la solution unique du p.c.i y ( x) = f ( x, y) , y( a) = α et ( yn )nn= N
=0 la suite des
approximations générée par la méthode d’Euler. Alors on a:
hM
| y( xn ) − yn | ≤ (exp( L( xn − a)) − 1) pour chaque n = 0, 1, · · · , N
2L
162
Preuve.
Pour n = 0, l0 inégalité est vérifiée puisque y( x0 ) = y0 = α .
Les équations (7.5.5) et (7.5.6) donnent:
1
| y( xn+1 ) − yn+1 | ≤ | y( xn ) − yn | + h| f ( xn , y( xn ) − f ( xn , yn )| + h2 | y”(θn )| (7.5.7)
2
Les hypothèses du théorème conduisent alors à la mojoration
h2 M
| y( xn+1 ) − yn+1 | ≤ | y( xn ) − yn |(1 + hL) +
2
h2 M
En appliquant le lemme avec zn = | y( xn ) − yn |, s = hL et t = il vient:
2
µ ¶
h2 M h2 M
| y( xn+1 ) − yn+1 | ≤ exp(((n + 1)hL) | y( x0 ) − y0 | + −
2hL 2hL
h2 M
Preuve. Identique à celle du théorème avec y( x0 ) − w0 = δ0 et t = + |δn |.
2
Remarque 7.5.4. La simplicité de la méthode d’Euler en fait un exemple pédagogi-
que d’introduction aux autres méthodes plus complexes. Cependant, un des critères
principaux de l’analyse numérique consiste à chercher des méthodes ayant l’ordre
de précision le plus élevé possible et comme la méthode d’Euler est d’ordre 1, son
utilisation se trouve limitée en pratique et on est amené à considérer des méthodes
plus précises. Trois directions principales permettent d’obtenir des méthodes d’ordres
élevés.
La première direction consiste à travailler avec des méthodes à un pas mais en
cherchant à atteindre des ordres élevés en utilisant un développement de Taylor
163
et en négligeant le terme d’erreur mais ces méthodes ont un handicap à cause des
dérivées susccessives de f .
Une deuxième possibilité est donnée par des choix appropriés de
φ f ( yn+k , · · · yn , xn ; h) dans l’équation (7.5.1), les méthodes de Runge-Kutta sont la
meilleure illustration de cette direction.
Enfin, une troisième direction est offerte par les Méthodes Linéaires à Pas Mul-
tiples (MLPM).
En se basant sur le critère de précision, on voit qu’on est obligé de chercher des
méthodes dont la performance est supérieure à celle d’Euler. On fait donc appel à
d’autres méthodes plus précises.
h2 00 hn hn+1 (n+1)
y( xn+1 ) = y( xn ) + hy0 ( xn ) + y ( xn ) + · · · + y(n) ( xn ) + y (ζn ),
2! n! (n + 1)!
164
Remarque 7.5.5. Au vu du critère de précision et bien que les méthodes de Taylor
paraissent faciles dans leur écriture, elles sont rarement utilisées dans la pratique à
cause des difficultés engendrées par le calcul des dérivées successives de f comme
fonction de deux variables. C’est pour cette raison qu’on cherche des méthodes
permettant d’atteindre un ordre élevé tout en évitant le calcul des dérivées succes-
sives de f . Au critère de précision s’ajoute le critère de coût.
l
y n+1 − y n = h ∑ bi ki
i =1
à !
l l
avec ki = f xn + ci h, yn + h ∑ ai j k j , i, = 1, 2, · · · , l et on suppose que ci = ∑ ai j , i, = 1, 2, · · · , l
j=1 j=1
Si ai j = 0 pour j > i alors:
à !
i
ki = f xn + ci h, yn + h ∑ ai j k j i, = 1, 2, · · · , l
j=1
on obtient ainsi:
k1 = f ( xn , yn )
k2 = f ( xn + c2 h, yn + c2 hk1 )
Remarques 7.5.6. 1. Les méthodes de R.K sont des méthodes à un pas, elles
peuvent-être écrite sous la forme générale yn+1 − yn = hφ f ( xn , yn , h), où
l
φ f ( x, y, h) = ∑ bi ki .
i =1
2. Les méthodes de R.K satisfont la condition aux racines, elles sont donc zéro-
stables. Par conséquent, pour étudier la convergence il suffit d’étudier la con-
sistance.
yn+1 = yn + h ( b1 k1 + b2 k2 + b3 k3 + b4 k4 )
165
en déterminant pour chaque ordre les valeurs possibles des paramètres.
k1 = f ( xn , yn )
k2 = f ( xn + c2 h, yn + c2 hk1 )
k3 = f ( xn + c3 h, yn + (c3 − a32 )hk1 + a32 hk2 )
.. .. ..
. . .
Les formes explicites des méthodes d’ordre 1,2 et 3 sont obtenues aprés differenti-
ation et identification.
Méthode d’ordre 1
En considérant yn+1 = yn + hb1 k1 avec k1 = f ( xn , yn ) on voit que le l’odre maximal
est 1 obtenu avec le paramètre b1 égal à 1 et qui donne la méthode d’Euler:
yn+1 − yn = h f ( xn , yn )
Méthodes d’ordre 2
En partant de yn+1 = yn + hb1 k1 + b2 k2 avec k1 = f ( xn , yn ) et
k2 = f ( xn + c2 h, yn + c2 k1 ), p = 2 est l’ordre maximal possible qu’on peut atteindre
1
si les paramètres b1 , b2 et c2 vérifient les équations b1 + b2 = 1 et b2 c2 = . Comme
2
on a deux équations pour trois inconnues, il en résulte une infinité de méthodes
explicites d’ordre 2 . On en donne quelques unes:
1
Méthode d’Euler modifiée: (b1 = 0, b2 = 1 et c2 = )
2
µ ¶
1 1
yn+1 − yn = h f xn + h, yn + k1
2 2
1
Méthode d’Euler améliorée: ( b1 = b2 = et c2 = 1)
2
h
yn+1 − yn = ( f ( xn , yn ) + f ( xn + h, yn + k1 ))
2
1 3 2
Méthode de Heun d’ordre 2 ( b1 = , b2 = et c2 = )
4 4 3
µ ¶
h 2 2
yn+1 = yn + f ( xn , yn ) + 3 f ( xn + h, yn + k1 )
4 3 3
k1 = h f ( xn , yn )
Méthodes d’ordre 3
En considérant yn+1 = yn + b1 k1 + b2 k2 + b3 k3 avec k1 = h f ( xn , yn ) ,
k2 = h f ( xn + c2 h, yn + c2 hk1 ) et k3 = h f ( xn + c3 h, yn + (c3 − a32 )k1 + a32 k2 ).
166
On obtient une famille de méthode d’ordre 3 si les paramètres vérifient les
équations
b1 + b2 + b3 = 1
1
b2 c2 + b3 c3 =
2
1
b2 c22 + b3 c23 =
3
1
b3 c2 a32 =
6
Deux représentants de cette famille de méthode d’ordre 3 sont :
1 3 1
Méthode de Heun d’ordre 3 (b1 = , b2 = 0, b3 = et c2 = et c3 =
4 4 3
2 2
, a32 = )
3 3
h
yn+1 = yn + (k1 + 3k3 )
4
k1 = h f ( xn , yn )
µ ¶
1 1
k2 = h f xn + h, yn + k1
3 3
µ ¶
2 2
k3 = h f xn + h, yn + k2
3 3
1 2 1 1
Méthode de Kutta d’ordre 3 (b1 = , b2 = , b3 = et c2 = et c3 = 1, a32 = 2)
6 3 6 2
h
yn+1 = yn + (k1 + 4k2 + k3 )
6
k1 = h f ( xn , yn )
µ ¶
1 1
k2 = h f xn + h, yn + k1
2 2
k3 = h f ( xn + h, yn − k1 + 2k2 )
Méthode d’ordre 4
Méthode de Runge-Kutta d’ordre 4
h
yn+1 = yn + (k1 + 2k2 + 2k3 + k4 )
6
k1 = h f ( xn , yn )
µ ¶
1 1
k2 = h f xn + h, yn + k1
2 2
µ ¶
1 1
k3 = h f xn + h, yn + k2
2 2
k4 = h f ( xn + h, yn + k3 )
167
Méthode de Runge-Kutta d’ordre 4
h
yn+1 = yn + (k1 + 3k2 + 3k3 + k4 )
8
k1 = h f ( xn , yn )
µ ¶
1 1
k2 = h f xn + h, yn + k1
3 3
µ ¶
2 1
k3 = h f xn + h, yn − k1 + k2
3 3
k4 = h f ( xn + h, yn + k1 − k2 + k3 )
25 1408 2197 1
yn+1 = yn + k1 + k3 + k4 − k5
216 2565 4104 5
avec
k1 = h f ( xn , yn )
µ ¶
1 1
k2 = h f xn + h, yn + k1
4 4
µ ¶
3 3 9
k3 = h f xn + h, yn + k1 + k2
8 32 32
µ ¶
12 1932 7200 7296
k4 = h f xn + h, yn + k1 − k2 + k3
13 2197 2197 2197
µ ¶
439 3680 845
k5 = h f xn + h, yn + k1 − 8k2 + k3 − k4
216 513 4104
µ ¶
1 8 3544 1859 11
k6 = h f xn + h, yn − k1 + 2k2 − k3 + k4 − k5
2 27 2565 4104 40
168
pour estimer l’erreur de troncature locale de la métohde de R-K d’ordre 5 :
13 2375 5 12 3
yn+1 = yn + k1 + k3 + k4 + k5 + k6
160 5984 16 85 44
avec
k1 = h f ( xn , yn )
µ ¶
1 1
k2 = h f xn + h, yn + k1
6 6
µ ¶
4 4 16
k3 = h f xn + h, yn + k1 + k2
15 75 75
µ ¶
2 5 8 5
k4 = h f xn + h, yn + k1 − k2 + k3
3 6 3 2
µ ¶
5 165 55 425 85
k5 = h f xn + h, yn − k1 + k2 − k3 + k4
6 64 6 64 96
µ ¶
12 4015 11 88
k6 = h f xn + h, yn + k1 − 8k2 + k3 − k4 + k5
5 612 36 255
µ ¶
1 8263 124 643 81 2484
k7 = h f xn + h, yn − k1 + k2 − k3 − k4 + k5
15 15000 75 680 250 10625
µ ¶
3501 300 297275 319 24068 3850
k8 = h f xn + h, yn + k1 − k2 + k3 − k4 + k5 + k7
1720 43 52632 2322 84065 26703
k
0
L( z( x); h) := ∑ (α j z(x + jh) − hβ j z (x + jh))
j=0
Si on suppose que z est une fonction qu’on peut dériver autant de fois qu’on
veut, on obtient:
169
Définition 7.5.9. La MLPM est dite d’ordre p si
C0 = C1 = · · · = C p = 0 et C p+1 6= 0.
7.6 Applications
La modélisation et la simulation par ordinateur ont fait de l’outil mathématique
un moyen important pour la compréhension, l’analyse et le contrôle de certaines
maladies. La formulation d’un modèle permet de tester des hypothèses, d’estimer
des paramètres, de construire des théories, de discuter des conjonctures, de for-
muler des scénarios prédictifs, de visualiser certaines sensibilités et de remédier
à l’impossibilité de certaines expériences pour coût élevé ou danger opérationnel.
La modélisation mathématique permet de mieux saisir les concepts d’épidémie, de
seuil et des nombres de contacts, de replacements ou de reproduction. Le modèle
mathématique fournit également le support nécessaire à la réalisation d’appareils
et d’équipements d’analyse et de contrôle dans le domaine médical. A ce propos,
nous renvoyons à une review récente de Hethcote [85] et aux 200 références citées
par l’auteur.
Dans ce chapitre, nous considérons trois types d’applications, la première traite
les maladies d’infection transmise de façon directe, la deuxième est consacrée à
la transmission par vecteur et la troisième traite l’effet de l’effort physique sur les
dynamiques d’insuline et de glucose.
Pour les modèles à transmission directe, nous nous limiterons aux cas discrets.
• S(t) : nombre des susceptibles à l’instant t (on désigne par susceptibles, les
individus qui peuvent avoir la maladie);
• E(t) : nombre des exposés à l’instant t (on désigne par exposés, les individus
qui sont infectés mais ne sont pas infectieux);
• I (t) : nombre des infectieux à l’instant t(les personnes qui sont déjà touchées
par la maladie et peuvent transmettre la maladie);
170
• R(t) : nombre des résistants à l’instant t (les personnes qui sont guéries avec
immunisation )
Chacune des classes précédentes est dite compartiment, et suivant les caractéristiques
de chaque maladie on a les types de modèles suivants : SI, SIS, SIR, SEIR, SEIRS,· · · )
Remarque 7.6.1. • En général les classes M et E sont souvent omises.
• une population est de taille fixe si le taux de naissance est égal au taux de
mortalité.
Exemple 7.6.1. Modèle SIS
C’est un modèle à deux compartiments utilisé essentiellement pour décrire les mal-
adies sexuellement transmissibles. Les guéris ne développent pas une immunité et
redeviennent susceptibles à la maladie: deux classes de sous populations suscep-
tibles et infectieux. La dynamique des sous populations est régie par le système
d’équations suivant:
Sn+1 = Sn (1 − λ .∆t.In ) + γ .∆t.In + µ .∆t(1 − Sn )
In+1 = In (1 − γ .∆t − µ .∆t + λ .∆t.Sn )
I0 + S0 = 1 avec S0 ≥ 0 et I0 > 0
où Sn et In sont les proportions d’individus susceptibles et infectieux, respective-
ment.
Comme Sn + In = 1 alors In+1 = In (1 − (γ + µ ).∆t + λ .∆t − λ .∆t.In ).
• si (λ − γ )∆t < 2 alors on a convergence vers le point fixe;
171
Le modèle SIS Le modèle SIS
0.7 0.9
0.8
0.6
0.7
0.5
0.6
Nombre des infectieux
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0 0
0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50
Temps Temps
0.9
1.2
0.8
1
0.7
Nombre des infectieux
0.6
0.8
0.5
0.6
0.4
0.3
0.4
0.2
0.2
0.1
0 0
0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50
Temps Temps
172
Exemple 7.6.2. Modèle SIR
Le modèle SIR est un modèle simple de transmission des maladies par contact entre
les individus ,ces derniers deviennent immunisés après une seule infection, une
troixième classe de sous populations est considérée: la classe R des isolés ou enlevés
”removed” ce sont les personnes qui, une fois infectés soit deviennent immunisés
soit ils meurent. Le modèle SIR est plus adapté pour les maladies d’enfants comme
la rougeole. Les équations qui régissent ce modèle sont :
Sn+1 = Sn (1 − λ .∆t.In ) + µ .∆t(1 − Sn )
I = In (1 − γ .∆t − µ .∆t + λ .∆t.Sn )
n+1
Rn+1 = Rn (1 − µ .∆t) + γ .∆t.In
S + I + R = 1 avec S > 0, I > 0 et R ≥ 0
0 0 0 0 0 0
40
Nombre des infectieux
15
30
10
20
5
10
0 0
0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50
Temps Temps
173
(SIR à deux populations)
80 160
70 140
Nombre des infectieux
50 100
40 80
30 60
20 40
10 20
0 0
5 10 15 20 25 30 35 40 45 50 5 10 15 20 25 30 35 40 45 50
Temps
Figure 7.4: SIR à deux populations Temps
174
Exemple 7.6.3. Modèle SEIR
Modèle SEIR avec λ indépendante du temps
Le modèle SEIR est une variante du modèle SIR où la période de latence est prise
en compte c.à.d une classe de la population est intermédiaire entre la classe des
susceptibles et les infectieux: la classe des exposés qu’ un susceptible traverse pour
passer à l ’état d’infection. Un autre paramètre est introduit qu’on note β et on
1
parle de période de latence . Les équations régissant ce modèle sont:
β
Sn+1 = Sn (1 − λ .In .∆t) + µ .∆t(1 − Sn )
En+1 = En + Sn λ .In .∆t − (β + µ ).∆t.En )
In+1 = In (1 − γ .∆t − µ .∆t) + β.∆t.En
Rn+1 = Rn (1 − µ .∆t) + γ .∆t.In
S0 + I0 + R0 + E0 = 1
175
Le modèle SEIR Le modèle SEIR
35 35
30
30
25
25
Nombre des infectieux
15
15
10
10
5
5
0
0 −5
0 50 100 150 200 250 300 0 50 100 150 200 250 300
Temps Temps
Le modèle SEIR
50
40
30
20
Nombre des infectieux
10
−10
−20
−30
−40
−50
0 50 100 150 200 250 300
Temps
176
ii) Pour β = 45, γ = 73, µ = 0.02, ∆t = 0.027, I0 = 0.001, E0 = 0.01 et S0 = 0.25
il y’a oscillation voir figure 7.7.
40
500
30
400
Nombre des infectieux
300 10
0
200
−10
100
−20
0 −30
0 50 100 150 200 250 300 0 50 100 150 200 250 300
Temps Temps
h
h ?
Rh h-
Figure 7.8: Schéma de transmission
177
Le modèle suppose un mixage homogène d’humain et de vecteur de telle sorte
que chaque piqure a la même probabilité d’être appliquée à un humain particulier.
En notant bs le taux moyen de piqures d’un vecteur susceptible, phv la probabilité
moyenne de transmission d’un infectieux humain à vecteur susceptible, le taux
d’exposition pour les vecteurs est donné par: (phv Ih bs )/ Nh .
en notant pvh la probabilité moyenne de transmission d’un infectieux vecteur à
un humain et Iv le nombre de vecteurs infectieux, le taux d’exposition pour les
humains est donné par: ( pvh Iv bi )/ Nh par conséquent:
• Le taux de contact adéquat d’humain à vecteurs est donné par: Chv = phv bs
La durée de vie humaine est prise égale à 25 000 jours (68.5 ans), et celle des
vecteurs est de: 4 jours. Les autres paramètres sont donnés dans le tableau 7.1
d’après [141].
178
Population vecteur
dSv = µv Nv − (µv + Chv Ih / Nh ) Sv
dt
dIv = (Chv Ih / Nh ) Sv − µv Iv
dt
Avec les conditions Sh + Ih + Rh = Nh , Sv + Iv = Nv , Rh = Nh − Sh − Ih et
Sv = Nv − Iv
Le système s’écrit
dSh
= µh Nh − (µh + p + Cvh Iv / Nh ) Sh
dt
dIh
= (Cvh Iv / Nh ) Sh − (µh + γh ) Ih
dt
dIv = C I / N ( N − I ) − µ I
hv h h v v v v
dt
Pour l’étude de stabilité et d’autres détails concernant les paramètres du modèle,
le lecture pourra consulter le papier cité en réferences. Dans le cadre restreint de ce
chapitre, nous discutons très brièvement l’output du modèle.
La diffuculté principale de la dengue provient du fait qu’elle est causée par qua-
tre virus différents et que l’ummunité provisoire acquise contre un virus ne protège
pas contre les autres virus, au contraire, un individu attaqué par un deuxième virus
court le danger d’évoluer vers le stade de Dengue Hemorragique.
La recherche de vaccin se complique par le fait que ce vaccin doit couvrir le
spectre des quatre virus et ceci constitue un problème.
Le modèle montre que la réduction du nombre de mosquitos (par insecticides
ou autre) n’est pas suffisante pour éradiquer l’épidémie de la Dengue. Elle peut
tout au plus retarder l’apprition de l’épidémie. Comme la découverte d’un vaccin
global n’est pas attendue dans le court terme, les auteurs suggèrent la combinaison
du contrôle des facteurs environnementaux et de vaccins partiels contre chaque
virus.
179
number of the infectious humans first epidemic
6000
N =30000
4000 v
N =50000
v
Nv=8000
2000
−2000
0 100 200 300 400 500 600
time
Nv=30000
4000
Nv=50000
population
Nv=8000
2000
−2000
700 750 800 850 900 950 1000
time
Figure 7.9:
180
Durant les dernières décennies, une littérature abondante a été consacrée à ce
sujet, le lecteur interessé pourra consulter les deux reviews récentes par Bellazi et
al (2001) et Parker et al. (2001).
Incorporant l’effet de l’effort physique sur la dynamique du glucose et de l’insuline,
Derouich et Boutayeb (2002) ont proposé le modèle suivant
dG (t) = −(1 + q2 ) X (t) G (t) + ( p1 + q1 )( Gb − G (t))
dt
dI (t)
= ( p3 + q3 )( I (t) − Ib ) + p2 X (t)
dt
où ( I (t) − Ib ) représente la différence entre l’insuline dans le plasma et l’insuline
de base (Gb − G (t)) représente la différence entre la concentration de glucose dans
le plasma et le glucose de base X (t) est l’insuline interstitiale, q1 , q2 ,q3 sont des
paramètres liés à l’effort physique.
Les auteurs ont discuté les résultats de ce modèle dans trois cas différents:
Dans les trois cas, le modèle met en évidence l’intérêt de l’effort physique à améliorer
la sensibilité de l’insuline mais le dernier cas reste le plus illustratif comme l’indique
la figure 7.10 où on voit qu’un diabétique non insulino-dépendant pourrait étre
amené à vivre continuellement avec 2g/ml de glucose dans le sang sans se rendre
compte des conséquences à long terme de cette overdose. Avec un effort physique,
le diabétique peut ramener la courbe de glucose au voisinage de la concentration
normale de 1g/ml.
181
glucose
300
250
200
glucose
150
100
50
0 20 40 60 80 100 120 140 160 180 200
time
0.04
interstitial insulin
plasma insulin
100
0.03
0.02
50
0.01
0 0
0 50 100 150 200 0 50 100 150 200
time time
182
aléatoires conduisant à des séries chronolgiques. D’autres vont défendre les modéles
déterministes sous-jaçant l’idée du chaos apparaissant comme du stochastique ou
un bruit blanc dans les séries chronologiques. Les travaux de May (1974, 75) et
plus particulièrement son excitant papier publié dans la revue Nature en 1976 sous
le titre de ”Simple mathematical models with very complicated dynamics” con-
stitueront des réferences pour les écologistes. Par la suite, Anderson et May vont
publier une multitude de papiers et de livres dévoués à l’épidémiologie des mal-
adies infectieuses.
Hoppensteadt (1975) est considéré par certains auteurs comme le premier à
avoir proposé une analyse mathématique des modèles avec structure d’âge.
Durant les dernières décennies, la littérature consacrée à la modélisation
épidémiologique a connu un essort considérable comme en atteste la dizaine de re-
views (Becker, 1978; Castillo-Chavez, 1989; Dietz, 1967-75-85; Wickwire, 1977; Het-
hcote, 1981-94-2000; Isham, 1988) et la trentaine de livres (Anderson & May, 1982-
91; Bailey, 1975-82; Bartlett, 1960; Castillo-Chavez, 1989; Diekman, 2000; Frauen-
thal, 1980; Grenfell and Dobson, 1995; Hethcote, 1984-92; Isham & Medeley, 1996;
Kranz, 1990; Lauwerier, 1981; Ludwig & Cooke, 1975; Mollison, 1991; Scott &
Smith, 1994). D’autres ouvrages et reviews traitent les modèles épidémiologiques
comme partie des modèles écologiques (Dietz K. (1975), Pielou (1977), Okubo (1980),
Kot (2000)).
183
7.8 Exercices
Exercice 7.8.1. On considère le modèle discret suivant:
Rn+1 = (1 − b) In + Rn
1. La population des susceptibles tend vers une limite S quand n tend vers
l’infini.
2. Si F = S/ S0 alors
Exercice 7.8.3. Chercher les solutions des systèmes differentiels et aux differences
0
suivants: y = Ay + ψ( x) , y(0) = α avec
à ! à ! à !
1 1 1 2
A= , ψ( x) = , α=
−1 1 x 1/2
à !
n
yn+4 − 6yn+3 + 14yn+2 − 16yn+1 + 8yn = φn avec yn ∈ R2 et φn = .
1
Exercice 7.8.4. Considérons le système aux differences :
yn+2 − 2µ yn+1 + µ yn = c, n = 0, 1, · · ·
avec yn , c ∈ Rm et 0 ≤ µ ≤ 1.
Montrer que yn converge vers c/(1 − µ ) quand n −→ ∞.
0 ≤ (1 + x)m ≤ exp(mx)
184
ii) Si s et t sont des réels positifs et ( zn )nn= k
=0 une suite vérifiant:
t
z0 ≥ − et zn+1 ≤ (1 + s) zn + t ∀ n = 1, · · · , k
s
Alors on a:
µ ¶
t t
zn+1 ≤ exp((n + 1)(1 + s)) + z0 − .
s s
0
B- Soit y( x) la solution unique du p.c.i y ( x) = f ( x, y) ,a ≤ x ≤ b , y( a) = α
et w0 , w1 , · · · , w N , les approximations vérifiant: w0 = α + δ0 et wn+1 = wn +
h f ( xn , wn ) + δn+1 ; n = 0, ..., N − 1
| f ( x, y) − f ( x, y∗ )| ≤ L| y − y∗ |∀ y, y∗ ∈ D
ii) il existe une constante M telle que:| y”( x)| ≤ M pour tout x ∈ [ a, b] .
iii) les perturbations verifient: |δn | < δ ∀n = 0, · · · , N
Montrer que:
µ ¶
1 hM δ
| y( xn ) − wn | ≤ + (exp( L( xn − a)) − 1) + |δ0 | exp L( xn − a)
L 2 h
Exercice 7.8.6. Soit A une matrice diagonalisable admettant les vecteurs propres V j
associés aux valeurs propres λ j avec Re(λ j ) < 0 pour tout j.
On considère les deux schémas d’Euler donnés par:
E1 : un+1 = un + hAun , u0 = η
E2 : un+1 = un + hAun+1 , u0 = η
185
Montrer que les solutions de ces deux schémas sont données par:
k
S1 : un = ∑ c j (1 + hλ j )n V j
j=1
k
S2 : un = ∑ c j ( 1 − h λ j )−n V j
j=1
Quelles conditions doit-on imposer à h pour que les solutions S1 et S2 tendent vers
zéro ou restent bornées quand n −→ ∞
k k
∑ α j yn+ j = h ∑ β j f n+ j
j=0 j=0
αk = 1 et |α0 | + |β0 | 6= 0
k
0
L( z( x); h) := ∑ (α j z(x + jh) − hβ j z (x + jh)) = C0 z(x) + C1 hz0 (x) + · · · + Cq hq z(q) (x) + · · ·
j=0
1. Montrer que:
k
C0 = ∑ α j = ρ(1)
j=0
k
C1 = ∑ ( jα j − β j ) = ρ0 (1) − σ (1)
j=0
k µ ¶
1 q 1 q−1
Cq = ∑ q
j αj −
(q − 1)
j β j , q = 2, 3, · · ·
j=0
186
Exercice 7.8.8. Chercher l’ordre des méthodes numériques suivantes:
h
1) yn+2 = yn + ( f n + 4 f n+1 + f n+2 )
6
4h
2) yn+4 = yn + ( 2 f n+1 − f n+2 + 2 f n+3 )
3
µ ¶
h 2 2
3) yn+1 = yn + f ( xn , yn ) + 3 f ( xn + h, yn + k1 ) , k1 = h f ( xn , yn )
4 3 3
h
4) yn+1 = yn + (k1 + 4k2 + k3 )
6
k1 = h f ( xn , yn )
µ ¶
1 1
k2 = h f xn + h, yn + k1
2 2
k3 = h f ( xn + h, yn − k1 + 2k2 )
h
yn+2 − ( 1 + α ) yn+1 + α yn = ((3 − α ) f n+1 − (1 + α ) f n )
2
où f n = f ( xn , yn ) et −1 ≤ α ≤ 1.
en supposant que y0 = 1 + ω h3
y1 = exp(h) + θ h3
187
3. On suppose que r1 est de la forme
h2
r1 = 1 + h + + O( h3 )
2
(a) Donner une expression analogue pour r2 .
(b) Etudier la convergence de yn dans le cas α = −1.
où f n = f ( xn , yn ) et f ( xn , y( xn )) = y0 ( xn ).
y0 ( x) = − y( x), y(0) = 1 ( P)
yn = C1 (r1 )n + C2 (r2 )n ,
r2 − exp(−h) exp(−h) − r1
avec C1 = et C2 = .
r2 − r1 r2 − r1
µ ¶
2 4 2 1/2 1 1 1 1 4
6. On donne 1 + h + h = 1 + h + h2 − h3 + h + O( h5 ).
3 9 3 6 18 216
1 1 1 4
En déduire que r1 = 1 − h + h2 − h3 + h + O(h5 ), r2 = −5 − 3h +
2 6 72
O( h2 ).
188
1 4
7. Montrer alors que C1 = 1 + O(h2 ) et C2 = − h + O( h5 ).
216
8. En considérant pour x fixé, nh = x, montrer que C1 (r1 )n −→ exp(− x) quand
n −→ ∞.
189
Annexe
Cette annexe est un guide d’initiation à MATLAB. MATLAB est un programme in-
teractif de calcul scientifique utilisable pour la résolution numérique de nombreux
problèmes mathématiques ou appliqués. En outre, MATLAB dispose de poten-
tialités graphiques importantes.
L’objectif est de permettre au débutant de se familiariser avec MATLAB qui est
avant tout un programme de calcul matriciel.
Toutefois, nous nous sommes limités aux méthodes les plus courantes décrites
dans les chapitres précédents, il est recommandé d’exécuter les exemples dans
l’ordre où ils apparaissent.
function [n]=mat_square(A)
[n,m] = size(A);
if n ~= m
disp(’ Error: the matrix should be squared’);
stop
end
return
Broyden
function [x,it]=broyden(x,Q,nmax,toll,f)
[n,m]=size(f); it=0; err=1;
fk=zeros(n,1); fk1=fk;
for i=1:n, fk(i)=eval(f(i,:)); end
while it < nmax & err > toll
s=-Q\fk; x=s+x; err=norm(s,inf);
if err > toll
for i=1:n, fk1(i)=eval(f(i,:)); end
Q=Q+1/(s’*s)*fk1*s’;
end
190
it=it+1; fk=fk1;
end
Choleski
Jacobi
function [xvect,err,nit]=jacobi(x0,n_max,tol,A,b)
P=diag(diag(A));
N=P-A; B=eye(size(A))-inv(P)*A; bb=P\b; r=b-A*x0; erreur=norm(r);
nit=1; xvect(:,1)=x0; err(1)=erreur;
return
Gauss-Seidel
function [xvect,err,nit]=gauss_seidel(x0,n_max,tol,A,b)
P=tril(A); F=-triu(A,1); B=inv(P)*F; bb=P\b; r=b-A*x0;
erreur=norm(r); nit=1; xvect(:,1)=x0; err(1)=erreur;
191
erreur = norm(r);
err(nit+1) =erreur;
nit = nit + 1;
end
return
Householder
function [H,Q]=houshess(A) n=max(size(A));
Q=eye(n); H=A; for
k=1:(n-2),
[v,beta]=vhouse(H(k+1:n,k));
I=eye(k);
N=zeros(k,n-k);
m=length(v);
R=eye(m)-beta*v*v’;
H(k+1:n,k:n)=R*H(k+1:n,k:n);
H(1:n,k+1:n)=H(1:n,k+1:n)*R;
P=[I, N; N’, R];
Q=Q*P;
end return
Horner
function [pnz,b] = horner(a,n,z)
b(1)=a(1);
for j=2:n+1,
b(j)=a(j)+b(j-1)*z;
end; pnz=b(n+1);
LU
function [L,U,P,Q] = LUpivtot(A,n)
P=eye(n); Q=P; Minv=P; for
k=1:n-1
[Pk,Qk]=pivot(A,k,n);
A=Pk*A*Qk;
[Mk,Mkinv]=MGauss(A,k,n);
A=Mk*A;
P=Pk*P;
192
Q=Q*Qk;
Minv=Minv*Pk*Mkinv;
end U=triu(A); L=P*Minv;
Gradient
function [x, error, niter, flag] = gradient(A, x, b, M, maxit,tol)
flag = 0;
niter = 0;
bnrm2 = norm( b );
if ( bnrm2 == 0.0 ),
bnrm2 = 1.0;
end
r = b - A*x;
error = norm( r ) / bnrm2;
if ( error < tol )
return,
end
for niter = 1:maxit
z = M \ r;
rho = (r’*z);
q = A*z;
alpha = rho / (z’*q );
x = x + alpha * z;
r = r - alpha*q;
error = norm( r ) / bnrm2;
if ( error <= tol ),
break,
end
end
if ( error > tol )
flag = 1;
end
return
Substitution
function [b]=backward_col(U,b)
[n]=mat_square(U);
193
l = min(diag(abs(U)));
if l == 0
disp(’Error: the matrix is singular’);
b = [];
break
end
for j = n:-1:2,
b(j)=b(j)/U(j,j);
b(1:j-1)=b(1:j-1)-b(j)*U(1:j-1,j);
end;
b(1) = b(1)/U(1,1);
return
Relaxation
Sturm
function [p]=sturm(dd,bb,x)
n=length(dd);
194
p(1)=1;
p(2)=d(1)-x;
for i=2:n
p(i+1)=(dd(i)-x)*p(i)-bb(i-1)^2*p(i-1);
end
return
Point fixe
function [xvect,xdif,fx,nit]=fixpoint(x0,nmax,toll,fun,phi)
err=toll+1;
nit=0;
xvect=x0;
x=x0;
fx=eval(fun);
xdif=[];
while (nit < nmax & err > toll),
nit=nit+1;
x=xvect(nit);
xn=eval(phi);
err=abs(xn-x);
xdif=[xdif; err];
x=xn;
xvect=[xvect;x];
fx=[fx;eval(fun)];
end;
Gradient conjugué
function [x, error, niter, flag] = conjgrad(A, x, b, P, maxit, tol)
flag = 0;
niter = 0;
bnrm2 = norm( b );
if ( bnrm2 == 0.0 ),
bnrm2 = 1.0;
end
r = b - A*x;
error = norm( r ) / bnrm2;
if ( error < tol )
return,
195
end
for niter = 1:maxit
z = P \ r;
rho = (r’*z);
if niter > 1
beta = rho / rho1;
p = z + beta*p;
else
p = z;
end
q = A*p;
alpha = rho / (p’*q );
x = x + alpha * p;
r = r - alpha*q;
error = norm( r ) / bnrm2;
if ( error <= tol ),
break,
end
rho1 = rho;
end
if ( error > tol )
flag = 1;
end
return
décomposition LU
196
end
end
return
Méthode de Newton
function [xvect,xdif,fx,nit]=newton(x0,nmax,toll,fun,dfun)
err=toll+1;
nit=0;
xvect=x0;
x=x0;
fx=eval(fun);
xdif=[];
while (nit < nmax & err > toll),
nit=nit+1;
x=xvect(nit);
dfx=eval(dfun);
if (dfx == 0), err=toll*1.e-10;
disp(’ Stop ’);
else,
xn=x-fx(nit)/dfx;
err=abs(xn-x);
xdif=[xdif; err];
x=xn;
xvect=[xvect;x];
fx=[fx;eval(fun)];
end;
end;
197
end
[L,U,P]=lu(Jxn);
else
for i=1:n, Fxn(i)=eval(F(i,:)); end
end
nit=nit+1; step=step+1; Fxn=-P*Fxn; y=forward_col(L,Fxn);
deltax=backward_col(U,y); x = x + deltax; err=norm(deltax);
if nit > nmax
disp(’ pas de convergence’);
break
end
end
198
Bibliographie
[1] Acevado M. (1994), Fundamentals of species interaction: Two and three communi-
ties, University of North Texas, Denton, Texas.
[2] Acevado M. (1994), Two species interaction, stage structure and environmental
factors, University of North Texas, Denton, Texas.
[3] Akerman E. et al. (1965), Model studies of blood glucose regulation, Bulletin of
mathematics and Biophysics 27: 21-37.
[4] Allen L. J. (1992), Some discrete-time SI, SIR, SIS, epidemic models, Lub-
bock,Texas.
[5] Allen L. J. (1994), Discrete-time population with age and structure I, II, Lub-
bock,Texas.
[8] Anderson R. M., May R. M. (1991), Infectious disease of humans: Dynamics and
control, Oxford University Press, Oxford.
199
[14] Bailey N. T. J. (1975), The mathematical theory of infectious diseases, 2nd ed.,
Hafner, NewYork.
[15] Bajpai et al. (1978), Numerical methods for engineers and scientists, John Wiley.
[21] Becker N. (1989), Analysis of infectious diseases data, Chapman and Hall, New
York.
[23] Berman A., Plemmons R. J. (1994), Non negative matrices in the mathematical
sciences, Classics in Applied Mathematics 9, SIAM, Philadelphia.
[24] Bernardelli H. (1941), Population waves, Journal of the Burma Research Soci-
ety 31, 1-18
[25] Bernouilli D. (1760), Essai d’une nouvelle analyse de mortalité causée par la petite
vérole, Académie Royale des Sciences, Paris, 1-45
[29] Boutayeb A., Twizell E. H. (1992), Numerical methods for the solution of special
sixth-order bvp, Intern. J. Computer Math. 45: 207 - 223.
200
[30] Boutayeb A., Kerfati A. (1994), Mathematical models in diabetology, Modelling,
Measurement & Control, C, AMSE Vol.44, N◦ 2: 53-63.
[31] Boutayeb A., Derouich M. (2002), Age strucured models for diabetes in East
Morocco, Mathematics and Computers in simulation 58: 215- 229,
[32] Boutayeb A., Chetouani A. (2003), Global extrapolation of numerical methods for
solving a parabolic problem with non local boundary conditions, Intern. J. Com-
puter Math 80, n◦ 6, 789-797.
[34] Brauer A. (1957), A new proof of theorems of Perron and Frobenius on nonnegative
matrices, Ducke Math. J. 24: 367-378.
[36] Bugl P. (1995), Differential equation matrices and models, Prentice Hall Inc. New
Jersey.
[38] Carnahan B., Luther H. and Wilkes J. (1969), Applied numerical methods, John
Wiley, New York.
[40] Caswell H. (2001), Matrix populations model, Sinauer associates, Inc., Sunder-
land, MA.
[44] Conte S. D., de Boor C. (1978), Elementary numerical analysis, Third edition,
McGraw-Hill International Book Company.
201
[45] Couffignal (1956), Résolution numérique des systèmes d’équations linéaires, Ey-
rolles.
[49] Dautray, Lions (1984), Analyse mathématique et calcul numérique pour les sci-
ences et les techniques, Masson.
[51] Della C. et al. (1970), On a mathematical model for the analysis of the glucose
tolerance curve, Diabetes 19: 445-449.
[53] Derouich M., Boutayeb A., Twizell E. H. (2003), A dengue fever model, Bio-
medical engeneering online 2:4.
[54] Derouich M., Boutayeb A. (2002), The effect of physical exercice on the dynamics
of glucose and insulin, Journal of Biomechanics 35 : 911-917.
[56] Dietz K. (1967), Epedemics and rumours: A survey, J.Roy. Statist. Soc. Ser.A,
130: 505-528.
[58] Dietz K., Schenzle D. (1985), Mathematical models for infectious disease statstics,
in A celebration of statstics Atkinson eds, Spring Verlag, New York.
[59] Djidjeli K., Twizell E. H., Boutayeb A. (1993), Numerical methods for special
nonlinear bvp of order 2m, J. Comp. Appl. Math. 47: 35 - 45
202
[60] Djidjeli K., Twizell E. H. (1991), Global extrapolations of numerical methods for
solving a third-order depersive partial differential equation, Intern. J. Computer
Math, Vol 41: 81-89.
[61] Durand (1972), Solution numérique des équations algébriques tome I, Masson.
[62] Durand (1972), Solution numérique des équations algébriques tome II, Masson.
[64] Evans D.J. (1984), Parallel SOR Iterative Methods, Parallel Computing 1: 3-18.
[65] Fisher R. A. (1930), The genitical theory of natural selection, Dover, New York.
[67] Frobenius G. (1908), Uber matrizen aus positiven elelmenten, S.B Preuss Akad.
Wiss Berlin, 471-76
[71] Gerschgorin S. (1931), Uber die Abgrenzung der eigenwerte einer matrix, Izv
Akad. Nauk SSSR Ser Mat. 7: 749-754.
[73] Golub G. H, Van Loan C. (1983), Matrix computations, Johns Hopkins Uni-
veristy Press, Baltimore.
[74] Golub G. H. (1989), Some history of the conjugate gradient and Lanczos algo-
rithms, Siam review, vol. 31.
[75] Golub G. H. (1983), Résolution numérique des grands systèmes linéaires, Eyrolles.
[76] Golub G.H. & Varga R.S. (1961), Chebychev Semi-Iterative Methods, Successives
Over- Relaxation Iterative methods, and second-Order Richardson Iterative Meth-
ods. Parts I & II, Numer. Math. 3: 147-56; 157-68.
203
[77] Gourdin A., Boumahrat M. (1988), Méthodes numériques appliquées, Office
Press Université Alger.
[78] Grenfell B. T., Dobson A. P. EDS. (1995), Ecology of infectious diseases in natural
populations, Cambridge University Press, Cambridge, UK.
[81] Hestenes, Stiefel (1952), Methods of conjugate gradients for solving linear sys-
tems, J. Res. Nat. Bur. Standards 49.
[84] Hethcote H. R. (1994), A thousand and one epidemic models, in frontiers in the-
oretical biology, Levin ed., Lectures notes in Biomath. 100 , Spring- Verlag,
Berlin.
[87] Himsworth H. P. & Ker R. B. (1939), Insulin sensitive and insulin insensitive
types of diabetes mellitus, Clinical Science4: 119-125.
[88] Hirsch M. W. & Smale S. (1974), Differential equations, Dynamical Systems, and
linear Algebra, Academic Press. Inc.
204
[91] Hirsch M. W. (1988), Systems of differential equations which are competitive or
cooperative: III. Competing species, Non linearity 1: 51-71.
[97] Horn R., Jonshon R. (1985), Matrix analysis, Cambridge University Press,
Cambridge, England.
[101] Isham V. (1988), Mathematical modeling of the transmission dynamics of HIV and
AIDS: A review, J.Roy. Statist.Soc. SEr. A, 151: 5-31.
[102] Isham V., Medeley G., EDS. (1996), Models for infectious human diseases, Cam-
bridge University Press, Cambridge, UK.
[103] Issaacson E., Keller H. B. (1966) Analysis of numerical methods, John Wiley.
[104] Jansen H. (1994), The Numerical modelling of measles dynamics, Brunnel, Uni-
versity, Uxbridge, England.
[105] Kendall D. G. (1949), Stochastic processus and population growth J, Roy. Statist.
Soc. B 11: 230-64.
205
[106] Kendall D. G. (1956), Deterministic and stochastic epidemics in closed populations,
Proc. 3rd Berkeley Symp. Math. Stat. Prob. Vol VI: 149-165.
[108] Kershau (1978), The incomplete Choleski conjugate gradient method for the itera-
tive solution of systems of linear equations, J. Comp. Phys. 26.
[111] Kranz J. ED. (1990), Epidemics of plant diseases: Mathematical analysis and mod-
elling, Spring-Verlag, Berlin.
[113] Lambert J. D.(1991), Numerical methods for ordinary differential systems, John
Wiley & Sons.
[114] Lascaux P., Théodor R. (1986), Analyse numérique matricielle appliquée à l’art de
l’ingénieur tome 1, Edition Masson.
[115] Lascaux P., Théodor R. (1987), Analyse numérique matricielle appliquée à l’art de
l’ingénieur tome 2, Edition Masson.
206
[121] Leslie P. H. (1948), Some further notes on matrices in population mathematics,
Biometrika, 35: 213-245
[126] Lotka A. J. (1922), The stabilty of the normal age distribution, Proc. Nat. Acad.
Sci. USA, 8: 339-345.
[127] Lotka A. J. (1931), The extinction of families, Journal of the Washington Acad-
emy of sciences: 377-380.
[128] Ludwig D., Cooke K. L. EDS. (1975), Epidemiology, SIMS Utah Conference
Proceedings, SIAM, Philadelphia.
[130] Lui Y. (1997), Numerical solution of the heat equation with nonlocal boundary
conditions, Technical Report Cambridge University.
[132] Marcus M., Minc H. (1964), A survey of matrix theory and matrix inequalities,
Allyn and Bacon , Boston.
[133] May R. M. (1976), Simple mathematical models with very complicated dynamics
nature, Lond. 261: 459-67
207
[136] Mendel G. J. (1865), Versuche uber Pflanzen-Hybriden. Verh, Naturforsch. Ver.
Brunn. 10.
[138] Minc H. (1988), Nonnegative matrices, John Wiley & sons, New York.
[139] Mollison D. (1996), Epidemic models: Their structure and relation data, Cam-
bridge University, Cambridge,UK.
[141] Newton E. A., Reiter P. (1992), A model of the transmission of dengue fever with an
evoultion of the impact of ultra-low volume (ULV) insecticide apllications on dengue
epidemics, Am J Trop Med Hyg, 47, 709-720.
[144] Ostrowski (1966), Solution of equations and systems of equations, 2nd edition,
Academic Press.
[145] Parker R.S et al.(2001), The intravenuous route toblood glucose control, IEEE
Engineering in Mediciner and biology 20: 65-73.
[146] Pearl R., Reed L. J. (1920), On the rate of growth of the population of the United
States since 1790, Proceeding of the National Academy of Sciences 6: 275-288.
[148] Pielou E. C. (1977), Mathematical ecology, John Wiely & Son New York.
[149] Plemmons R.J. (1986), A parallel block iterative Scheme Applied to Computations
in Structural Analysis, SIAM J Alg. and Disc Methods 7: 337-347.
[150] Qin Mengzhao (1983), Difference schemes for the dispersive equation, Comput-
ing, 31: 261-267.
[151] Rabinowtz (1970), Numerical methods for nonlinear algebraic equations, Gordon
and Breach.
208
[152] Reed M.(1987), Msc lectures notes, Brunel University Middx U.K.
[153] Reid J. (1970), On The method of conjugate gradient method of large sparse systems
of linear equations, Proc. of Oxford conf of Inst of Math. appl.
[154] Ross R. (1911), The prevention of Malaria, 2nd ed. , Murray, London.
[156] Scott M. E., Smith G. EDS. (1994), Parasitic and infectious diseases, Academic
Press, San Diego.
[159] Serge G. et al. ( 1973), Modelling blood glucose and insulin kinetics in normal
diabetic, Diabetes 22: 94-99 and obese subjects.
[160] Sibony M., Mrdon J. Cl. (1984), Analyse numérique I et II, Hermann éditeurs.
[161] Smilatova K., Sujan S. (1991), Dynamical methods in biological science: A math-
ematical treatement, Ellis Horwood.
[162] Stein P. & Rosenberg R. L.(1948), On the solution of linear simultaneous equation
by iteration, J. London Math. Soc. 23:111-118.
[163] Stewart, Jensen (1975), Solutions numériques des problèmes matriciels, Eyrolles,
(1975).
[167] Twizell E. H. (1988), Numerical methods with applications in the Biomedical sci-
ences, Ellis Horwood, Chichester and John Wiley & Sons, New york.
[168] Twizell E. H., Boutayeb A. (1991), Numerical methods for the solution of special
and general sixth-order bvp, with application to Bénard eigenvalue problems, Proc.
Roy. Soc. Lon. 431: 433 - 450.
209
[169] Varga R. (2000), Matrix iterative analysis 2nd edition, Springer, Berlin-Tokyo.
[170] Varah J. M. (1972), On the solution of block tridiagonal systems araising from
certain finite difference equations, Math. Comp. 26: 859-868.
[172] Vignes (1980), Algorithmique numérique tome II, équations et systèmes linéaires,
Technip.
[173] Von Foerster H. (1959), Some remarks on changing populations in the kinetics of
cellular proliferation, Math. Comp. 26: 382-407.
[174] Watson (1980), Approximation theory and numerical methods, John Wiley.
[176] White R.E. (1989), Multisplitting with different weighting schemes, SIAM J.
Matrix Anal. Appl. 10: 481-493.
[177] Wickwire K. (1977), Mathematical models for the controle of pests and infectious
diseases: A survey, Theoret. Population Biol., 11: 182-238.
[178] Wiggins S. (1990), Introduction to applied nonlinear systems and chaos, Spring-
Verlag, New York.
[180] Wilkinson J. H. (1988), The algebraic eigenvalue problem, Clarendon Press Ox-
ford.
[181] Yicang Z., Xiwen L. (1996), The dynamics and application of the nonlinear matrix
model of the disabled population, Science College, XiánJiaotang University.
210
Notes de solutions des exercices
Chapitre 1
exercice 1.6.1
i) ¯ ¯
¯ ¯ ¯ ¯
a x
∑i ¯∑ j i j j ¯ ∑i ∑ j ¯ ai j x j ¯ ¯ ¯
¯ ai j ¯ .
k Ak1 = sup ¯¯ ¯
¯ ≤ sup ¯ ¯
∑ j ¯x j ¯
≤ max
i
∑
x6=0 ¯∑ j x j ¯ x6=0 j
k Aei k1 = ∑ | ai j |,
j
iii) A> A est une matrice symétrique définie positive, considérons une base or-
thonormée (e1 , · · · , en ) de vecteurs propres de A> A.
h A> Ax, xi
k Ak2 = sup
x6=0 h x, xi
∑ λi | xi |2
i
= sup
x6=0 ∑ | xi |2
i
≤ λn = ρ( A> A).
211
2. evident.
3.
k Axk∗ k BAxk∞
k Ak∗ = sup = sup .
x6=0 k xk∗ x6=0 k Bx k∞
5. la jeme colone de A−1 est (2 j−2 , 2 j−3 , · · · , 1, · · · , 0)> 1 est situé à la jeme posi-
tion.
k Ak1 = n et k A−1 k1 = 2n−1 .
6. ¯ ¯
¯ ¯
¯ ¯
k Axk∞ = max ¯∑ ai j xi ¯ ≥ | aii | − ∑ | ai j | max | x j | ≥ δ k xk∞ .
j ¯ i ¯ i6= j j
On prend y = A−1 x.
7. k Ak∞ = 4 + ε, le 1 donne λ4 ≤ 4 + ε.
On prend dans le 2, B = I.
Chapitre 2
exercice 2.12.1
1.
(2) ai1
ai j = ai j − a1 j , i, j = 2, · · · , n.
a11
(2) a
aii = aii − i1 a1 j , i = 2, · · · , n.
a11
2.
¯ ¯ ¯ ¯ ¯ ¯
¯ (2) ¯ ¯ ai1 ¯ ¯ ai1 ¯
¯ a j j ¯ > | a j j | − ¯¯ a1 j ¯ > | a j j | − ¯¯ a1 j ¯¯
a11 ¯ i∑ 6= j a11
¯ ¯ ¯ ¯
¯ (2) ai1 ¯ ¯ ai1 ¯
¯
> ∑ ¯a j j + ¯
a1 j ¯ + | a1 j | − ¯ ¯ a1 j ¯¯
i6= j
a 11 a 11
¯ ¯ ¯ ¯
¯ (2) ¯ ¯ a1 j ¯
¯
> ∑ ¯ a j j ¯ − ∑ | ai1 | ¯ ¯¯ + | a1 j |
i6= j i6= j
a11
¯ ¯
¯ (2) ¯
> ∑ ¯ ai j ¯ .
i6= j
212
3. Il n’est pas nécessaire de pivoter
exercice 2.12.2
1
1. Si α 6= alors I − α ab> est inversible.
b> a
¡ ¢−1 A−1 ab> A−1
A + ab> = A−1 − où λ = b> A−1 a.
1+λ
¡ ¢−1
2. I + mk e>
k = I − mk e>
k .
Chapitre 3
exercice 3.8.1
i) A−1 = ( I − R)−1 B.
iii) x̂ − x = A−1 r.
2. Calcul facile.
° ° 1
3. (a) ( I − M)−1 b( x) est contractante car °( I − M)−1 °∞ ≤ 3 et kbk∞ ≤ .
4
24
(b) k = .
33
(c) On peut améliorer la vitesse de convergence en prenant εb( x) avec ε <<
1.
213
1. Utiliser le théorème 3.3.4.
exercice 3.8.6
I il suffit qu’il existe une norme matricielle telle que k I − α Ak < 1, ceci est
assuré si on prend k.k∞ .
II 1) 1 + λiα est valeur propre de I − α A par suite pour que la méthode con-
2
verge, il faut et il suffit que α < .
λ1
2
2) α∗ = .
λ1 + λn
exercice 3.8.7
2. θ = 0 et ω = 1 donnent Jacobi.
r = ω = 1 donnent Gauss-Seidel.
r = ω donne relaxation.
3. il suffit d’écrire.
det(α D − β L − U ) = det(α D − βµ L − µ −1 U ).
exercice 3.8.8
1. On suppose ici que les ai sont positifs. Vous pouvez le montrer directement
en regroupant dans la somme h A(h) x, xi des termes positifs, ou utiliser un
théorème qui assure qu’une matrice est définie positive si et seulement si ces
éléments diagonaux sont positifs.
2. Voir chapitre 3.
exercice 3.8.9
214
2. (a) T = (rI + H2 )−1 (rI − H1 )(rI + H1 )−1 (rI − H2 ), c = (rI + H2 )−1 (rI −
H1 )(rI + H1 )−1 (rI − H2 )b + (rI + H2 )−1 b.
¯ ¯ ¯ ¯
¯r − λj ¯ ¯r − µj ¯
¯
(b) ρ( T ) ≤ max ¯ ¯ max ¯ ¯ , où λi et µi désignent les valeurs pro-
j r + λ ¯ j ¯r + µ ¯
j j
pres de H1 et H2 .
(c) xk converge vers x∗ solution de ( I + T )−1 ( I − T ) x∗ = b.
¯ ¯ µ¯ ¯ ¯ ¯¶
¯ (r − x ) 2 ¯ ¯ (r − a ) 2 ¯ ¯ (r − b ) 2 ¯ √
(d) min max ¯ ¯ ¯ = min max ¯ ¯ , ¯ ¯ , et ropt = ab.
r ≥0 a ≤ x ≤b (r + x ) 2 ¯ r ≥0 a ≤ x ≤b ¯ (r + a ) 2 ¯ ¯ (r + b ) 2 ¯
exercice 3.8.10
a − ib
(b) λ = .
1 − a + ib
(c) On a x H L1 x + x H L> H −1/2 ( L + L> ) D −1/2 x = 2a, or A = D − L −
1 x = x D 1 1
L est définie positive, donc x D /2 Ax H D −1/2 > 0 ainisi 1 − 2a > 0 et
> H − 1
Chapitre 4
√
exercice 4.4.1 les points d’équilibre sont 0 et s−1 α si α > 0.
Si α > 0 alors 0 n ’est pas un point d’attraction.
Si α = 0 alors 0 est un point d’attraction si s > 1.
√
Le point x∗ = s−1 α est un point d’attraction si s < 1.
exercice 4.4.2
√
5−1
1. α = , β = 1 − α.
2
2. Posons L1 = b − a, x1 = a, x4 = b x2 = x0+ L1 β, x3 =x4 − L1 β.
x1 x2
x2 x3
Si f ( x3 ) < f ( x2 ) on remplace par
x x − αβ L
3 4 1
x4 x4
215
x1 x1
x2 x1 + αβ L1
Si f ( x3 ) > f ( x2 ) on remplace par
x3 x2
x4 x3
x1 + x2
min = .
x∈[ a,b] 2
exercice 4.4.6
1. Φ( x) = x3 − x2 − 1, |Φ0 ( x)| ≥ 1.
1 1
2. Φ( x) = 1 + + 2 , |Φ0 ( x)| < 1.
x x
exercice 4.4.3
1. Soit λ > 0.
√
Si 1 + b > 0 alors il existe un seul point d’équilibre P∗ = 2+b
λ, il est stable si
1
λ< .
1+b
Si b = −1, le système est stationnaire.
1
Si 1 + b < 0 alors il existe deux points d’équilibre P∗ = 0 et P∗∗ = √ , où
r−1
λ
1
r = −(1 + b), 0 est stable si λ < , P∗∗ est stable si −2 < b < −1.
−(1 + b)
1
2. Les points d’équilibre sont P∗ = 0 et P∗∗ = log λ .
α
0 est stable si 0 < λ < 1.
P∗∗ est stable si e−2 < λ < 1.
exercice 4.4.9
216
exercice 4.4.10
1. Affirmatif.
exercice 4.4.14
à r ! à r !
4rE 4rE
K r+ r2 − K r− r2 −
K K
a) Les points d’esuilbre sont P∗ = et P∗∗ = .
2r 2r
Le point P∗ est stable et P∗∗ est instable.
rK
b) Si E > , le système n’admet pas de points d’équilibre et P(t) est strictement
4
décroissante.
K (r − E )
c) Les points d’esuilbre sont P∗ = 0 et P∗∗ = , si r > E.
r
Le point P∗ est instable et P∗∗ est stable.
à r !
4rE
K (r − 1 ) + (r − 1 ) 2 −
K
c) Les points d’esuilbre sont P∗ = et
à r ! 2r
4rE
K (r − 1 ) − (r − 1 ) 2 −
K 4KE
P∗∗ = , si (r − 1)2 > .
2r r
4KE
Le point P∗ est stable si 0 < (r − 1)2 − < 4 et P∗∗ est instable.
r
exercice 4.4.7
1
a) < α < 51.
101
1
b) R = 51, r = .
101
217
c) Le nombre de racines négatives = 0 − 2k, donc une seule possibilité: 0.
exercice 4.4.4. Il suffit de tracer les graphes des fonctions f et g (les points d’equlibre
sont l’intersection des 2 graphes) et de dériver les conditions de stabilité en les im-
β x2
posant à la fonction Φ( x) = x + rx − .
α + x2
exercice 4.4.2
exercice 4.4.4
2.
exercice 4.4.13
1. (a) P( y) = g0 ( yn ) y + ( g( yn ) − g0 ( yn ) yn ).
y f −1 ( xn )
(b) P( y) = + ( xn − ).
f 0 ( f −1 ( xn )) f 0 ( f −1 ( xn ))
(c) Méthode de Newton.
2. Méthode de la sécante.
g00 ( yn ) 2 g00 ( yn ) 2
3. (a) P( y) = y + ( g( yn ) − g00 ( yn ) yn ) y + g( yn ) + yn − g0 ( yn ) yn .
2 2
f 00 o f −1 1
(b) g00 = − 00 −1 3 , g0 = 0 −1 .
(f of ) (f of )
f 000 (α ) f 00 (α )
(c) C = − + .
6 f 0 (α ) 2( f 0 (α ))2
Chapitre 5
1. a) J ( xk )sk = F ( xk ), xk+1 = xk + sk .
α uv>
b) P−1 = I + .
1 − α v> u
B−1α uv> Bk−1
c) Bk−+11 = Bk−1 − k .
1 + α v> Bk−1 u
218
d) Bk−+11 yk = sk .
2. a), b) ,c) définissent la méthode des directions conjuguées décrite das le chapitre
cinq.
Chapitre 6
Chapitre 7
exercice7.7.1
µ ¶
1 1
si α 6= 1 et β 6= − l’ordre p = 1 et C2 = (α − 1) β +
2 2
1 1
si α 6= 1 et β = − l’ordre p = 2 et C3 = (α − 1)
2 12
µ ¶
1 1
si α = 1 et β 6= − l’ordre p = 2 et C3 = − β +
2 2
1 1
si α = 1 et β = − l’ordre p = 3 et C4 = −
2 12
219
exercice7.7.2 On cherchera une solution particulière du système differentiel (resp.
du système aux differences) de la forme:
à ! à à !!
ax + b an + b
φ( x) = resp.ψn =
cx + d cn + d
1. α0 = −5, α1 = 4, β0 = 2, β1 = 4, ordre 3.
exrecice 7.7.4
c
Par passage à la limite yn −→ .
1−µ
exercice 7.7.5
exercice 7.7.8
1. Méthode d’ordre 1.
2.
3.
220
Chapitre 8
problème 8.0.1
1. Prendre M = I et N = I − A.
3. (a) a > 1.
0 −a −a
(b) J = a 0 − a , J converge si a < 1/2.
a a 0
0 0 0
(c) G = − a 0 0 , ρ( G ) = 0, si a 6= 0 la méthde de Gauss Seidel
−a −a 0
converge plus vite que la méthode de Jacobi.
Problème 8.0.2
3. Si Dk = 0 alors Ck = A−1 .
Problème 8.0.3
exercice 1
b b
0 − −
a a
b b
1. M J =
− a 0 −
a
b b
− − 0
a a
221
2b
2. Méthode convergente si k M J k∞ = < 1 donc a > 2b.
a
a + 2b
3. A est inversible car elle est à diagonale dominante, C ( A) ≤ .
a − 2b
exercice 2
1. Méthode d’ordre 1.
2. Méthode d’ordre 3.
exercice 2
3. (a) Si µ est valeur propre de M(0, 1) et λ est valeur propre de M(0, ω) alors
λ = 1 − ω + ωµ .
(b) Si la méthode de Jacobi converge alors |µi | < 1 pour tout i, A étant
symetrique on peut ranger ses valeurs propres par ordre croissant µ1 ≤
1 1
µ2 ≤ · · · ≤ µn , donc ≤ ∀i.
1 − µ1 1 − µi
La méthode (∗) converge si −1 < 1 − ω + ωµi < 1 et ceci est vérifié si
2
0<ω< .
1 − µ1
222
Index
A C
Acacias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 Carrée . 6, 22, 24, 31, 33, 47–49, 51–53, 56,
Age . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15–17 73–75, 77, 86, 142
Akerman . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 Castillo-Chavez . . . . . . . . . . . . . . . . . . . . . . 182
Algorithme . . . 25, 37, 38, 41, 42, 46, 47, 53, Caswell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
75, 86, 91, 109, 120, 123, 124, 134 Cauchy . . . . . . . . . . . . . . . . . . . 73, 83, 97, 153
Anderson . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 Cayley-Hamilton . . . . . . . . . . . . . . . . . . . . 138
Applications . . . . . . . . . . . . 15, 102, 125, 169 Chebyshev . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
Approximation . . . . . . . . . . . . . . . . . . . . 82, 85 Choleski . . . . . . . . . . . . . . . . . . 34, 35, 46, 190
Ciarlet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Axelsson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Classique . . . . . . . . . . . . . . . . . . . . . . 53, 60, 78
Collatz . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21, 72
B Comparaison . . . . . . . . . . . . . . . . . . . . . . . . . 61
Bailey . . . . . . . . . . . . . . . . . . . . . . . . . . 181, 182 Conditionnement . . . . . . . . . . . 36, 46, 75, 77
Bartlett . . . . . . . . . . . . . . . . . . . . . . . . . 181, 182 Conolly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
Becker . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 Consistance . . . . . . . . . . . . . . . . . . . . . 159, 164
Bellazi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 Convergence . . . . . . . . 14, 53, 55, 58, 59, 74,
Belt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 77, 78, 81–85, 108–110, 117–120,
Berman . . . . . . . . . . . . . . . . . . . . . . . . . . . 21, 72 137, 142
Convergente . . . . . . . . . . . . . . . . . . . 51, 75, 82
Bernardelli . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Conway . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
Bernouilli . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
Cooke . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
Beverton-Holt . . . . . . . . . . . . . . . . . . . . . . . 107
Cramer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
Bolie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
Creuse . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37, 47
Bolker . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175
Bonnet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 D
Boutayeb . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 Décomposition . . . 44, 46, 75, 142, 143, 195
Brauer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Définie positive . . . . . . . . . . . . . . . . . . . . . . . 74
Broyden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 Définie positive 9, 34, 41, 46, 58, 59, 74, 75,
broyden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189 78, 123
Budan-Fourier . . . . . . . . . . . . . . . . . . . . . . . . 97 Danilevski . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
223
Della . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 Gauss . 26, 29–33, 40, 47, 48, 51, 55–57, 60,
Dengue . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176 62, 74, 76–79, 190
Derouich . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 Givens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
Descartes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 Globale . . . . . . . . . . . . . . . . . . . . . . . 82, 83, 118
Descente . . . . . . . . . . . . . . . . 25, 120, 122, 123 Glucose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180
Diabète . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 Golub . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
diagonalisable . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Gradient . . . . . 121, 123, 124, 134, 192, 194
Diekman . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 Graeffe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
Dietz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 Grenfell . . . . . . . . . . . . . . . . . . . . . . . . 175, 182
Directe . . . . . . . . . . . . . . . . . . . . . . . 24, 50, 137 Grover . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
Direction . . . . . . . . . . . . . . . . . . . . . . . 121, 123
Dobson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 H
H-matrice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
E Hahn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
Ecologiques . . . . . . . . . . . . . . . . . . . . . . . . . . 125 Haner . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
Elémentaire . . . . . . . . . . 26, 27, 29, 31, 32, 98 Hassel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
Epidémiologiques . . . . . . . . . . . . . . . . . . . . 169 Hassell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
Equivalente . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 hermitienne . 7, 9, 11, 46, 52, 58, 59, 74, 75,
Ergodicité . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 78
Erreur . . . . 36, 37, 51, 81, 82, 109, 110, 113 Hessenberg . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
Euler . . . . . . . . . 106, 161–163, 165, 168, 184 Hessien . . . . . . . . . . . . . . . . . . . . . . . . . 123, 124
Evans . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Hethcote . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
Exemple . . . . . . . . . . . . . . . . . . . 12, 13, 27, 30, Heun . . . . . . . . . . . . . . . . . . . . . . . . . . . 165, 166
35, 36, 44, 82–84, 86, 97, 98, 100, Hilbert . . . . . . . . . . . . . . . . . . . . . . . . 13, 35, 37
116, 118, 121, 122, 138, 141 Himsworth . . . . . . . . . . . . . . . . . . . . . . . . . . 179
Exercices . . . . . . . . . . . . . . . . . . 22, 48, 73, 133 Hirsch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
Holling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
F Hoppensteadt . . . . . . . . . . . . . . . . . . . . . . . . 182
Factorisation . . . . . . . . . . . . 32–36, 39–42, 47 Horn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Fehlberg . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 Horner . . . . . . . . . . . . . . . . . . . . . . . . . . 91, 191
Fibonacci . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 Householder . . . . 21, 26, 35, 36, 49, 72, 191
Frankel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Hudson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
Frauenthal . . . . . . . . . . . . . . . . . . . . . . . . . . 182
Frobenis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 I
Frobenius . . . . . . . . . 14, 21, 26, 72, 139, 140 Imprimitive . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
Indirecte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
G Induite . . . . . . . . . . . . . . . . . . . . . 12, 51, 52, 58
Gantmakher . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Inférieure . . . . . . . 22, 25, 33, 34, 39, 77, 142
224
Insuline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 Longueur . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
Inversible . . 8, 12, 13, 24, 26, 30–34, 36, 37, Lotka . . . . . . . . . . . . . . . . . . . . . . . . . . . 132, 181
39–41, 50, 53, 56, 58, 75–78, 117, LR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
119, 139, 142 LU . . . . . . . . . . . . . . . 32–34, 47, 75, 191, 195
Irréductible . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Ludwig . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
irréductible . . . . . . . . . . . . . . . . . . . . . . . . 16, 39 Ludwing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
Isham . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
Isolé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 M
Itération . . . . . 51, 53, 73, 74, 81, 83, 84, 115 M-matrice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
Méthode . . . . . . . . . . . . . . . . . . . . . . 24, 25, 29–
J 33, 36, 40, 45–47, 50–53, 55–57,
Jacobi . 50, 53, 60–63, 74–78, 145–147, 190 59–63, 74, 75, 77, 78, 86, 91, 115,
Johnson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 117–124, 137–139, 142, 145–147
Jones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 Malthus . . . . . . . . . . . . . . . . . . . . . . . . 102, 106
Jordan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Marcus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . 14, 15
K
Kahan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 Maximum . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
Kantorovich . . . . . . . . . . . . . . . . . . . . . . . . . 119 May . . . . . . . . . . . . . . . . . . . . . . . 107, 132, 182
Karlin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Maynard . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
Keller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 McKendrick . . . . . . . . . . . . . . . . . . . . . . . . . 181
Ker . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 Medeley . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
Kermack . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181 Mehmke . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
Kot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132, 182 Meurant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
Kranz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 Mewborn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Krylov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 Minc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Minimum . . . . . . . . . . . . . . . . 45, 78, 120–122
L Modèles . . . . . . . . . . . . . . . . . . . . . . . . 125, 169
Lambert . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 Moindres carrées . . . . . . . . . . . . . . . . . . . . . . 45
Lauwerier . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 Mollison . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
Lawton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 Morison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
Leonard . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 Multi-pas . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
Leslie . . . . . . . . . . . . . . . . . . . . . . . . . . 15, 16, 21
Liapunov . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 N
Linéaire . . . . . . . . . . . . . . . . . . . . . . . 24, 47, 50 Newton . 85, 86, 96, 99, 109, 113, 116–119,
Lipschitz . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 122, 133, 196
Locale . . . . . . . . . . . . . . . . . . . . . . . . 84, 85, 118 Nicholson . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
Logistique . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 Non carrée . . . . . . . . . . . . . . . . . . . . . . . . 43, 44
Logofet . . . . . . . . . . . . . . . . . . . . . . . . . . 21, 132 Normale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
225
Norme . 6, 11, 12, 22, 36, 51, 52, 58, 75, 77, Rectangulaire . . . . . . . . . . . . . . . . . . 43–46, 48
115, 122 Regula falsi . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
Numérique . . . . . . . . 24, 37, 81, 91, 101, 121 Relaxation . . 51, 57–59, 61–63, 75, 78, 118,
193
O Ricker . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
Okubo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 Rosenberg . . . . . . . . . . . . . . . . . . . . . . . . . 71, 72
Oliveira . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 Ross . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
Opération . . . . . . . . . . . . . . . . . . . . . 24, 25, 47 Rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
Optimisation . . . . . . . . . . . . . . . . . . . . . . . . 119 Runge-Kutta . . . . . . . . . . . 163, 164, 166, 167
Orthogonale . . . . 7, 26, 35–37, 46, 145, 147
Ostrowski . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 S
Schaefer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
P Schur . . . . . . . . . . . . . . . . . . . . . . . . . . 9, 10, 12
Parker . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 Scott . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
Perlis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 Scuda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
Permutation . . . . . . . . . . . . 7, 27, 31–33, 146 SEIR . . . . . . . . . . . . . . . . . . . . . . . . . . . 170, 174
Perron . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14, 21 SEIRS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
Perturbation . . . . . . . . . . . . . . . . . . . . . . 36, 37 Semblable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
Pielou . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 Serge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
Pivot . . . . . . . . . . . . . . . . . . . . . . . . . . 29–33, 47 Sherman . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
Plemmons. . . . . . . . . . . . . . . . . . . . . . . . . 21, 72 SI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
Point fixe . . . . . . . . . . . . . . . 82, 115, 116, 194 Singulière . . . . . . . . . . . . 9, 22, 36, 44, 46, 48
Polynôme8, 60, 81, 91, 94, 96, 97, 100, 101, SIR . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170, 172
138, 139
SIS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
Primitive . . . . . . . . . . . . . . . . . . . . . . . . . . . 7, 14
Smilatova . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
Produit scalaire . . . . . . . . . . . . . . . . . . . . . . . . 6
Smith . . . . . . . . . . . . . . . . . . . . . . . . . . 136, 182
Puissance . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
Sous-matrice . . . . . . . . . . . . . . . . . . . . . . 33, 34
Southwood . . . . . . . . . . . . . . . . . . . . . . . . . . 109
Q
Spectre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
QR . . . . . . . . . . . . . . . . . . . . . . . . . 47, 143, 144
Stable . . . . . . . . . . . . . . . . . . . . . . . . . . . 36, 137
Quasi-dominante . . . . . . . . . . . . . . . . . . . . . . . 8
Stationnaire . . . . . . . . . . . . . . . . . . . . . . . . . 156
R Steffenson . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
Réductible . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Stein . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71, 72
Résolution 24, 25, 30, 36, 40, 41, 44, 46, 50, Stetter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
63, 122, 124, 139 Stieljes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7, 72
Racine . 8, 9, 60, 62, 64, 65, 78, 86, 94–101, Stochastique . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
110, 113, 138 Structure . . . . . . . . . . . . . . . . . . . . . . . . . 37, 46
Rayon spectral . . . . . 6, 8, 22, 52, 59, 74, 137 Sturm . . . . . . . . . . . . . . . . . . 99, 100, 114, 193
226
Subordonnée . . . . . . . . . . . . . . . . . . . 12, 22, 36 Y
Sujan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 Young . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
Supérieure 14, 22, 25, 29, 31, 33, 39, 77, 96, Yule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
142, 143
Z
Symétrique . . . . . . . 7, 26, 34, 46, 74, 98, 145
Ziegler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
Système . . 14, 24, 25, 29–31, 34, 36, 38, 40,
44–47, 49, 50, 73, 74, 77, 78, 116,
124, 137, 139
T
Taylor . . . . . . . . . . . . . . . 81, 96, 120, 162–164
Triangulaire . . 22, 25, 29, 31–34, 38, 39, 46,
76, 77, 137, 142, 143
Tridiagonale . . 39, 49, 60, 61, 73, 75, 77, 78,
113
Troncature . . . . . . . . . . . . . . . . . 159, 167, 168
U
Unitaire . . . . . . . . . . . . . . . . . . . . . 7, 35, 43, 48
V
Valeur propre . . . . . . . . . . . . . . . . . . . . 7–9, 11,
14, 22, 23, 60–62, 74, 76–79, 114,
137, 138, 141–143, 147
Van Loan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
Varga . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21, 72
Vecteur propre . . . . 8, 14, 137, 139, 141, 147
Verhulst . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
Verner . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
Vitesse . . . . . . 14, 53, 58, 74, 77, 78, 81, 108
Volterra . . . . . . . . . . . . . . . . . . . . . . . . 131, 132
von Mises . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
W
Wegstein . . . . . . . . . . . . . . . . . . . . . . . . . . 84, 85
White . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
Whittaker . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
Wickwire . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
227