Introduction

L’analyse numérique est un domaine des mathématiques appliquées plus connue actuel-
lement sur le nom de modélisation et calcul scientifique.
Dans ce domaine on s’intéresse à la mise en équation de problèmes concrêts et à la définition
de méthode numérique puisssantes pour la résolution.
Pour la mise en équation ( modélisation) le mathématicien a besoin des spécialistes du do-
maine où intervient le problème ( biologie, mécanique, transport, météo, etc..).
1) la modélisation
Grâce aux théories et aux sciences existantes, on arrive à mettre en équation un problème
donné. c’est une étape trés importante, en effet une modélisation hasardeuse d’un problème
peut n’avoir aucun lien avec le problème réel à résoudre.
2)Existence
Le modèle mathématique obtenu est généralement trés complexe. l’existence et éventuellement
l’unicité de la solution sont difficiles à montrer. cela constitue une première tâche du mathématicien.
3)Approximation
L’existence et l’unicité sont montrées, on passe (procède) alors à l’approximation de la so-
lution par des méthodes de discrétisation. les méthodes les plus utilisées (Différences finies,
éléments fins, volumes finis, transformations de Fouriers, ondelettes, méthodes spectrales,etc.)
4)Calcul
Aprés la discrétisation du problème on l’écrit sous forme de progamme dans un langage donné
( C,C++, Java, Python, Matlab,Scilab,...)
5) Interprétation
En collaborationavec les spécialistes de l’origine du problème traité, on procède à l’interprétation
des résultats numériques obtenus et éventuellement quand cela est possible à la comparaison
avec des essais expérimentaux.

1
Université de Thiès Dpt Maths

Draft,I.MBAYE 2 2018-2019
Chapitre 1

Equations algébriques

1.1 Le problème
Soit une fonction réelle continue d’une variable réelle f : I = [a, b] → R telle que f (x) = 0.
l’objectif est de trouver la racine α de l’équation
(E) : f (x) = 0.
Certaines situations peuvent se produire ( pas de solution, une solution unique ou plusieurs
solutions)

1.2 Solution itérative

*) pour une fonction non linéaire quelconque la solution analytique est rarement dispo-
nible.
*) les méthodes numériques pour approcher une racine α sont en générales itératives
*) si la solution est itérative le problème consiste à construire une suite (xk )k∈N telle que
limk→∞ xk = α
*) En notant g0 = f le problème est lié à la recherche des extréma d’une fonction g
(optimisation) car les zéro de f sont les extréma de g.

1.3 Critères d’arrrêt des itérations

Pour faire cesser le calcul des termes successifs de xk destinés à approcher la limite α
solutio de l’équation (E). Nous devons choisir un critère qui entraine l’arrêt des itérations :
celui-ci provoquera souvent une sortie de boucle dans les algorithmes utilisés. Les critères
d’arrêt peuvent être trés divers. la liste présentée est loin d’être exhaustive :
*) Majoration par une valeur seuil  de l’écart |xk+1 − xk | entre les termes consécutifs de
la suite xk .
*) Majoration par une valeur seuil  de la variation relative |xk+1 −xk |
|xk+1 |
des termes de la
suite xk .

3
Université de Thiès Dpt Maths

*) La majoration de |f (xk )| par une valeur seuil 

*) La majoration du nombre d’itération par une valeur N agissant comme un garde
fou destiner à éviter les temps de calcul excessif et la perte de contrôle probable de
l’évolution des erreurs d’arrondi qui risque d’en découler.

1.4 Objectifs
Nous présentons les méthodes les plus classiques :
*) La méthode des approximations succcessives( point fixe)
*) La méthode de Newton
*) la méthode de dichotomie ou bissection
‘

1.5 La méthode des approximations succcessives

Soit à résoudre l’équation (E) : f (x) = 0 où f ∈ C 0 ([a, b]). On suppose que l’équation

f (x) = 0 ⇔ x = g(x).

La fonction g devant vérifier certaines hypothèses (x = g(x) est dit problème du point fixe).
A partir de x0 ∈ [a, b] on construit la suite xk par xk+1 = g(xk ). cet algorithme est dit
méthode des approximations successives (figure).
Théorème 1 : Supposons que g est définie et continue sur I = [a, b] et que :
a) g(I) ⊂ I
b) g est contractante i.e

∃L ∈ [0, 1[, |g(x) − g(y)| ≤ L|x − y|, ∀(x, y) ∈ I 2 .

‘ Alors ∀x0 ∈ I, la suite xk converge vers α unique solution de l’équation x = g(x).

Preuve : montrons que la suite xk est de cauchy.
soit m > n deux entiers naturels. Alors
m−1
X m−1
X
|xn − xm | = | xk − xk+1 | ≤ |xk − xk+1 |.
k=n k=n

On utilise le fait que xk ∈ I, ∀k et

|xk − xk+1 | = |g(xk−1 ) − g(xk )| ≤ L|xk−1 − xk | . . . ≤ Lk |x0 − x1 |

, d’où
1 − Lm−n
|xn − xm | ≤ |x0 − x1 |Ln →0
1−L
lorsque n → ∞. Donc xk est une suite de cauchy, d’où la convergence vers α ∈ I. Comme g
est continue on a : g(α) = α.

Draft,I.MBAYE 4 2018-2019
Université de Thiès Dpt Maths

Remarque : α ainsi trouvé est unique solution de l’équation g(x) = x. En effet s’il existe
deux solutions α1 et α2 on aurait

|α1 − α2 | = |g(α1 ) − g(α2 )| ≤ L|α1 − α2 |

ce qui implique :
|α1 − α2 | < |α1 − α2 |
car L < 1 et |α1 − α2 | =6 0 ce qui est absurde.
Remarque : Une condition suffisante pour que la fonction g soit contractante est que g
soit dérivable sur I et que supx∈I |g 0 (x)| = L < 1. pour la preuve on applique le théorème des
accroissements finis sur I.
Théorème 2 : On suppose que les hypothèses du théorème 1 sont vérifiées. Alors on a :

Ln
|xn − α| ≤ |x0 − x1 | .
1−L

preuve d’aprés la preuve du théorème 1 on a

1 − Lm−n
|xn − xm | ≤ |x0 − x1 |Ln
1−L
passons à la limite à l’infini m → ∞ dans l’inégalité. On obtient alors la majoration d’ erreur
suivante
Ln
|xn − α| ≤ |x0 − x1 | .
1−L

Définition : On dit que la convergence de la suite xk vers α est d’ordre p ∈ N∗ si et seulement

si ∃ap > 0 tel que
|xn+1 − α|
lim = ap .
n→∞ |xn − α|p

Théorème 3 : la méthode des approximations successives est d’ordre 1 dans le cas général.
Preuve : On fait un DL à l’ordre 1 autour de α de g(xk ) ; on a

1 00
xk+1 = g(xk ) = g(α) + g 0 (α)(xk − α) + g (yk )(xk − α)2
2!
où α < yk < xk . on a alors
xk+1 − α
= g 0 (α) + g 00 (yk )(xk − α)
xk − α

ce qui implique :
|xk+1 − α|
= |g 0 (α) + g 00 (yk )(xk − α)|
|xk − α|

Draft,I.MBAYE 5 2018-2019
Université de Thiès Dpt Maths

lorsque k tend vers l’infini on obtient :

|xk+1 − α|
lim p
= |g 0 (α)|.
k→∞ |xk − α|

dans le cas général. On déduit que la convergence de la méthode des approximations succes-
sives est d’ordre 1.
Remarque : (ordre de convergence) soit g(x) = x ⇔ f (x) = 0. Si g (1) (s) = g (2) (s) = . . . =
g (r) (s) = 0 alors la méthode est au moins d’ordre r + 1.

1.6 Méthode de Newton

Soit à résoudre l’équation (E) : f (x) = 0 où f ∈ C 0 ([a, b]). ON suppose que l’équation
f (x) = 0 ⇔ x = g(x). où
f (x)
g(x) = x − 0 .
f (x)
On pose
f (xk )
xk+1 = xk −
f 0 (xk )
avec x0 donné, cet algorithme est dit méthode de Newton.
Théorème 4 : On suppose que f s’annule en α. On suppose que dans un voisinage de α,
la fonction f est classe C 2 et que f 0 (x) 6= 0 pour tout x ∈ Vα . Alors la méthode de Newton
converge vers α et la convergence est au moins d’ordre 2.
Preuve : puis que g 0 (α) = 0, il existe donc un voisinage Vα tel que |g 0 (α)| < 1. En appliquant
le théorème 1 on a donc la convergence dés que x0 ∈ Vα .
En outre, on fait un DL de f (α) autour de xn on a alors :

1
0 = f (α) = f (xn ) + (α − xn )f 0 (xn ) + (α − xn )2 f 00 (yn )
2!
d’où
xn+1 − α f 00 (yn )
= ,
(xn − α)2 2f 0 (xn )
on a finalement :
|xn+1 − α| f 00 (α)
lim = | |.
n→∞ |xn − α|2 2f 0 (α)
Donc la méthode est au moins d’ordre 2.
Théorème 5 : (Estimation d’erreur) On suppose que les hypothèses du théorème 4 sont
satisfaites. On pose
M2 = sup |f 00 (x)|,
x∈I

m1 = inf |f 0 (x)|
x∈I

Draft,I.MBAYE 6 2018-2019
Université de Thiès Dpt Maths

et
M2
C=
2m1
On a alors l’estimation d’erreur :
1 n
|xn − α| ≤ [C|x0 − α|]2 ,
C
∀n ∈ N.
textbfPreuve : D’aprés le théorème 4 On a : |xn+1 − α| ≤ C|xn − α|2 , ∀n ∈ N. Nous allons
démontrer par récurrence le résultat suivant :
1 n
|xn − α| ≤ [C|x0 − α|]2 ,
C
. Pour n = 0 on a : |x0 − α| ≤ |x0 − α| vrai. Supposons que la relation est vraie au rang n et
montrons qu’elle est vraie au rang n + 1 i.e
1 n+1
|xn+1 − α| ≤ [C|x0 − α|]2 ,
C
On a
|xn+1 − α| ≤ C|xn − α|2 ,
d’aprés l’hypothèse de récurrence on a alors, le résultat.
Exercice : soit f (x) = x3 + x − 1000 = 0 sur I = [9, 10]
A) x = 1000 − x3 ,
1000
B) x = x2
− x1 ,
1
C) x = (1000 − x) 3
Parmi les méthodes A, B, C laquelle converge ?
Sous matlab on peut utiliser la fonction fzero() ou fsolve() pour trouver la racine approchée
de l’équation f (x) = 0 sur [a, b]. En effet on écrit les lignes de commandes suivantes :
>> f=inline(’x. ∧ 3 + x − 1000’)
>> r=fzero(f,9)

on obtient la racine r = 9.9667 de l’équation non linéaire f (x) = x3 + x − 1000 = 0.

Exercice :
réaliser une telle étude graphique et établir la convergence ou la divergence de la suite dans
les cas suivants :
√
1) un+1 = 1 + un avec u0 = 1 ,
p
2) un+1 = 1 + u2n avec u0 = 1,
3
3) un+1 = 1+2u2n
avec u0 = 0.5,
2
4) un+1 = 1+u2n
avec u0 = 0.5.

Draft,I.MBAYE 7 2018-2019
Université de Thiès Dpt Maths

1.7 Méthode de Dichotomie

On suppose que f (x) = 0 admet une solution unique α ∈ I et que f ∈ C 0 (I). l’idée de la
dichotomie est de localiser la solution α dans des intervalles dont la taille se divise par deux
à chaque étape. il est claire que la convergence est assurée par la définition de la méthode.
La localisation de la solution utilise le théorème des valeurs intermédiaires. On procède de la
manière suivante :
*) On a : f (a) × f (b) < 0,
*) on pose c = a+b 2
,
*) si f (a) × f (c) < 0 alors α ∈]a, c[,
*) si f (a) × f (c) = 0 alors α = a ou α = c,
*) si f (c) × f (b) < 0 alors α ∈]c, b[,
*) on réitère le procédé jusqu’à atteindre la précision souhaitée.
‘

1.8 Notion d’accélération de la convergence

Définition : Soit xn et zn deux suites qui converge vers α. On dit que zn converge plus
vite que xn si
zn − α
lim | | = 0.
n→∞ xn − α

Théorème 6 : Si xn est d’ordre p et zn d’ordre r tel que r > p alors zn converge plus vite
que xn .

Draft,I.MBAYE 8 2018-2019
Chapitre 2

Interpolation Polynômiale

2.1 Interpolation de Lagrange

Soient x0 , x1 . . . , xn n + 1, points de la droite réelle distincts deux à deux et y0 , y1 . . . , yn ,
n + 1 autres points de la droite réelle. L’interpolation consiste à trouver une fonction f
telle que f (xi ) = yi pour tout i ∈ {0, 1, . . . , n}. En outre si f est un polynôme on parle
d’interpolation polynômiale.
Proposition 1 :Soient x0 , x1 . . . , xn n + 1, points de la droite réelle distincts deux à deux
et y0 , y1 . . . , yn , n + 1 autres points de la droite réelle. Il existe un unique polynôme P de
degré inférieur où égal à n tel que P (xi ) = yi . Ce polynôme ainsi trouvé est appelé polynôme
d’interpolation aux points xi pour i ∈ {0, 1, . . . , n}.
Preuve (méthode directe) :
Soit P (x) = a0 + a1 x + . . . + an xn un polynôme de degré n tel que  P (xi ) = yi pour tout 
1 x0 . . . xn0
 1 x1 . . . xn 
1 
i ∈ {0, 1, . . . , n}, On a le système d’équation suivant : AX = b où A =  .. .. ..  .,

..
 . . . . 
1 xn . . . xnn
 
a0
 a1  Q
X =  ..  et , Ce système admet une solution unique si et seulement si det(A) = i6=j (xi −
 
 . 
an
xj ) 6= 0 (Déterminant de Vendermone), ce qui est le cas car les xi sont distincts deux à deux
pour i ∈ {0, 1, . . . , n}. Unicité : soient P1 et P2 deux polynômes distincts tels que P1 (xi ) = yi
et P2 (xi ) = yi , on a alors P1 (xi ) = P2 (xi ) ce qui implique que (P1 − P2 )(xi ) = 0 d’où P1 − P2
est un polynôme de degré inférieur ou égal à n qui a n + 1 racines ce qui est absurde donc
P1 = P2 .
Preuve (méthode constructive) :
On définit les polynômes de lagrange li pour tout i vérifiant li (xj ) = δij on obtient alors
(x − x0 ) . . . (x − xi−1 )(x − xi+1 ) . . . (x − xn )
li (x) = .
(xi − x0 ) . . . (xi − xi−1 )(xi − xi+1 ) . . . (xi − xn )

9
Université de Thiès Dpt Maths

On définit le polynôme suivant

n
X
P (x) = yi li (x)
i=0

qui vérifie P (xi ) = yi pour tout i et qui est unique. Ce polynôme ainsi défini est appelé
polynôme d’interpolation de Lagrange aux points xi .
Exercice : Déterminer le polynôme d’interpolation de Lagrange P aux points x0 = −1, x1 =
0, x2 = 1 tels que P (x0 ) = 2, P (x1 ) = 1, P (x2 ) = 3.
1
Exercice : Soit la fonction f (x) = 1+x 2 . Déterminer le polynôme d’interpolation de Lagrange

P interpolant f points x0 = −1, x1 = 0, x2 = 1.

Théorème 1 : Soient P un polynôme interpolant f aux points xi pour i ∈ {0, 1, . . . , n},
avec f ∈ C n+1 ([a, b]). Alors on a l’erreur d’interpolation définie par pourQtout x ∈ [a, b], il
existe ψ ∈]a, b[ tel que f (x) − P (x) = (n+1)!1
f (n+1) (ψ)W (x), avec W (x) = ni=0 (x − xi ).
Preuve

2.2 Interpolation d’Hermite

Proposition 1 : Soient x0 , x1 . . . , xn n + 1, points de la droite réelle distincts deux à
deux, y0 , y1 . . . , yn , n + 1 autres points de la droite réelle et z0 , z1 . . . , zn , n + 1 autres points
de la droite réelle. Il existe un unique polynôme P de degré inférieur où égal à 2n + 1 tel
que P (xi ) = yi et P 0 (xi ) = zi . Ce polynôme ainsi trouvé est appelé polynôme d’interpolation
d’Hermine aux points xi pour i ∈ {0, 1, . . . , n}.
Preuve (méthode constructive) : le polynôme s’écrit :
n
X n
X
P (x) = ψi (x)yi + φi (x)zi ,
i=0 i=0

où
ψi (x) = [1 − 2(x − xi )li0 (x)]li2 (x)
et
φi (x) = (x − xi )li2 (x).
Exercice : Déterminer le polynôme d’interpolation d’Hermite P aux points x0 = −1, x1 =
0, x2 = 1 tels que P (x0 ) = 0, P (x1 ) = 1, P (x2 ) = 3 et P 0 (x0 ) = 2, P 0 (x1 ) = 1, P 0 (x2 ) = 3 .
1
Exercice : Soit la fonction f (x) = 1+x 2 . Déterminer le polynôme d’interpolation d’Hermite

P interpolant f points x0 = −1, x1 = 0, x2 = 1.

Théorème 2 : Soient P un polynôme interpolant f aux points xi pour i ∈ {0, 1, . . . , n},
avec f ∈ C 2n+2 ([a, b]). Alors on a l’erreur d’interpolation définie par pour tout Q x ∈ [a, b], il
1
existe ψ ∈]a, b[ tel que : f (x) − P (x) = (2n+2)! f (2n+2) (ψ)W 2 (x), avec W (x) = ni=0 (x − xi ).
Preuve

Draft,I.MBAYE 10 2018-2019
Chapitre 3

Intégration numérique ou Quadratures

Rb
Soit f ∈ C 0 ([a, b]). On cherche la valeur de l’intégrale I = a
f (x)dx. Pour cela on utilise
la formule d’approximation suivante :
Z b n
X
f (x)dx ≈ wi f (xi )
a i=0

où les wi sont les poids de quadratures et les xi sont les noeuds. Différentes méthodes seront
étudiées.

3.1 Formules de quadratures

3.1.1 Méthode du rectangle
Rb
Soit f ∈ C 0 ([a, b]), on cherche la valeur de l’intégrale I = a f (x)dx par la méthode du
rectangle à gauche, à droite et du point milieu respectivement :
Z b
f (x)dx ≈ f (a)(b − a),
a

Z b
f (x)dx ≈ f (b)(b − a)
a

et
Z b
a+b
f (x)dx ≈ f ( )(b − a)
a 2
.
les méthodes du rectangle sont à un point.

11
Université de Thiès Dpt Maths

3.1.2 Méthode des trapèzes

Rb
Soit f ∈ C 0 ([a, b]), on cherche la valeur de l’intégrale I = a f (x)dx par la méthode des
trapèzes définie par : Z b
(b − a)
f (x)dx ≈ (f (a) + f (b)) .
a 2

La méthode des trapèzes est à deux points a et b.

3.1.3 Méthode de Simpson

Rb
Soit f ∈ C 0 ([a, b]), on cherche la valeur de l’intégrale I = a f (x)dx par la méthode de
Simpson définie par :
Z b
(b − a) a+b
f (x)dx ≈ (f (a) + 4f ( ) + f (b)).
a 6 2

La méthode de Simpson est à trois points a, a+b2R

et b.
1 x
Exercice : Calculer la valeur de l’intégrale I = 0 e dx avec la méthode du rectangle (gauche
et droite), la méthode du point milieu, la méthode des trapèzes et la méthode de Simpson.
Comparer les résultats avec la valeur exacte de I. conclure

3.1.4 Définition et Propositions

Définition 1 : (Ordre d’une méthode) Une méthode est au moins d’ordre n si pour tout
polynôme P de degré inférieur ou égal à n l’erreur de quadratures
Z b n
X
E(P ) = P (x)dx − wi f (xi ) = 0.
a i=0

Proposition 1 : La méthode du rectangle gauche et à droite sont d’ordre R b 0.

En effet : - soit f (x) = P (x) = 1 on a alors Iq = (b − a), d’autre part a P (x)dx = (b − a),
d’où E(P ) = 0 ;
Rb 2 2)
-soit f (x) = P (x) = x on a alors Iq = a(b − a), d’autre part a P (x)dx = (b −a 2
, d’où
E(P ) 6= 0.
Donc la méthode est d’ordre 0.
idem pour la méthode des rectangles à droite.
Proposition 2 : La méthode du poin milieu gauche est d’ordre 1. Rb
En effet : - soit P (x) = 1 on a alors Iq = (b − a), d’autre part a P (x)dx = (b − a), d’où
E(P ) = 0 ;
2 2) Rb 2 2)
-soit P (x) = x on a alors Iq = a+b2
(b − a) = (b −a2
, d’autre part a P (x)dx = (b −a 2
, d’où
E(P ) = 0.
2 (b2 −a2 ) Rb
- soit f (x) = P (x) = x2 on a alors Iq = (a+b)4
(b − a) = 2
, d’autre part a
P (x)dx =

Draft,I.MBAYE 12 2018-2019
Université de Thiès Dpt Maths

(b3 −a3 )
3
d’où E(P ) 6= 0
,
Donc la méthode est d’ordre 1.
Proposition 3 : La méthode des trapèzes est d’ordre 1. Rb
En effet : - soit f (x) = P (x) = 1 on a alors Iq = (b − a), d’autre part a P (x)dx = (b − a),
d’où E(P ) = 0 ;
2 2) Rb 2 2)
-soit f (x) = P (x) = x on a alors Iq = a+b
2
(b − a) = (b −a
2
, d’autre part a P (x)dx = (b −a
2
,
d’où E(P ) = 0.
Rb 3 3)
- soit f (x) = P (x) = x2 on a alors Iq = f raca2 + b2 2(b − a), d’autre part a P (x)dx = (b −a3
,
d’où E(P ) 6= 0
Donc la méthode est d’ordre 1.
Proposition 4 : La méthode de Simpson est d’ordre 3.
Preuve (voir proposition 1,2,3)

3.1.5 Erreur d’intégration

Proposition 1 : Soit f ∈ C 1 ([a, b]). La méthode du rectangle gauche (à droite ) a pour
erreur Z b
(b − a)2
Eg = f (x)dx − f (a)(b − a) = f 0 (η)
a 2
et Z b
(b − a)2
Ed = f (x)dx − f (b)(b − a) = −f 0 (η) .
a 2

En effet, en utilisant le DL de f autour de a à l’ordre 1, on a f (x) = f (a) + (x − a)f 0 (ηx )

Rb Rb Rb
puis en intégrant sur [a, b] on a a f (x)dx = a f (a)dx + a (x − a)f 0 (ηx )dx et en utilisant
Rb Rb
la formule de la moyenne on obtient le résultat a f (x)dx = f (a)(b − a) + f 0 (η) a (x − a)dx
Rb 2
On en déduit l’estimation d’erreur suivante pour Eg = a f (x)dx − f (a)(b − a) = f 0 (η) (b−a)
2
la méthode des rectangle à gauche. Pour la méthode du rectangle à droite on fait le même
raisonnement autour de b.
Proposition 2 : Soit f ∈ C 2 ([a, b]). La méthode du point mileu a pour erreur d’intégration
Z b
a+b 00 (b − a)3
Em = f (x)dx − f ( )(b − a) = f (η) .
a 2 24

En effet, en utilisant le DL de f autour de a+b 2

à l’ordre 2, on a f (x) = f ( a+b 2
) + (x −
a+b 0 a+b 1 a+b 2 00
2
)f ( 2
) + 2!
(x − 2
) f (ηx ). puis en intégrant sur [a, b] et en utilisant la formule de la
moyenne on obtient le résultat.
Proposition 3 : Soit f ∈ C 2 ([a, b]). La méthode des trapèzes a pour erreur d’intégration
Z b
f (a) + f (b) (b − a)3
Et = f (x)dx − ( )(b − a) = −f 00 (η)
a 2 12
.
En effet en utilisant le théorème 1 du chapitre 2 pour n = 1 on a alors f (x) − P (x) =

Draft,I.MBAYE 13 2018-2019
Université de Thiès Dpt Maths

1 00
Rb Rb Rb
2!
f (ψ)(x − a)(x − b). En intégrant sur [a, b] on a : a f (x)dx − a P (x)dx = 2!1 a f 00 (ψ)(x −
Rb
a)(x − b)dx et en utilisant la formule de la moyenne on obtient le résultat : a f (x)dx −
Rb Rb
( f (a)+f
2
(b)
)(b − a) = 2!1 f 00 (η) a (x − a)(x − b)dx. On en déduit l’erreur : Et = a f (x)dx −
3
( f (a)+f
2
(b)
)(b − a) = −f 00 (η) (b−a)
12
.
4
Proposition 4 : Soit f ∈ C ([a, b]). La méthode de Simpson a pour erreur d’intégration
b
(b − a) (b − a)5
Z
a+b
Es = f (x)dx − (f (a) + 4f ( ) + f (b)) = −f (4) (η) .
a 6 2 2880

En effet, en utilisant le théorème 1 du chapitre 2 pour n = 2 on a alors f (x) − P (x) =

1 (4)
4!
f (ψ)(x − a)(x − a+b 2
)2 (x − b). en intégrant sur [a, b] et en utilisant la formule de la
moyenne on obtient le résultat :

3.2 Formules composites ou répétées

On subdivise l’intervalle [a, b] en N intervalle de longueur h = b−a
N
. On pose xi = a + i ∗ h
pour tout i ∈ {0, 1, . . . , N } appelés points du maillage ou de la discrétisation. On a aussi
xi+1 − xi = h, x0 = a et xN = b.

3.2.1 Rectangle composite

Sur chaque intervalle [xi , xi+1 ] pour i ∈ {0, 1, . . . , N − 1}, on applique la formule du
rectanagle, ensuite on somme sur i de 0 à N − 1. On a alors :
Z xi+1
f (x)dx ≈ f (xi )(xi+1 − xi ),
xi

ce qui entraine la formule du rectangle à gauche composite suivante :

Z b N
X −1
f (x)dx ≈ h f (xi );
a i=0

de même pour la formule du rectangle à droite composite suivante :

Z b N
X −1
f (x)dx ≈ h f (xi+1 );
a i=0

de même pour la formule du point milieu composite suivante :

Z b N −1
X h
f (x)dx ≈ h f (xi + ).
a i=0
2

Draft,I.MBAYE 14 2018-2019
Université de Thiès Dpt Maths

3.2.2 Trapèze composite

Sur chaque intervalle [xi , xi+1 ] pour i ∈ {0, 1, . . . , N − 1}, on applique la formule des
trapèzes, ensuite on somme sur i de 0 à N − 1. On a alors
Z xi+1
f (xi ) + f (xi+1 )
f (x)dx ≈ (xi+1 − xi ).
xi 2
Ce qui entraine
Z b N
X −1
f (x)dx ≈ h(f (a) + 2 f (xi ) + f (b)).
a i=1

3.2.3 Simpson composite

Sur chaque intervalle [xi , xi+1 ] pour i ∈ {0, 1, . . . , N − 1}, on applique la formule des
trapèzes, ensuite on somme sur i de 0 à N − 1. On a alors
Z xi+1
(xi+1 − xi h
f (x)dx ≈ (f (xi ) + 4f (xi + ) + f (xi+1 ).
xi 6 2
Ce qui entraine
Z b N −1 N −1
h X X h
f (x)dx ≈ (f (a) + 2 f (xi ) + +4 f (xi + ) + f (b)).
a 6 i=1 i=0
2

3.2.4 Erreur d’intégration des formules composites

Proposition 1 : Soit f ∈ C 1 ([a, b]). La méthode du rectangle composite gauche (à droite
) a pour erreur
Z b
h
Eg = f (x)dx − f (a)(b − a) = (b − a) ∗ f 0 (η)
a 2
et Z b
h
Ed = f (x)dx − f (b)(b − a) = −(b − a) ∗ f 0 (η) .
a 2
2
En effet, Sur chaque intervalle [xi , xi+1 ] pour i ∈ {0, 1, . . . , N −1} on a Egi = f 0 (η) (h)2 , ensuite
2
on somme de i = 0 à i = N − 1, on obtient alors Eg = N ∗ f 0 (η) h2 ; et remplaçant N par
N = b−a h
, on obtient : Eg = (b − a)f 0 (η) h2
Remarques : |Eg | = O(h) et |Ed | = O(h).
Proposition 2 : Soit f ∈ C 2 ([a, b]). La méthode du point mileu composite a pour erreur
d’intégration Z b
a+b (h)2
Em = f (x)dx − f ( )(b − a) = (b − a)f 00 (η)
a 2 24
.
Preuve : (voir la preuve de la proposition 1).

Draft,I.MBAYE 15 2018-2019
Université de Thiès Dpt Maths

Remarque : |Em | = O(h2 ).

Proposition 3 : Soit f ∈ C 2 ([a, b]). La méthode des trapèzes composite a pour erreur
d’intégration
Z b
f (a) + f (b) (h)2
Et = f (x)dx − ( )(b − a) = −(b − a)f 00 (η) .
a 2 12

Preuve : (voir la preuve de la proposition 1).

Remarque : |Et | = O(h2 ).
Proposition 4 : Soit f ∈ C 4 ([a, b]). La méthode de Simpson a pour erreur d’intégration
b
(b − a) (h)4
Z
a+b
Es = f (x)dx − (f (a) + 4f ( ) + f (b)) = −(b − a)f (4) (η) .
a 6 2 2880

Preuve : (voir la preuve de la proposition 1).

Remarque : |Es | = O(h4 ).

3.2.5 convergence des formules composites

R b Proposition 1 : Les formules composites convergent vers la valeur I de l’intégrale
a
f (x)dx .
En effet, On a |E| = O(hp ) avec p ∈ N∗ , donc pour h → 0 i.e (N grand) |E| tend vers zéro
par conséquent Iq → I.
Calcul de la valeur approchée de l’intégrale sous matlab R1
On utilise la fonction quad(). Par exemple pour calculer I = 0 xex dx on écrit sous MatLab :
>> f=inline(’x.*exp(x)’)
>> I=quad(f,0,1)

on obtient le résultat suivant : I = 1.0000

Draft,I.MBAYE 16 2018-2019
Chapitre 4

Résolution numérique des équations

différentielles

4.1 Position du problème

Soit f ∈ C 0 ([t0 , t0 + T ] × R) à valeurs dans R. On cherche une fonction t 7→ y(t) ∈
C 1 ([t0 , t0 + T ]) et qui vérifie l’équation différentielle avec donnée initiale en t0 suivante :
 0
y (t) = f (t, y(t)), ∀t ∈ [t0 , t0 + T ]
(ED) :
y(t0 ) =η
Théorème 1 : On suppose que f est lipschitzienne par rapport à y : ∃L > 0 tel que
∀(y, y ∗ ) ∈ R, ∀t ∈ [t0 , t0 + T ] on ait : |f (t, y) − f (t, y ∗ )| ≤ L|y − y ∗ |. Alors il existe une unique
fonction t 7→ y(t) ∈ C 1 ([t0 , t0 + T ]) solution de l’(ED).
Preuve :
Exemple 1 : y 0 = t + y facile à résoudre et y 0 = exp(t) + sin(y) cette équation différentielle
admet une solution mais il est difficile à résoudre à la main. L’objectif dans ce chapitre est
de définir des méthodes d’approximation pour résoudre ce genre de problème. Dans tout
ce chapitre on considère les notations suivantes. On subdivise l’intervalle [t0 , t0 + T ] en N
T
intervalle de longueur h = N et on pose ti = t0 + ih où i ∈ {0, 1, . . . , N − 1}. On pose aussi
yi une approximation de y au point ti .

4.1.1 Méthode d’Euler

Sur chaque intervalle [ti , ti+1 ], on intègre l’(ED), on obtient alors
Z ti+1 Z ti+1
0
y (t)dt = f (t, y(t))dt.
ti ti
Rt
Pour retrouver la méthode d’Euler Explicite (progressive) on calcule tii+1 f (t, y(t))dt par la
Rt
méthode du rectangle à gauche on a alors tii+1 f (t, y(t))dt = hf (ti , y(ti )) . On définit alors

17
Université de Thiès Dpt Maths

la méthode d’Euler explicite par :


yi+1 = yi + hf (ti , yi ), ∀i ∈ {0, 1, . . . , N − 1}
y0 =η
Rt
Pour retrouver la méthode d’Euler implicite (rétrograde) on calcule tii+1 f (t, y(t))dt par la
Rt
méthode du rectangle à droite on a alors tii+1 f (t, y(t))dt = hf (ti+1 , y(ti+1 )) . On définit
alors la méthode d’Euler implicite par :

yi+1 = yi + hf (ti+1 , yi+1 ), ∀i ∈ {0, 1, . . . , N − 1}
y0 =η

4.1.2 Méthode du trapèze où de Crank-Nicolson

R ti+1
on calcule ti
f (t, y(t))dt par la méthode du trapèze on a alors
Z ti+1
h
f (t, y(t))dt = (f (ti , y(ti )) + f (ti+1 , y(ti+1 ))).
ti 2

On définit alors la méthode du trapèze par :

yi+1 = yi + h2 (f (ti , yi ) + f (ti+1 , yi+1 )), ∀i ∈ {0, 1, . . . , N − 1}



y0 =η

4.1.3 Méthode de Runge Kutta d’ordre 2

R ti+1
On calcule ti
f (t, y(t))dt par la méthode du point milieu on a alors
Z ti+1
ti+1 + ti ti+1 + ti
f (t, y(t))dt = h(f ( , y( )).
ti 2 2
Rt
On fait le DL autour de xi de y à l’ordre 1, on a alors tii+1 f (t, y(t))dt = h(f (ti + h2 , y(ti + h2 ))
or y(ti + h2 ) = y(ti ) + h2 y 0 (ti ) + . . .
On obtient la méthode de Runge Kutta d’ordre 2 :

yi+1 = yi + h(f (ti + h2 , yi + h2 f (ti , yi )), ∀i ∈ {0, 1, . . . , N − 1}



y0 =η

On pose yi+ 1 = yi + h2 f (ti , yi ) On obtient la schéma équivalent suivant :

2


 y0 =η
yi+ 1 = yi + h2 f (ti , yi )
2
yi+1 = yi + h(f (ti + h2 , yi+ 1 ), ∀i ∈ {0, 1, . . . , N − 1}

2

Draft,I.MBAYE 18 2018-2019
Université de Thiès Dpt Maths

4.1.4 Méthode de Runge Kutta d’ordre 2

Elle est basée sur l’utilisation de la méthode de Simpson. Donc le schéma sera défini par :


 y0 =η
 y i+ 12 = yi + h2 f (ti , yi )



y i+ 1 = yi + h2 f (ti + h2 , y i+ 1 )
2 2
y i+1 = yi + hf (ti + h2 , y i+ 1 )



 2
h h h

 y
i+1 = yi + 6 (f (ti , yi ) + 2f (ti + 2 , y i+ 1 ) + 2f (ti + 2 , y i+ 1 ) + f (xi+1 , y i+1 )
2 2

4.2 Généralisation
Définition 1 : On appelle schéma à un pas toute méthode d’approximation du type :

yi+1 = yi + hΦ(ti , yi , h), ∀i ∈ {0, 1, . . . , N − 1}
(M )
y0 =η
où Φ est une fonction de trois variables.
Définition 2 : La méthode (M ) est dite consistante si :
y(ti+1 − y(ti )
| − Φ(ti , y(ti ), h)|
h
tend vers zéro lorsque h tend vers zéro pour tout i.
Définition 3 : La méthode (M ) est dite d’ordre p s’il existe une constante C indépendante
de h telle que
y(ti+1 − y(ti )
| − Φ(ti , y(ti ), h)| ≤ Chp
h
pour tout i.
Définition 4 : la méthode (M ) est dite stable s’il existe deux constantes indépendantes de
h, C1 et C2 telles que :
max |yi − zi | ≤ C1 |y0 − z0 | + C2 max ki
i i

Définition 5 : la méthode (M ) est convergente si

max |yi − y(ti )| → 0
i

lorsque i → 0
Théorème 1 : Une méthode à un pas stable et consistant est convergente. Si de plus, elle
est d’ordre p alors
|ei | = |yi − y(ti )| ≤ Chp
où C est une constante indépendante de h.
Théorème 2 : Une condition suffisante de stabilité est que Φ soit lipschitzienne par rapport
à sa deuxième variable avec une constante de lipschitz indépendante de h :
|Φ(t, y, h) − Φ(t, y ∗ , h)| ≤ C|y − y ∗ |.

Draft,I.MBAYE 19 2018-2019
Université de Thiès Dpt Maths

corollaire 1 : une méthode d’ordre 1 est consitante.

4.2.1 Problème de cauchy linéaire

Soit
y 0 (t) = −λy(t), ∀t ∈ [t0 , t0 + T ]

(ED) :
y(t0 ) =η
avec λ un réel strictement positif. Stabilité de la méthode d’Euler explicite
soit le schéma d’Euler explicite suivant :

yi+1 = yi − hλyi , ∀i ∈ {0, 1, . . . , N − 1}
y0 =η
Lemme 1 : Le schéma est conditionnellement stable si et seulement si h < λ2 .
En effet, yi+1 = yi − hλyi implique yi = η(1 − hλ)i . Pour que cette suite géométrique converge
il faut et il suffit que |1 − hλ| < 1, ce qui entraine la condition de stabilité suivante : h < λ2 .
Stabilité de la méthode d’Euler implicite soit le schéma d’Euler explicite suivant :

yi+1 = yi − hλyi+1 , ∀i ∈ {0, 1, . . . , N − 1}
y0 =η
Lemme 2 : Le schéma est inconditionnellement stable.
1 1
En effet, yi+1 = yi −hλyi+1 implique yi+1 = 1+hλ yi ,d’où yi = η( 1+hλ )i . Cette suite géométrique
1
converge quelque soit h car sa raison 1+hλ < 1. Donc, le schéma d’Euler explicite est incon-
ditionnellement stable.

Draft,I.MBAYE 20 2018-2019
Université de Thiès Dpt Maths

Résolution numérique dune EDO sous matlab

 0
y = −3y, ∀t ∈ [0, 1]
y0 = 10

On définit d’abord la fonction :

function [ dy ] = exemple(t,y )
dy=-3.*y ;
end
ensuite le M-file suivant :
tmin = 0;
tmax = 1;
y0 = 10;
[t y] = ode45(0 exemple0 , [tmin tmax], y0)
plot(t, y)
grid
xlabel(0 t0 )
ylabel(0 y 0 )
10

5
y

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
t

Draft,I.MBAYE 21 2018-2019
Université de Thiès Dpt Maths

Draft,I.MBAYE 22 2018-2019
Chapitre 5

Systèmes d’équation linéaire

Les systèmes d’équation linéaire sont rencontrés dans de nombreux problème tels que la
résolution des équations aux dérivées partielles, des équations différentielles ordinaires.

5.1 Introduction
Un système d’équation linéaire s’écrit sous la forme :


 a11 x1 + a12 x2 + . . . + a1n xn = b1
a21 x1 + a22 x2 + . . . + a2n xn = b2




... .


 ... .
... .




an1 x1 + an2 x2 + . . . + ann xn = bn


Matriciellement le système s’écrit :

AX = b

où :
 
a11 a12 . . . a1n

 a21 a22 . . . a2n 

 . . ... . 
A= 

 . . ... .  
 . . ... . 
an1 an2 . . . ann

23
Université de Thiès Dpt Maths

 
b1

 b2 

 . 
b= 

 . 

 . 
bn
et  
x1

 x2 

 . 
X= 

 . 

 . 
xn
où A et b sont les matrices des données et X le vecteur des inconnues. ce système admet une
solution unique si et seulement si det(A) 6= 0.
Remarque : la résolution des systèmes d’équations linéaires avec des matrices triangulaires
ou diagonales sont trés faciles.
Compléments sur les matrices
* Une matrice A est symétrique si et seulement si At = A ;
* A est une matrice tridiagonale si et seulement si |i − j| > 1, aij = 0 ;
* Une matrice A est dite triangulaire supérieure si et seulement si pour tout i > j,
aij = 0 ;
* Une matrice A est dite triangulaire inférieure si et seulement si pour tout i < j,
aij = 0 ;
* Une matrice A est dite symétrique définie positive si et seulement si pour tout X 6= 0
on a X t AX > 0 ;
Une matrice A est dite à diagonale dominante stricte si et seulement si |aii | >
* P
i6=j |aij|  
x1
 x2 
 
 . 
 : kXk∞ = maxi |xi | ; kXk1 = n |xi | et kXk2 =
P
* Normes d’un vecteur X =   .  i=1
 
 . 
xn
pPn
2
i=1 xi ;
* Normes d’une matrice carrée A = (aij )i,j : kAk = supX6=0 kAXk
p
kXk
;kAk 2 = ρ(A∗ A) ;
Pn Pn
kAk1 = maxj i=1 |aij | et kAk∞ = maxi j=1 |aij | ;
* kABk ≤ kAkkBk  
1
Exercice : Soit le vecteur X suivant :X =  −2  calculer les normes du vecteur avec
4

Draft,I.MBAYE 24 2018-2019
Université de Thiès Dpt Maths

Matlab.
En effet on utilise la fonction norm(). on a alors :
>> X=[1 ;-2 ;4]
>> N1=norm(X,1)
>> N2=norm(X,2)
>> Ninf=norm(X,inf )
et on a les résultats suivants : N 1 = 7, N 2 = 
4.5826 et
 N inf = 4.
1 2
Exercice : Soit la matrice A suivante :A = calculer les normes de la matrice A
0 4
avec Matlab.
En effet on utilise la fonction norm(). on a alors :
>> A=[1 2 ;0 4]
>> N1=norm(A,1)
>> N2=norm(A,2)
>> Ninf=norm(A,inf )
et on a les résultats suivants : N 1 = 6, N 2 = 4.4954 et N inf = 4.

On utilisera les méthodes directes et les méthodes itératives pour la résolution des systèmes
d’équation linéaire.
Méthodes directes On appelle méthode directe de résolution d’un système d’équation
linéaire un algorithme qui permet d’obtenir la solution de ce système Une méthode en un
nombre fini d’opérations telles que (+,-,/, extraction de racines). Parmi les méthodes directes
on peut citer : la méthode de Cramer, la méthode de Pivot de Gauss, la méthode LU , la
méthode de Cholesky.
Méthodes itératives On les oppose aux méthodes directes. La méthode est itérative si on
peut construire une suite Xk qui converge vers X solution du système linéaire lorsque n tend
vers l’infini. Parmi les méthodes itérations on peut citer : la méthode de Jacobi, la méthode
de Gauss-Seidel, la méthode du gradient et la méthode du gradient conjugué.

5.2 Méthodes directes

5.2.1 Méthode LU
Soit à résoudre le système AX = b.
Elle consiste à décomposer la matrice A sous la forme A = LU où L est une matrice trian-
gulaire inférieure et U une matrice triangulaire supérieure. Cette décomposition s’appuie sur
la méthode de pivot de Gauss. On a l’équivalence suivante :

AX = b ⇔ LY = b et U X = Y

.
Proposition : Soit une matrice carrée A d’ordre n. La factorisation LU de A avec lii = 1
pour tout i ∈ {1, 2 . . . , n} existe et est unique si et seulement si les sous matrices principales

Draft,I.MBAYE 25 2018-2019
Université de Thiès Dpt Maths

Ai d’ordre i ∈ {1, 2 . . . , n − 1} sont inversibles.

Preuve  
2 −1 0
Exercice : Soit A =  −1 2 −1  . Résoudre le système d’équation suivant : AX = b
  0 −1 2
1
avec b =  3  .
4  
1 1 1
Exercice : Soit A =  2 3 4  . Déterminer la décomposition LU de A. En déduire la
2 5 7  
1
solution du système d’équation suivant : AX = b avec b =  2 .
1
 
−2 1 −1 1
 2 0 4 −3 
Exercice : Soit A =   −4 −1 −12 9  . Déterminer la décomposition LU de A. En


−2 1 1 −4
 
1
 2 
déduire la solution du système d’équations AX = b avec b =   1 .


−1
Résoudre d’un système linéaire Ax = b sous matlab
on utilise la fonction la inv() où la fonction linsolve()
>> A=[2 -1 0 ; -1 2 -1 ; 0 -1 2];
>> b=[1 ; 3 ; 4];
>> X=inv(A)*b
On peut remplacer la ligne de commande X=inv(A)*b par X = A \ b ou par X =
3.2500
linsolve(A, b) et on obtient le même résultat X = 5.5000
4.7500
Décomposition A = LU sous matlab on utilise la fonction lu() en effet :
>> A=[2 -1 0 ; -1 2 -1 ; 0 -1 2];
>> [L U]=lu(A);
on obtient
 le résultat suivant :
1 0 0
L=  −0.5000 1 0  et
 0 −0.6667 1 
2.000 −1.0000 0
U = 0 1.500 −1.0000  .
0 0 1.3333

Draft,I.MBAYE 26 2018-2019
Université de Thiès Dpt Maths

5.2.2 Méthode de Cholesky

Soit à résoudre le système AX = b. On suppose que A est symétrique définie positive.
Proposition : Une condition nécessaire et suffisante pour qu’une matrice A soit symétrique
définie positive est qu’il existe une matrice B triangulaire inférieure telle que BB t = A.
si les coefficients diagonaux de B sont choisis strictement positifs la décomposition est alors
unique
Preuve
On peut dans ce cas appliquer la méthode LU .
Algorithme de construction de B
√ a1j
B11 = a11 et Bj1 =
B11
pour 2 ≤ j ≤ n ; v
u
u i−1
X
2
Bii = taii − Bik
k=1
et Pi−1
aij − k=1 Bik Bjk
Bji =
Bii
pour tout i + 1 ≤ j ≤
n. 
9 6 3
Exercice : Soit A =  6 20 6  . Montrer que A est symétrique définie positive ensuite
3 6 3
t
appliquer l’algorithme précédent pour trouver la matrice B tele que A = BB
 . 
  4 0 12 −6
1 −2 0  0 1 2 1 
Même questions pour les matrices suivantes : A =  −2 8 −6  ; A =   12 2 49 −4
.

0 −6 25
−6 1 −4 51
Sous matlab pour retrouver la matrice triangulaire inférieure B de la matrice symétrique
définie positive A. On utilise la fonction chol(,’lower’). On les commandes Matlab sui-
vantes :
>> A=[9 6 3 ; 6 20 6 ;3 6 3]
>> B=chol(A,’lower’)
>> Bt=chol(A,’upper’)
   
3 0 0 3 2 1
On obtient les résultats suivants : B =  2 4 0  et Bt =  0 4 1 
1 1 1 0 0 1

5.3 Méthodes itératives

On suppose que la matrice A peut s’écrire sous la forme
A=M −N

Draft,I.MBAYE 27 2018-2019
Université de Thiès Dpt Maths

où M est une matrice facilement inversible.

On peut écrire le système AX = b sous la forme d’un point fixe. En effet :

AX = b ⇔ (M − N )X = b ⇔ X = M −1 N X + M −1 b

.
On
 peut définir la méthode itérative suivante :
X0 donnée
Xk+1 = M −1 N Xk + M −1 b, ∀k ∈ {0, 1, . . . , N }
La convergenge de cette méthode dépend de la matrice itérative M −1 N . En effet, on a

ek+1 = Xk+1 − X

en rempaçant
Xk+1 = M −1 N Xk + M −1 b
et
X = M −1 N X + M −1 b
dans l’expression de ek+1 on obtient alors

ek+1 = M −1 N (Xk − X) = M −1 N ek

ce qui entraine
ek = (M −1 N )k e0
. On a
lim ek = 0 ⇔ lim (M −1 N )k = 0
k→∞ k→∞

Théorème 1 : On dit que la suite des puissances Ak −→ 0 lors que k −→ 0 si et seulement

si ρ(A) < 1 ( où ρ(A) = maxλ∈σ(A) |λ|).
Preuve

Calacul des valeurs propres de A on utilise la fonction eig() en effet :

>> A=[2 -1 0 ; -1 2 -1 ; 0 -1 2];
>> [V]=eig(A);
on obtient les valeurs propres de A suivantes :
0.5858
V = 2.0000 .
3.4142
Calacul du rayon spectral de la matrice A
on utilise la fonction vrho() en effet :
>> A=[2 -1 0 ; -1 2 -1 ; 0 -1 2];
>> rho=vrho(A);
on obtient le rayon spectral de A qui est égal à : rho = 3.4142

Draft,I.MBAYE 28 2018-2019
Université de Thiès Dpt Maths

5.3.1 Méthode de Jacobi

On pose A = D − E − F avec D la matrice diagonale de A, −E la partie triangulaire
strictement inférieure de A et −F la partie triangulaire strictement supérieure de A. On
suppose que aii 6= 0 pour tout 1 ≤ i ≤ n. On choisit M = D. On obtient alors

Xk+1 = D−1 (E + F )Xk + D−1 b

.On définit l’agorithme suivant :

X0 donnée
Xk+1 = D−1 (E + F )Xk + D−1 b, ∀k ∈ {0, 1, . . . , N }
appelé la méthode de Jacobi. On pose BJ = D−1 (E + F ) la matrice itérative de Jacobi.

5.3.2 Méthode de Gauss-Seidel

On se fait les mêmes hypothèses que dans la méthode de Jacobi. On choisit M = D − E
on obtient alors
Xk+1 = (D − E)−1 F Xk + (D − E)−1 b
.On définit l’algorithme suivant :
X0 donnée
Xk+1 = (D − E)−1 F Xk + (D − E)−1 b, ∀k ∈ {0, 1, . . . , N }
appelé la méthode de Gauss-Seidel. On pose BGS = (D − E)−1 F la matrice itérative de
Gauss-Seidel.  
3 0 4
Exercice : Pour la matrice A suivante : A =  7 4 2 
−1 1 2
Déterminer la matrice itérative de Jacobi BJ et la matrice itérative de Gauss-Seidel BGS .
Conclure.
En effet : On utilise les fonctions diag() pour déterminer la matrice D, tril(A,-1) pour
déterminer la partie triangulaire inférieure stricte, triu(A,1) pour déterminer la partie tri-
angulaire supérieure stricte.
>> A=[3 0 4 ; 7 4 2 ; -1 1 2];
>> D=diag(diag(A))
>> E=-tril(A,-1)
>> F=-triu(A,1)
>> BJ=inv(D)*(E+F)
>> rhoBJ=vrho(BJ)
On obtient
 les résultats
 suivants
 :   
3 0 0 0 0 0 0 0 −4
D =  0 4 0 , E =  −7 0 0 , F =  0 0 −2 ,
0 0 2 1 −1 0 0 0 0
0 0 −1.3333
BJ =  −1.7500 0 −0.5000  et rhoBJ = 1.1251 > 1. La méthode de Jacobi ne
0.5000 −0.5000 0

Draft,I.MBAYE 29 2018-2019
Université de Thiès Dpt Maths

converge pas.
>> A=[3 0 4 ; 7 4 2 ; -1 1 2];
>> D=diag(diag(A))
>> E=-tril(A,-1)
>> F=-triu(A,1)
>> BGS=inv(D-E)*F
>> rhoBGS=vrho(BGS)
On obtient
 les résultats
 suivants
 :   
3 0 0 0 0 0 0 0 −4
D =  0 4 0 , E =  −7 0 0 , F =  0 0 −2 ,
0 0 2 1 −1 0 0 0 0
0 0 −1.3333
BGS =  0 0 1.833  et rhoBGS = 1.5833 > 1. La méthode de Gauss-Seidel ne
0 0 −1.5833
converge pas.  
−3 3 −6
Exercice : Pour la matrice A suivante : A =  −4 7 −8 
5 7 −9
Déterminer le rayon spectral de BJ et le rayon spectral de de BGS . Conclure.
en effet on a les résultats suivants :
rhoBJ = 0.8133 < 1 la méthode de Jacobi converge et rhoBGS = 1.1111 > 1 la méthode
de Gauss-Seidel ne converge pas.
Remarque : la méthode de Jacobi peut converge pour un système alors que la méthode de
Gauss-Seidel ne converge pas et inversement.  
4 1 1
Exercice : Pour la matrice A suivante : A =  2 −9 0 .
0 −8 −6
Déterminer le rayon spectral de BJ et le rayon spectral de de BGS . Conclure.
en effet on a les résultats suivants :
rhoBJ = 0.4438 < 1 la méthode de Jacobi converge et rhoBGS = 0.0185 < 1 la méthode
de Gauss-Seidel converge .
Remarque : la méthode de Gauss-Seidel converge plus rapidement que la méthode de Jacobi
car rhoBGS < rhoBJ < 1.  
7 6 9
Exercice : Pour la matrice A suivante : A =  4 5 −4 .
−7 −3 8
Déterminer le rayon spectral de BJ et le rayon spectral de de BGS . Conclure.
en effet on a les résultats suivants :
rhoBJ = 0.6411 < 1 la méthode de Jacobi converge et rhoBGS = 0.7746 < 1 la méthode
de Gauss-Seidel converge .

Draft,I.MBAYE 30 2018-2019
Université de Thiès Dpt Maths

Code de la méthode de Jacobi

clc
A=[ 2 -1 0 ;-1 2 -1 ;0 -1 2]
b=[1 ;3 ;4]
D=diag(diag(A))
E=-tril(A,-1)
F=-triu(A,1)
BJ=inv(D)*(E+F)
B=inv(D)*b
X0=[1 ; 0 ;0]
X1=BJ*X0 + B ;
eps=input(’entrer eps :’)
k=0 ;
while norm(X1-X0,inf ) > eps
X0=X1 ;
X1=BJ*X0 + B ;
k=k+1 ;
end
disp(X1)
k

5.3.3 Théorèmes de convergence de la méthode de Jacobi et Celle

de Gauss-Seidel
Théorème 1 : Si A est une matrice à diagonale dominante stricte, les méthodes de Jacobi
et de Gauss Seidel sont convergentes.
Preuve
Théorème 2 : Si A est une matrice symétrique définie positive , la méthode de Gauss Seidel
converge de manière monotone pour la norme k.kA .
Preuve
Exercice : Vérifier la convergence de la méthode de Jacobi et de Gauss-Seidel pour les
matrices
 :       
3 0 4 −3 3 −6 4 1 1 7 6 9
A1 =  7 4 2  ; A2 =  −4 7 −8  ; A3 =  2 −9 0  ; A4 =  4 5 −4  ;
−1 1 2 5 7 −9 0 −8 −6 −7 −3 8

5.3.4 Méthode SOR

On se fait les mêmes hypothèses que dans la méthode de Jacobi. On choisit
1
M = (D − ωE)
ω
et
1
N= ((1 − ω)D + ωF )
ω
Draft,I.MBAYE 31 2018-2019
Université de Thiès Dpt Maths

avec ω ∈]0, 1[. On obtient alors

1 1
Xk+1 = ( D − E)−1 N Xk + ( D − E)−1 b.
ω ω
On
 définit l’algorithme suivant :
X0 donnée
Xk+1 = ( ω1 D − E)−1 N Xk + ( ω1 D − E)−1 b, ∀k ∈ {0, 1, . . . , N }
appelé la méthode SOR (la méthode sur relaxation
 successive.)
7 6 9
Exercice : Pour la matrice A suivante : A =  4 5 −4 .
−7 −3 8
1
Déterminer le rayon spectral de BSOR et le rayon spectral de de BBSOR avec ω = 2

5.3.5 Méthode du gradient

5.3.6 Méthode du gradient conjugué

Draft,I.MBAYE 32 2018-2019
Chapitre 6

Résolution numérique des problèmes

aux limites : la méthode des
différences finies

6.1 Exemples et motivation

Nous citons quelques modèles mathématiques traduisant des phénomènes physiques et
mécaniques.
Déformation d’un fil élastique

Soit le modèle M1 suivant trouver u : Ω =]a, b[→ R :

−u00 (x) + c(x)u(x) = f (x), ∀x ∈]a, b[ (6.1)
u(a) = ga , (6.2)
u(b) = gb (6.3)
(6.4)
Où c(x) ∈ L∞ (]a, b[) un coéfficient qui dépend des propriétés du matériau (module de Young,
de l’inertie et du poids) et f ∈ L2 (]a, b[) représente la charge extérieure. Ce problème modélise
la déformation u d’un fil élastique de longueur b − a avec (a < b) fixé (ga = gb = 0) à ses
extrémités x = a et x = b. 6.1 est l’équation de conservation et 6.2 représente les conditions
aux limites.
Diffusion de la chaleur
Soit le modèle (M2 ) suivant trouver T :]a, b[×R∗+ → R tel que :

∂T (x, t) ∂ 2 T (x, t)
−α = f (x, t), ∀x ∈]a, b[ , ∀t > 0 (6.5)
∂t ∂x2
T (x, 0) = T0 (x), ∀x ∈]a, b[, (6.6)
T (a, t) = Ta (t), ∀t > 0, (6.7)
T (b, t) = Tb (t), ∀t > 0 (6.8)

33
Université de Thiès Dpt Maths

Où T est la température de la barre métallique de longueur l([a, b]) = b−a, α > 0 le coéfficient
de diffusion thermique de la barre. Ce problème modélise la diffusion de la chaleur le long
d’une barre métallique [a, b].
Déformation d’une poutre encastrée
Soit le modèle (M3 ) suivant trouver u :]a, b[→ R :

u(4) (x) + b(x)u(x) = f (x), ∀x ∈]0, l[ (6.9)

u(a) = u(b) = 0, (6.10)
u0 (a) = u0 (b) = 0 (6.11)
Où b(x) = EpEIps(x) , Es module de réaction du sol, Ep le module d’élasticité du pieu, Ip (x)
l’inertie du pieu, f ∈ L2 (]a, b[) et 0 < b ∈ L∞ (]a, b[). Ce problème modélise la déformation u
d’une poutre encastrée de longueur b − a.
Tous ces problèmes sont des appelés ”problèmes aux limites” et leur résolution analylique est
trés complexe et généralement impossible. ıous allons dans ce chapitre proposer une méthode
d’approximation à savoir la méthode des différences finies pour trouver la solution approchée
des ces types de problèmes.

6.2 Principe de la méthode des différences finies

6.2.1 Le quotient différentiel
La méthode est basée sur l’approximation des dérivées par leur quotient différentiel. En
effet soit u une fonction dérivable au point x on a alors
u(x + h) − u(x)
u0 (x) = lim
h→0 h
. Pour h ≪ 0 on pose u0 (x) approximativement égale à son quotient différentiel u(x+h)−u(x)
h
i.e
u(x + h) − u(x)
u0 (x) '
h
pour h ”trés petit” (h ≪ 0). On utilise le développement limité de u autour de x à l’ordre
2. on a alors
h2
u(x + h) = u(x) + hu0 (x) + u(2) (η)
2!
avec x < η < x + h, on obtient alors
u(x + h) − u(x) h
− u0 (x) = u(2) (η) → 0
h 2!
lorsque h → 0. Ce qui entraine l’approximation suivante
u(x + h) − u(x)
u0 (x) '
h
Draft,I.MBAYE 34 2018-2019
Université de Thiès Dpt Maths

pour h ”trés petit”.

Approximation de la dérivée première Avec le même principe on obtient les approxi-
mations suivantes de la dérivée première :
le quotient différentiel progressif

u(x + h) − u(x)
u0 (x) '
h
, le quotient différentiel rétrograde

u(x) − u(x − h)
u0 (x) '
h
et le quotient différentiel centré

u(x + h) − u(x − h)
u0 (x) '
2h
pour h ”trés petit” (h ≪ 0).
Approximation de la dérivée seconde
Avec la même idée développée précédemment on a l’approximation de la dérivée seconde
par :
u(x + h) − 2u(x) + u(x − h)
u00 (x) '
h2
pour h ”trés petit” (h ≪ 0).
Approximation de la dérivée quatrième
Avec la même idée développée précédemment on a l’approximation de la dérivée qutrième
par :
u(x − 2h) − 4u(x − h) + 6u(x) − 4u(x + h) + u(x + 2h)
u(4) (x) '
h4
pour h ”trés petit” (h ≪ 0).

6.2.2 Discrétisation ou maillage du domaine

On subdivise l’intervalle [a, b] en N + 1 intervalles de longueur h = Nb−a
+1
, appelé pas de la
discrétisation où du maillage. On pose xi = a + ih pour tout i ∈ {0, 1, 2 . . . N + 1} les points
du maillage avec x0 = a et xN +1 = b.

6.3 Application de la méthode sur le modèle (M1)

Pour chaque point intérieur xi avec i ∈ {1, 2, . . . N }. On applique l’approximation sui-
vante :
u(xi + h) − 2u(xi ) + u(xi − h)
u00 (xi ) ' ;
h2
Draft,I.MBAYE 35 2018-2019
Université de Thiès Dpt Maths

On pose ui la valeur approchée de u au point xi . On a alors

ui+1 − 2ui + ui−1
u00 (xi ) ' .
h2
L’équation (6.1) devient
ui+1 − 2ui + ui−1
− + ui = fi
h2
avec fi la valeur approchée de f au point xi . Les conditions aux bords sont discrétisées par :
u(a) = u0 = ga et u(b) = uN +1 = gb . On a le schéma à trois points suivant :

ui+1 − 2ui + ui−1

− + ci ui = fi , ∀i ∈ {1, 2, . . . N } (6.12)
h2
u0 = ga (6.13)
uN +1 = gb , (6.14)

approchant le modèle (M1 ). On dit usellement qu’on a discrétisé le problème par une méthode
de différences finies utilisant le ”schéma à trois points” de la dérivée seconde.
Matriciellement, Le schéma s’écrit : Ah Uh = Fh où les coéfficients de Ah sont définis par :

 aii = h22 + ci , ∀i ∈ {1, 2, . . . N }
aij = −1 h2
, ∀(i, j) ∈ ({1, 2, . . . N })2 et |i − j| = 1 ,
aij = 0, ailleurs


 
u1

 u2 

 . 
Uh =  

 . 

 . 
uN
et
f1 + hga2
 

 f2 

 . 
Fh =  .

 . 

 . 
fN + hgb2
Pour déterminer la solution discrète Uh , il suffit donc de résoudre le système d’équation
linéaire Ah Uh = Fh .
Proposition 1 : Supposons c ≥ 0. La matrice Ah est symétrique définie positive.
Preuve
Proposition 2 : Supposons la solution du problème continu 6.1-6.2 est de classe C4 ([a, b]).
Alors le schéma 6.12 - 6.14 est consistant d’ordre deux pour la norme k.k∞ .

Draft,I.MBAYE 36 2018-2019
Université de Thiès Dpt Maths

Preuve
Théorème 1 : On suppose c ≥ 0. Si la solution du problème continu 6.1-6.2 est de classe
C4 ([a, b]). Alors e schéma 6.12 - 6.14 d’approximation par différences finies est convergent
2
d’ordre 2 pour la norme k.k∞ . Plus précisément, on a : kUh − π(u)k ≤ h96 supx∈[a,b] |u(4) (x)|.
Preuve

6.4 Application de la méthode sur le modèle (M2)

On suppose t ∈ [0, T ]. On subdivise l’intervalle [0, T ] en M + 1 intervalles de temps de
longueur ∆t = MT+1 . On pose tn = n∆t pour tout n ∈ {0, 1, 2 . . . M + 1}.
On subdivise l’intervalle [a, b] en N + 1 intervalles de longueur h = Nb−a +1
, appelé pas de la
discrétisation où du maillage. On pose xi = a + ih pour tout i ∈ {0, 1, 2 . . . N + 1} les points
du maillage avec x0 = a et xN +1 = b.
On note Tin la valeur approchée de T au point xi à l’instant tn . Ainsi l’équation 6.5 devient :
Tin+1 − Tin un − 2uni + uni−1
− α i+1 + ui = fin
∆t h2
pour (i, n) ∈ {1, 2, . . . N } × ({0, 1, 2, . . . M + 1}).
la condition initiale 6.6 devient :Ti0 = T0 (xi ), les conditions aux bords 6.7-6.8 deviennent
respectivement T0n = Tan , TNn +1 = Tbn pour tout n ∈ {0, 1, 2, . . . M + 1}. On obtient finalement
le schéma explicite suivant :

Tin+1 −Tin T n −2Tin +Ti−1

n


 ∆t
− α i+1 h2 = fin , ∀(i, n) ∈ {1, 2, . . . N } × ({0, 1, 2, . . . M })
Ti0 = T0 (xi ), ∀i ∈ {0, 1, 2, . . . N + 1}


 T0n = Tan , ∀n ∈ {0, 1, 2, . . . M + 1}
TNn +1 = Tbn , ∀n ∈ {0, 1, 2, . . . M + 1}


et le schéma implicite suivant :

n+1
Tin+1 −Tin Ti+1 −2Tin+1 +Ti−1
n+1


 ∆t
−α h2
= fin+1 , ∀(i, n) ∈ {1, 2, . . . N } × ({0, 1, 2, . . . M })
Ti0

= T0 (xi ), ∀i ∈ {0, 1, 2, . . . N + 1}

 T0n = Tan , ∀n ∈ {0, 1, 2, . . . M + 1}
TNn +1 = Tbn , ∀n ∈ {0, 1, 2, . . . M + 1}


Ces schémas sont linéaires à deux niveaux.

Matriciellement, le schéma explicite s’écrit :
T n+1 = AT n + ∆tF n
avec A une matrice d’ordre N dont les coéfficients sont définis par :

 aii = 1 − 2c, ∀i ∈ {1, 2, . . . N }
aij = c, ∀(i, j) ∈ ({1, 2, . . . N })2 et |i − j| = 1 ,
aij = 0, ailleurs


Draft,I.MBAYE 37 2018-2019
Université de Thiès Dpt Maths

α∆t
avec c = h2
,
T1n+1
 

 T2n+1 

 . 
T n+1 = ,

 . 

 . 
n+1
TN

la condition initiale
 
T0 (x1 )

 T0 (x2 ) 

0
 . 
T = ,

 . 

 . 
T0 (x2 )

et
f1n + hα2 Tan
 

 f2n 

n
 . 
F = ∆t  .

 . 

 . 
fNn + hα2 Tbn

Avec la méthode implicite la matrice d’itération est l’inverse de A.

La stabilité du schéma est équivalente à kAn k ≤ K, pour tout n ≥ 0. ce qui veut dire que la
suite des puissances de A est bornée.
Lemme 1 : Le schéma explicite est stable en norme L∞ si et seulement si la condition
CF L : α∆t
h2
≤ 12 . Le schéma implicite est inconditionnellement stable en norme L∞ .
On peut également définir le θ-schéma ainsi :

Tin+1 −Tin T n −2T n +T n T n+1 −2T n+1 +Ti−1
n+1



 ∆t
− α(1 − θ) i+1 h2i i−1 − αθ i+1 i
h2
= (1 − θ)fin + θfin+1 ,

 ∀(i, n) ∈ {1, 2, . . . N } × ({0, 1, 2, . . . M })
Ti0 = T0 (xi ), ∀i ∈ {0, 1, 2, . . . N + 1} .
T0n = Tan , ∀n ∈ {0, 1, 2, . . . M + 1}




TNn +1 = Tbn , ∀n ∈ {0, 1, 2, . . . M + 1}



Remarque : Pour θ = 12 , on obtient le schéma de Crank-Nicolson.

Lemme 2 : Le schéma explicite est consistant précis à l’ordre 1 en temps et 2 en espace. De
plus, on choisit de garder le rapport α∆t
h2
= 16 constant, alors ce schéma est précus à l’ordre 2
en temps et 4 en espace.

Draft,I.MBAYE 38 2018-2019
Université de Thiès Dpt Maths

6.4.1 Etude de la stabilité en norme L2

On injecte dans le schéma un mode de Fourier :

Tin = An (k) exp(2jπkxi )

avec xi = ih et on en déduit la valeur du facteur d’amplificateur A(k).

On appelle condition de stabilité de Von-Neumann l’inégalité A(k) ≤ 1, ∀k ∈ Z.
Exercice : Montrer que le θ-schéma est stable en norme L2 inconditionnellement si 12 ≤ θ ≤ 1
et sous la condition CF L : 2(1 − 2θ)α∆t ≤ h2 si 0 ≤ θ ≤ 12 .
Résolution numérique d’un problème pratique soit le problème (P ) suivant

1
−u00 (x) + u(x) = P, ∀x ∈]0, 1[ (6.15)
EI
u(a) = 0, (6.16)
u(b) = 0 (6.17)
(6.18)

quitraduit la déformation d’une poutre en flexion fixée à ses extrémités, où le module de
young E = 40 l’inertie I = 10 et le poids P = 20.
on va résoudre le système Ah Uh = Fh . Dans un M-file on écrit ses lignes de code :

a = input(0 entrer la borne inf erieure a :0 );

b = input(0 entrer la borne superieure b :0 );
N = input(0 entrer le nombre dintervalles :0 );
E = input(0 entrer le module de young :0 );
I = input(0 entrer l inertie :0 );
P = input(0 entrer le poids :0 ); h = (b − a)/(N + 1);
F h = P ∗ ones(N, 1);
Ah = gallery(0 tridiag 0 , N, −1/(h ∧ 2), 2/(h ∧ 2) + 1/(E ∗ I), −1/(h ∧ 2));
U h = Ah \ F h
U = [0; U h; 0];
xi = a : h : b;
plot(xi, U )
grid
xlabel(0 Longueur de la poutre0 )
ylabel(0 Def ormation u de la poutre0 )
On obtient la solution du problème (P ) dans la figure suivante :

Draft,I.MBAYE 39 2018-2019
Université de Thiès Dpt Maths

2.5

2
Deformation u de la poutre

1.5

0.5

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Longueur de la poutre

Résolution numérique d’un problème pratique soit le problème (P ) suivant :

∂T (x, t) ∂ 2 T (x, t)
− = 0, ∀x ∈]0, 1[ , ∀t > 0 (6.19)
∂t ∂x2
T (x, 0) = exp(−x), ∀x ∈]0, 1[, (6.20)
T (0, t) = 0, ∀t > 0, (6.21)
T (0, t) = 0, ∀t > 0 (6.22)

on va résoudre le système T n+1 = AT n . Dans un M-file on écrit ses lignes de code :

a = input(0 entrer la borne inf erieure a :0 );

b = input(0 entrer la borne superieure b :0 );
N = input(0 entrer lenombre d intervalles :0 );
dt = input(0 entrer le pas de temps :0 );
alpha = input(0 entrer alpha :0 );
T 0 = inline(0 exp(−x)0 );
h = (b − a)/(N + 1);
xi = a + h : h : b − h;
T 0 = T 0(xi0 );
c = (alpha ∗ dt)/h2
if c <= 0.5
A = gallery(0 tridiag 0 , N, c, 1 − 2 ∗ c, c);
T = [];
f orn = 1 : 10
T (:, n) = (An ) ∗ T 0;
end

Draft,I.MBAYE 40 2018-2019
Université de Thiès Dpt Maths

plot(xi, T (:, 1), 0 r0 , xi, T (:, 5), 0 g 0 , xi, T (:, 10), 0 b0 )

grid xlabel(0 Longueur de la barre0 )
ylabel(0 temperature T 0 )
legend(0 T 10 , 0 T 50 , 0 T 100 )
else
disp(0 pas de stabilite0 )
end

On obtient la solution du problème (P ) dans la figure suivante :

1
T1
T5
0.9 T10

0.8

0.7
temperature T

0.6

0.5

0.4

0.3

0.2

0.1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Longueur de la barre

Draft,I.MBAYE 41 2018-2019