Introduction au calcul scientifique

Licence 2ème année - PC/SPI

Notes de cours

Petite intro ici. Ces notes de cours s’inspirent très largement du livre de Quarte-
roni, Salieri et Gervasio [1].

I Résolution d’équations
Objectif. Trouver les valeurs de x pour lesquelles une certaine fonction f : R → R
s’annule (f (x) = 0). Autrement dit, il s’agit de trouver les racines de f .
• La résolution exacte n’est pas toujours possible.
Par exemple, si f est un polynôme quelconque de degré supérieur à 4,
La situation se complique encore si f n’est pas un polynôme !
• Dans ce type de cas : résolution numérique approchée,
• Quelle méthode utiliser ? Converge-t-elle ? Avec quelle rapidité ?
• Qualité de l’approximation obtenue ?

Méthode générale. On utilise une méthode itérative qui, à partir d’une ou

plusieurs valeurs initiales (disons x0 ), construit une suite x1 , ..., xn qui, si tout se
passe bien, converge vers une racine de f .

I.1 Quelques exemples types

Fond d’investissement.
• Un client dépose v euros dans une banque au début de chanque année,
• Il retire M à la fin de la n-ième année,
• Quel est le taux annuel moyen T de cet investissement ?
n
X 1+T
M =v (1 + T )k = v [(1 + T )n − 1].
T
k=1

1−xn+1
Pn i
Rappel : k=0 x = 1−x .

1
Dans ce cas, T est une racine de la fonction f suivante :
1+T
f (T ) = M − v [(1 + T )n − 1].
T
Équation d’état d’un gaz réel.
• On considère n moles de gaz à une température T et une pression p,
• Quel est le volume V occupé par ce gaz ?
• Équation de van der Waals (Johannes Diderik van der Waals, 1837 – 1923)
pour un gaz réel (6= gaz parfait) :
  n 2 
p+a (V − nb) = nRT ,
V
Ici, a et b sont des coefficients qui dépendent du gaz considéré, et R est la constante
des gaz parfaits (R ≈ 8.3144621 J·mol−1 ·K−1 ).
Nous obtenons une équation polynomiale de degré 3 en V . Une méthode analy-
tique est possible, mais fastidieuse (méthode de Cardan).

Mécanique Statique.
• On considère le système de quatre barres a1 , a2 , a3 , et a4 suivant :
y
a3

a4
a2
¯
®

a1 x

• Pour un angle β donné, nous voulons déterminer la (les ?) valeur(s) de α.

• Il suffit d’introduire l’angle γ entre l’horizontale et la barre a3 , de projeter
les barres a2 , a3 , et a4 sur les axes x et y, et d’éliminer γ entre les deux
équations obtenues (en utilisant l’identité cos2 (γ) + sin2 (γ) = 1). On obtient
l’équation de Freudenstein (Ferdinand Freudenstein, 1926 – 2006) :
a1 a1 a21 + a22 − a23 + a24
cos(β) − cos(α) − cos(β − α) = − ,
a2 a4 2a2 a4

2
 où ai est la longueur de la barre ai . On cherche donc les racines de la fonction
a1 a1 a2 + a22 − a23 + a24
f (x) = cos(β) − cos(x) − cos(β − x) + 1 .
a2 a4 2a2 a4
• Ce problème n’admet des solutions que pour certaines valeurs de β,
• Il y a plusieurs solutions pour certaines valeurs de β,
• Il n’y a pas d’expression explicite de ces solutions dans le cas général.

Dynamique des populations.

• Appelons xk le nombre d’individus (ex : bactéries) à la génération k.
• La valeur de xk+1 est une fonction de xk :

xk+1 = φ(xk ). (1)

• En général, φ(x) = xR(x) où R(x) modélise la vitesse d’évolution de la

population considérée, et peut être choisie de différentes manières.
Modèle de Malthus (Thomas Malthus, 1766 – 1834) :

R(x) = r ; r > 0.

Modèle de croissance avec ressources limitées (Pierre François Verhulst,

1804 – 1849). Dans ce cas, la croissance de la population est limitée par
la disponibilité des ressources :
r
R(x) = ; r > 0, K > 0.
1 + xK
Modèle prédateurs-proies avec saturation. Ce modèle prend en compte une
population antagoniste :
rx
R(x) = .
1 + (x/K)2
• Un état stationnaire est une taille de population x? pour laquelle

x? = φ(x? ),

ou, de manière équivalente R(x? ) = 1. Trouver la valeur de x? revient à

trouver les racines de f (x) = R(x) − 1. Remarquons ici que (1) est un
exemple d’itération du point fixe.

3
I.2 Dichotomie
• Considérons une fonction f continue sur un intervalle [a, b],
• Si f (a)f (b) < 0 (c’est-à-dire si f (a) et f (b) sont de signes opposés), alors f
a au moins une racine dans ]a, b[ (Théorème des valeurs intermédiaires).
• Supposons que cette racine est unique. La méthode de dichotomie consiste à
générer une suite d’intervalles ]ak , bk [ de plus en plus petits, qui contiennent
cette racine. Le centre xk de ces intervalles converge vers la racine de f .

Poser a0 = a, b0 = b, k = 0.
Tant que k < n, faire
Poser xk = (ak + bk )/2.
Si f (xk ) = 0, alors
Retourner xk .
Sinon
Si f (ak )f (xk ) < 0,
Poser ak+1 = ak , bk+1 = xk .
Si f (xk )f (bk ) < 0,
Poser ak+1 = xk , bk+1 = bk .
Incrémenter k de 1.

• La méthode s’arrète lorsque k = n ou si l’on trouve f (xk ) = 0 pour k < n.

Exemple. Avec f (x) = x2 − 1 et l’intervalle initial ] − 0.25, 1.25[, on obtient :

I0 =] − 0.25, 1.25[, x0 = 0.5,

I1 =]0.5, 1.25[, x1 = 0.875,
I2 =]0.875, 1.25[, x2 = 1.0625,
I3 =]0.875, 1.0625[, x3 = 0.96875.

• L’erreur à l’étape k est majorée par la demi-longueur de l’intervalle ]ak , bk [,

• Il est possible d’implémenter un test d’arrêt utilisant cette majoration. Il
suffit d’effectuer les itérations tant que ε < bk − ak où ε est une valeur de
l’erreur maximale souhaitée.
• On peut montrer que pour une valeur de ε donnée, il suffit d’effectuer k

4
 itérations pour obtenir ε < bk − ak où :
 
b−a
k > log2 − 1.
ε
Ceci sera montré en séance de TD.

I.3 Méthode de Newton

• Si f est dérivable, nous pouvons exploiter les valeurs de sa dérivée.
• Nous allons faire comme si f est localement (approximativement) affine.
• La fonction y tangente de f en xk s’écrit :

y(x) = f (xk ) + f 0 (xk )(x − xk ).

• Si on approche f autour de xk par y, on peut confondre la racine de y avec

une racine de xk . Cette racine (appelons-la xk+1 ) s’écrit :
f (xk )
xk+1 = xk − . (2)
f 0 (xk )
• On ne peut écrire cette égalité que si f 0 (xk ) 6= 0 !
• Puisque f n’est en général pas affine, on utilise (2) pour construire une suite
x1 , ..., xn à partir d’une valeur initiale x0 . Ce sont les valeurs de cette suite
qui nous permettent d’approcher une racine de f .
• En pratique, le choix de x0 est déterminant pour que la suite xk converge
vers une racine de f . On peut par exemple choisir x0 en examinant le graphe
de f ou en effectuant quelques itérations préliminaires de la méthode de
dichotomie.
• Voici un pseudo-code pour cette méthode. Remarquons qu’il est inutile de
définir une variable xi par itération puisque seule la dernière de ces valeurs
nous intéresse. Ce pseudo-code évite cela en n’utilisant que la variable x dont
la valeur est modifiée à chaque itération. En outre, il faut introduire une
variable supplémentaire y qui contient la valeur précédente, afin de pouvoir
écrire le test d’arrêt (on effectue des itérations tant que |xi+1 − xi | > ε
c’est-à-dire tant que |x − y| > ε).

Poser y = x0 .
Poser x = y − f (y)/f 0 (y).

5
 Tant que |y − x| > ε faire
Poser y = x.
Poser x = x − f (x)/f 0 (x).
Fin tant que
Renvoyer x.

Exemple. Reprenons la fonction de l’exemple précédent, f (x) = x2 − 1. Dans ce

cas, f 0 (x) = 2x. Avec le terme initial x0 = 2, on obtient :
x20 −1 2
x1 = x0 − 2x0 = 2 − 2 2−1
2 = 54 = 1.25,
x21 −1
x2 = x1 − 2x1 = 54 − 25/16−1 41
10/4 = 40 = 1.025,
x22 −1 41 1681/1600−1 3281
x3 = x2 − 2x2 = 40 − 82/40 = 3280 ≈ 1.000305.

• Si la suite xk converge vers une racine simple α de f (c’est-à-dire si f 0 (α) 6= 0)

et si f est deux fois dérivable, alors on a le résultat suivant :
xk+1 − α f 00 (α)
lim = 0 .
k→∞ (xk − α)2 2f (α)
On dit alors que la méthode de Newton a une convergence quadratique. Ce
résultat se prouve en utilisant la formule de Taylor Lagrange à l’ordre 2 :

0 f 00 (ξ)
f (α) = f (xk ) + f (xk )(α − xk ) + (α − xk )2 .
2
Rappelons que, dans cette formule, ξ est dans l’intervalle ouvert de bornes
α et xk . Puisque f (α) = 0, et selon (2), on obtient :
xk+1 − α f 00 (ξ)
= 0 .
(xk − α)2 2f (xk )
Il suffit alors de passer à la limite quand k → ∞.
• Si la suite xk converge vers une racine α de f de multiplicité m (c’est-à-
dire si f 0 (α) = 0, ..., f (m−1) (α) = 0 et f (m) (α) 6= 0), alors la convergence
de la méthode n’est plus quadratique, mais linéaire. On peut retrouver une
convergence quadratique en modifiant (2) de la manière suivante :
f (xk )
xk+1 = xk − m .
f 0 (xk )
• Ceci necessite de connaître m a priori...
• Comme pour la dichotomie, nous pouvons faire un test d’arrêt qui prend en

6
 compte l’erreur maximale voulue. En effet, la différence entre deux itérations
successives fournit un estimateur correct de l’erreur. Ainsi, on peut effectuer
les itérations (2) jusqu’à ce que

|xk+1 − xk | < ε.

• On peut aussi faire un test d’arrêt sur le résidu (c’est-à-dire sur la valeur de
f (xk )). On arrête alors les itérations lorsque |f (xk )| < ε, mais ce test n’est
satisfaisant que lorsque |f 0 (x)| ≈ 1 au voisinage de la racine. Si |f 0 (x)|  1
alors on sous-estime l’erreur et si |f 0 (x)|  1, alors on sur-estime l’erreur.

I.4 Méthode de point fixe

• On observe que la suite xk+1 = cos(xk ) avec x0 = 1 converge vers 0.73908513.
• Si la suite définie pour un certain x0 par :

xk+1 = φ(xk ),

converge, alors elle converge vers α tel que α = φ(α). On appelle cette mé-
thode la méthode du point fixe. La fonction φ est appelée fonction d’itération.
• On dit que α est un point fixe de φ. C’est aussi une racine de f (x) = x−φ(x).
• Toutes le fonctions n’ont pas de point fixe. Par exemple φ(x) = exp(x).
• Le méthode ne converge pas forcément même si φ admet un point fixe.

Remarque. La méthode de Newton peut être formulée comme une méthode de

point fixe avec la fonction d’itération
f (x)
φ(x) = x − .
f 0 (x)

I.5 Théorèmes de convergence

Théorème du point fixe.
i. Si φ est continue sur [a, b] et si φ(x) ∈ [a, b] pour tout x ∈ [a, b], alors cette
fonction admet au moins un point fixe.
ii. De plus, si φ est L-Lipschitzienne sur [a, b] avec L < 1, alors φ admet un
unique point fixe α dans [a, b] et la méthode du point fixe converge vers α

7
 pour tout x0 ∈ [a, b]. On dit que φ est L-Lipschitzienne sur [a, b] lorsque :

|φ(y) − φ(x)| ≤ L|y − x|,

pour tout x et y dans [a, b].

Preuve.
i. Théorème des valeurs intermédiaires avec g(x) = φ(x) − x.
ii. Contradiction avec L < 1 si il y a deux points fixes. En outre, puisque φ est
L-Lipschitzienne, on prouve :
|xk − α|
≤ Lk .
|x0 − α|

• En pratique, il est difficile de trouver un intervalle [a, b] qui satisfait les

hypothèses de la proposition. Si x0 est suffisamment proche d’un point fixe
de φ, on a le résultat local suivant :

Théorème d’Ostrowski.
• On considère un point fixe α d’une fonction φ dérivable au voisinage de α.
On suppose que φ0 est continue.
• Si |φ0 (α)| < 1, alors la méthode du point fixe converge vers α dès que x0 est
sufisamment proche de α. De plus :
xk+1 − α
lim = φ0 (α).
k→∞ xk − α

• Dans le cadre des deux résultats précédents, la convergence est au moins

linéaire : lorsque k est assez grand, l’erreur d’une étape à la suivante est
multipliée par une constante plus petite que 1 (L ou |φ0 (α)|) qui ne dépend
pas de k. Ce coefficient s’appelle le coefficient de convergence asymptotique,
et la convergence est d’autant plus rapide que de coefficient est petit.
• Dans le cadre du théorème d’Ostrowski, si |φ0 (α)| > 1, alors la méthode du
point fixe ne peut pas converger. En effet, si xk s’approche suffisamment de
α, on aura alors
xk+1 − α
> 1,
xk − α
et l’erreur ne pourra pas tendre vers 0.

8
 • Rappelons que la méthode de Newton est une méthode de point fixe avec :
f (x)
φ(x) = x − .
f 0 (x)
Et que la méthode converge de manière quadratique vers α lorsque f 0 (α) 6= 0.
En fait, plus généralement :

Proposition.
• On considère un point fixe α d’une fonction φ deux fois dérivable au voisinage
de α. On suppose que φ00 est continue et que la méthode du point fixe converge
vers α pour un certain x0 .
• Si φ0 (α) = 0, alors :
xk+1 − α 1 00
lim = φ (α).
k→∞ (xk − α)2 2

• D’après la proposition suivante, si φ0 (α) = 0, l’erreur de la méthode du point

fixe est correctement approchée par |xk+1 − xk | lorsque xk est suffisamment
proche de α. Notons que dans le cas de la méthode de Newton, on a bien
φ0 (α) = 0, lorsque α est une racine simple de f , ce qui justifie le critère
d’arrêt que nous avons énoncé dans ce cas.

Proposition (Test d’arrêt).

Si φ est dérivable et φ0 est continue alors :
α − xk 1
lim = .
k→∞ xk+1 − xk 1 − φ0 (α)

I.6 Méthode de Newton-Horner

• La méthode de Newton-Horner permet d’approcher les racines d’un poly-
nôme de degré quelconque.
• Elle combine la méthode de Newton avec l’algorithme de Horner, qui permet
de calculer efficacement la valeur d’un polynôme en un point.

Décrivons d’abord la méthode de Horner. Considérons, par exemple, le polynôme

p(x) = 1 + 2x + 3x2 + 4x3 + 5x4 .

Pour calculer p4 (x), il faut calculer les monômes

• a0 (aucune opération)

9
 • a1 ×x (une multiplication)
• a2 ×x×x (deux multiplications)
• a3 ×x×x×x (trois multiplications)
• a4 ×x×x×x×x (quatre multiplications)

Il faut donc 10 multiplications. Pour sommer les 5 monômes, il faut 4 additions.

Nous pouvons faire mieux en utilisant une écriture différente du polynôme p.
Factorisons d’abord x de la façon suivante :

p(x) = 1 + x(2 + 3x + 4x2 + 5x3 ).

Nous pouvons faire une factorisation similaire à l’intérieur de la parenthèse, c’est-

à-dire pour le polynôme 2 + 3x + 4x2 + 5x3 :

p(x) = 1 + x(2 + x(3 + 4x + 5x2 )).

Le polynôme 3 + 4x + 5x2 peut lui aussi être réécrit de la même façon :

p(x) = 1 + x(2 + x(3 + x(4 + 5x))).

En utilisant cette écriture, nous devons calculer :

• b3 = 4 + 5×x (une multiplication et une addition)
• b2 = 3 + b3 ×x (une multiplication et une addition)
• b1 = 2 + b2 ×x (une multiplication et une addition)
• b0 = 1 + b1 ×x (une multiplication et une addition)

À la fin de ce calcul, on remarque que b0 est égal à p(x). Il ne nous a fallu que 4
multiplications et 4 additions. Cette remarque se généralise pour des polynômes
quelconques. Considérons le polynôme
n
X
pn (x) = ak xk . (3)
k=0

Pour calculer la valeur de pn (x) en utilisant l’écriture (3), il est nécéssaire de faire
n sommes et n(n + 1)/2 produits (car pour calculer le monôme de degré i, il faut
i multiplications). Cependant, pn (x) peut être réécrit

pn (x) = a0 + x (a1 + x (a3 + ... + x (an−1 + an x))) . (4)

En utilisant cette écriture, pour calculer pn (z) il suffit de faire n additions et n

10
multiplications. Pour cela, il s’agit de poser bn = an et de calculer, récursivement

bk = ak + bk+1 z, (5)

lorsque k varie de n − 1 à 0. On obtient b0 = pn (z). Cet algorithme est la méthode

de Horner et il s’écrit de la façon suivante en pseudo-code :

Poser bn = an .
Pour k = n − 1, n − 2, ..., 0, faire
Poser bk = ak + bk+1 z.
Fin pour

Une remarque importante est que tous les bk calculés par la méthode dépendent
de z. Rappelons en particulier que b0 est égal à pn (z).

Cette méthode permet de factoriser pn de la façon suivante.

Proposition. Si z est une racine de pn , alors :

pn (x) = (x − z)qn−1 (x; z),

où n
X
qn−1 (x; z) = bk xk−1 .
k=1

Preuve. D’après (5), pour tout k dans {0, 1, ..., n − 1},

ak = bk − bk+1 z.

Par conséquent,
ak xk = bk xk − bk+1 zxk .

On obtient alors l’expression suivante de pn (x) :

n−1
X n−1
X
n k
pn (x) = an x + bk x − bk+1 zxk .
k=0 k=0

Rappelons que bn = an et que b0 = pn (z) = 0. Par conséquent

n−1
X
pn (x) = xqn−1 (x; z) − z bk+1 xk .
k=0

11
Il suffit maintenant de réindexer la somme restante pour s’aperçevoir que
n−1
X
bk+1 xk = qn−1 (x; z).
k=0

La proposition découle de ces deux dernières égalités. 

• D’après cette proposition, les racines de pn sont z et les racines de qn−1 (en
tant que polynôme de la variable x).
• Ainsi, nous pouvons trouver toutes les racines de pn par la méthode suivante,
appelée méthode de déflation :

Pour m = n, n − 1, ..., 1, faire

Trouver une racine rm de pm par une méthode d’approximation,
Calculer qm−1 (x; rm ) par la méthode de Horner,
Poser pm−1 (x) = qm−1 (x, rm ).
Fin pour

La méthode de Newton-Horner consiste à utiliser la méthode de Newton comme

méthode d’approximation pour calculer rm , et à utiliser la propriété suivante pour
exprimer la dérivée de pn dans l’équation d’itération :

Proposition. p0n (x) = qn−1 (x; x)

Preuve. Si z = x, (5), devient

ak = bk − bk+1 x.

Par conséquent, si k > 0,

kak xk−1 = kbk xk−1 − kbk+1 xk .

Sommons ces égalités lorsque k varie de 1 à n − 1 :

n−1
X n−1
X n−1
X
k−1 k−1
kak x = kbk x − kbk+1 xk .
k=1 k=1 k=1

Remarquons qu’en réindiçant la deuxième somme,

n−1
X n−1
X
k
kbk+1 x = (k − 1)bk xk−1 + (n − 1)bn xn−1 .
k=1 k=1

12
Par conséquent,
n−1
X n−1
X
k−1
kak x = bk xk−1 − (n − 1)bn xn−1 .
k=1 k=1

Si nous ajoutons nbn xn−1 aux deux membres de cette égalité, nous obtenons, dans
le membre de droite p0n (x) et dans le membre de gauche, qn−1 (x; x). 

• Ainsi, les itérations de la méthode de Newton s’écrivent :

pn (xk ) pn (xk )
xk+1 = xk − 0
= xk − .
pn (xk ) qn−1 (xk ; xk )
• En appliquant la méthode de Newton, nous allons donc obtenir les coefficients
de qn−1 en calculant pn (xk ) par la méthode de Horner, et nous les utiliserons
ensuite pour calculer qn−1 (xk ; xk ) de nouveau par la méthode de Horner.

Exemple. Considérons le polynôme p2 (x) = x2 − 3x + 2, de racines 1 et 2.

Efectuons une itération de la méthode de Newton en partant de x0 = 3.
• Calculons p2 (3) par la méthode de Horner. Ici a2 = 1, a1 = −3 et a0 = 2 :
b2 = a2 = 1 ; b1 = a1 + 3b2 = 0 ; b0 = a0 + 3b1 = 2.
• Dans ce cas p02 (3) = q1 (3; 3) = 3b2 + b1 = 3, donc
2 7
x1 = 3 − = .
3 3

Exercice. Considérons le polynôme

p3 (x) = x3 − 6x2 + 11x − 6

a) Appliquez la méthode de Newton-Horner en partant de x0 = 0. Vers quel

nombre α la méthode semble t-elle converger ?
b) Calculez q2 (x; α).
c) Déterminez les racines de q2 (x; 2).
d) Déterminez toutes les racines de p3 .

Exercice. Considérons le polynôme

p3 (x) = x3 − 7x − 6

a) Appliquez la méthode de Newton-Horner en partant de x0 = 0. Vers quel

13
 nombre α la méthode semble t-elle converger ?
b) Calculez q2 (x; α).
c) Déterminez les racines de q2 (x; 2).
d) Déterminez toutes les racines de p3 .

II Systèmes linéaires
Considérons le système hydraulique suivant :

x1

5
x3

x2

• Le flux qui passe dans chaque conduite est indiqué à côté de la conduite,
• Les flux positifs vont dans le sens de la flêche,
• Les flux négatifs vont dans le sens inverse à la flêche.
• Que valent les flux x1 , x2 et x3 ? On écrit la conservation du flux aux (trois)
noeuds afin de trouver le système :

 x1 + x2
 = 5,
x1 − x3 = 3,

x2 + x3 = 2.


• Est-ce que ce système admet des solutions ? Oui.

• Est-ce que ce système admet plusieurs solutions ? Oui.

14
II.1 Cas général
Considérons un système de m équations linéaires avec n variables (ou inconnues)
x1 , ..., xn :
a1,1 x1 + · · · + a1,n xn = b1 ,
.. .. .. ..
. . . .
am,1 x1 + · · · + am,n xn = bm ,

Objectif. Trouver tous les n-uplets (x1 , ..., xn ) qui satisfont ce système.

II.2 La méthode du pivot de Gauss sur des exemples

Les règles de manipulation d’un système d’équations sont les suivantes :

• Nous pouvons multiplier une équation par une constante non nulle,
• Nous pouvons ajouter un certain nombre de fois une équation à une autre,
• Nous pouvons échanger deux équations.

Objectif. Se ramener à un système avec (autant que possible) une variable par
ligne et par colonne.

Exemple. Considérons le système :


 x + 4y − z = 6,

2x + 4y + 3z = 19, (6)

3x + 5y − 2z = 7.


Pour faire disparaitre le monôme 2x dans la première colonne, il suffit de retran-

cher 2 fois la première ligne à la deuxième :

 x + 4y − z = 6,

− 4y + 5z = 7,

3x + 5y − 2z = 7.


Il est possible ensuite de faire disparaitre le monome 3x de la première colonne

15
en retranchant 3 fois la première ligne à la troisième :

 x + 4y − z = 6,

− 4y + 5z = 7,

− 7y + z = −11.


Ajoutons maintenant la troisième ligne à la première et retranchons 5 fois la

troisième ligne à la deuxième :

 x − 3y
 = −5,
31y = 62,

− 7y + z = −11.


Divisons d’abord la deuxième ligne par 31, puis Ajoutons 3 fois cette deuxième
ligne à la première et 7 fois cette deuxième ligne à la troisième :

x
 = 1,
y = 2,

z = 3.


L’unique solution du système (6) se lit ici directement : c’est le triplet (1, 2, 3).

Remarque. Ces calculs peuvent être fait sur la matrice augmentée du système,
c’est-à-dire sur le tableau de nombre obtenu en enlevant les variables, et en rem-
plaçant la colonne des signes = par un trait vertical. C’est de cette manière là que
nous procéderons. Par exemple, pour le système ci-dessus :
 
1 4 -1 6
2 4 3 19 
 

3 5 -2 7

Méthode. La méthode du pivot de Gauss consiste à effectuer les opérations

utilisées plus haut, mais dans un ordre bien précis, en particulier pour le choix de
la colonne. Voilà trois premières règles qui définissent cet ordre :

• On traite les colonnes de gauche à droite. La colonne en cours de traitement

s’appelle colonne du pivot,
• Lorsqu’on traite une colonne, on élimine tous les termes de la colonne sauf

16
 un avant de passer à la colonnee suivante.
• Lorsqu’on traite une colonne, les éliminations se font toujours avec la même
ligne. Cette ligne s’appelle la ligne du pivot. On change de ligne du pivot à
chaque fois qu’on change de colonne du pivot.

Reprenons notre exemple avec ces règles :

 
1 4 -1 6
2 4 3 19 
 

3 5 -2 7
On retranche 2 fois la première ligne à la deuxième, et 3 fois la première ligne à
la troisième. Pour le traitement de la première colonne, la ligne du pivot est la
première ligne :
 
1 4 -1 6
0 -4 5 7
 

0 -7 1 -11
Pour le traitement de la deuxième colonne, la ligne du pivot est la deuxième ligne.
On divise cette ligne par −4. Ceci fait apparaitre un 1 directeur en deuxième
colonne-deuxième ligne. Les 1 directeurs se trouvent au croisement de la colonne
et de la ligne du pivot (on remarque donc que le 1 en première ligne, première
colonne est aussi un 1 directeur) :
 
1 4 -1 6
0 1 -5/4 -7/4 
 

0 -7 1 -11
On retranche 4 fois la deuxième ligne à la première, et on ajoute 7 fois la deuxième
ligne à la troisième :
 
1 0 4 13
0 1 -5/4 -7/4
 

0 0 -31/4 -93/4
Divisons maintenant la troisième ligne par −31/4 pour faire apparaitre un 1
directeur (la dernière ligne du pivot sera la troisième) :

17
  
1 0 4 13
0 1 -5/4 -7/4 
 

0 0 1 3
Enfin, retranchons 4 fois la troisième ligne à la première et ajoutons là 5/4 de fois
à la deuxième :
 
1 0 0 1
0 1 0 2
 

0 0 1 3
Méthode (suite). Les trois règles énoncées ci-dessous ne suffisent pas pour dé-
crire la méthode du pivot de Gauss de manière complète. En effet, dans la méthode
du pivot de Gauss,

• La ligne du pivot de la première colonne est toujours la première ligne,

• La ligne du pivot des autres colonnes est toujours la ligne immédiatement
sous la précédente ligne du pivot.

Pour faire en sorte que ces règles soient vérifiées, nous allons devoir éventuellement
échanger des lignes ! Considérons un nouvel exemple :



 x2 + 2x3 + 3x4 = 9,

 x + 2x
1 2 + 3x3 + 4x4 = 15,


 x1 + x2 + x3 + x4 = 6,

x + x
1 2 + 2x3 + x4 = 8.
La matrice augmentée de ce système s’écrit :
 
0 1 2 3 9
 
1 2 3 4 15 
 
1 1 1 1 6 
 
1 1 2 1 8
Commençons par échanger la première et la deuxième ligne. En général, si il y
a un 0 là où le pivot doit être placé, on l’échange avec la première des lignes
suivantes qui n’a pas de 0 dans la colonne à traiter. Si toutes les lignes suivantes
ont un 0 dans cette colonne, on considère que cette colonne est traitée et on passe

18
à la suivante. Il n’y aura alors pas de 1 directeur dans cette colonne !
 
1 2 3 4 15
 
0 1 2 3 9 
 
1 1 1 1 6 
 
1 1 2 1 8
Retranchons la première ligne à la troisième et à la quatrième :
 
1 2 3 4 15
 

 0 1 2 3 9 


 0 -1 -2 -3 -9 

0 -1 -1 -3 -7
Retranchons deux fois la deuxième ligne à la première et ajoutons-là à la troisième
et à la quatrième :
 
1 0 -1 -2 -3
 
0 1 2 3 9 
 
0 0 0 0 0 
 
0 0 1 0 2
Il faut maintenant échanger la troisième et la quatrième ligne pour placer la ligne
du pivot en troisième ligne :
 
1 0 -1 -2 -3
 
0 1 2 3 9 
 
0 0 1 0 2 
 
0 0 0 0 0
Ajoutons enfin la troisième ligne à la première ligne et retranchons-là deux fois à
la deuxième ligne :
 
1 0 0 -2 -1
 

 0 1 0 3 5 


 0 0 1 0 2 

0 0 0 0 0
La quatrième colonne est traitée puisque nous ne pouvons plus fabriquer de ligne

19
du pivot. Ici le système admet plusieurs solutions (c’est le cas quand le nombre
de 1 directeurs est inférieur au nombre d’inconnues). Si on réécrit cette matrice
augmentée comme un système d’équations, on obtient :

 x1 = −1 + 2x4 ,

x2 = 5 − 3x4 ,

x3 = 2.


L’inconnue x4 est devenue un paramètre, c’est-à-dire que, pour chaque valeur de

x4 (dans R), il y a une solution (x1 , x2 , x3 , x4 ) du système obtenue en exprimant x1
et x2 en fonction de x4 . La variable x3 , elle, sera toujours égale à 2. Par exemple,
pour x4 = 2,

(x1 , x2 , x3 , x4 ) = (3, −1, 2, 2),

est une solution du système. Pour x4 = −5,

(x1 , x2 , x3 , x4 ) = (−11, 38, 2, −5),

est une autre solution. Pour obtenir toutes les solutions du système, il faut faire
varier les paramètres dans tout R. On écrit l’ensemble des solutions de la manière
suivante :
      


 x1 −1 2 



      
 x2   5   −3 
 ; x4 ∈ R

S=  =  + x4 




 x3   2 
  

 0 
 


 
 x4 0 1 

Les exemples précédents étaient des systèmes carrés. Appliquons la méthode à un

système rectangulaire, où le nombre de lignes est inférieur au nombre de colonnes :



 x1 + 2x2 − 2x3 + 2x4 + x5 = 1,

 x
1 + x2 + x3 + x4 + x5 = 3,


 2x1 + x2 + x3 + x5 = 3,

 3x
1 + 2x2 + 2x3 + x4 + 2x5 = 6.
C’est un système linéaire de 4 équations à 5 inconnues, dont la matrice augmentée
s’écrit :

20
  
1 2 -2 2 1 1
 

 1 1 1 1 1 3 


 2 1 1 0 1 3 

3 2 2 1 2 6
Retranchons 1 fois la première ligne de la deuxième, 2 fois de la troisième et 3
fois de la quatrième :
 
1 2 -2 2 1 1
 

 0 -1 3 -1 0 2 


 0 -3 5 -4 -1 1 

0 -4 8 -5 -1 3
Multiplions la deuxième ligne par −1.
 
1 2 -2 2 1 1
 
0 1 -3 1 0 -2 
 
0 -3 5 -4 -1 1 
 
0 -4 8 -5 -1 3
Retranchons 2 fois la deuxième ligne de la première et ajoutons-là 3 fois à la
troisième et 4 fois à la quatrième :
 
1 0 4 0 1 5
 
0 1 -3 1 0 -2 
 
0 0 -4 -1 -1 -5 
 
0 0 -4 -1 -1 -5
Divisons la troisième ligne par −4 :
 
1 0 4 0 1 5
 
0 1 -3 1 0 -2 
 
0 0 1 1/4 1/4 5/4
 
0 0 -4 -1 -1 -5
Retranchons 4 fois la troisième ligne de la première et de la quatrième et ajoutons-
là 3 fois à la deuxième :

21
  
1 0 0 -1 0 0
 

 0 1 0 7/4 3/4 7/4


 0 0 1 1/4 1/4 5/4

0 0 0 0 0 0
La quatrième et la cinquième colonne correspondent à un paramètre chacune.
C’est-à-dire qu’il y a deux paramètres pour ce systèmes (x4 et x5 ). L’ensemble
des solutions s’écrit :
        

 x1 0 1 0 

 
x 7/4 −7/4 −3/4

         


  2        

      
S =  x3  =  5/4  + x4  −1/4  + x5  −1/4  ; x4 , x5 ∈ R
       
 
x 0 1 0

         


  4        

 
x5 0 0 1
 

Considérons maintenant le cas d’un système rectangulaire dont le nombre d’équa-

tions est supérieur au nombre de variables :



 x1 + x2 = 0,

 2x1 + x2 = −1,


 3x1 + 2x2 = −1,

 −5x1 − 3x2 = 2.
C’est un système de quatre équations à deux inconnues, de matrice augmentée :
 
1 1 0
 
2 1 -1 
 
3 2 -1 
 
-5 -3 2
Retranchons 2 fois la première ligne de la deuxième et 3 fois de la troisième, et
ajoutons-là 5 fois à la quatrième :
 
1 1 0
 
0 -1 -1 
 
0 -1 -1 
 
0 2 2
Multiplions la deuxième ligne par −1 :

22
  
1 1 0
 

 0 1 1 


 0 -1 -1 

0 2 2
Retranchons la deuxième ligne de la première et 2 fois de la quatrième, et ajoutons-
là à la troisième :
 
1 0 -1
 

 0 1 1 


 0 0 0 

0 0 0
L’ensemble des solutions s’écrit :
( ! !)
x1 −1
S= =
x2 1

Un système peut ne pas avoir de solutions. Considérons par exemple le système



 x + y + 2z = 5,
−x + y = 6,

2x + y + 3z = 1.


La matrice augmentée de ce système s’écrit :

 
1 1 2 5
 -1 1 0 1 
 

2 1 3 6
Ajoutons la première ligne à la deuxième et retranchons 2 fois cette ligne à la
troisième :
 
1 1 2 5
0 2 2 6 
 

0 -1 -1 -5
Divisons la deuxième ligne par 2 :

23
  
1 1 2 5
0 1 1 3 
 

0 -1 -1 -5
Retranchons la deuxième ligne de la première et ajoutons-là à la troisième :
 
1 0 1 2
0 1 1 3 
 

0 0 0 -2
Nous ne pouvons plus continuer. La méthode s’arrête. La troisième ligne de cette
matrice augmentée se lit 0 = −2. Dans ce cas (c’est-à-dire, en général, lorsque
le membre de gauche d’une équation est nul et que son membre de droite est
non-nul), le système n’admet pas de solutions !

Recapitulatif. Lorsqu’on résoud un système linéaire (de n équations à m in-

connues aui n’est pas forcément carré puisque n et m peuvent être différents), la
méthode du pivot de Gauss permet de transformer la matrice augmentée initiale
en une matrice dite échelonnée réduite, dans laquelle :
• Les colonnes qui contiennent un 1-directeur ont des zéros partout ailleurs,
• Chaque 1-directeur se trouve dans la ligne qui est immédiatement sous le
1-directeur précédent (qui n’est pas forcément dans la colonne précédente),
• Les colonnes qui ne contiennent pas de 1-directeur ont des zéros sous le
1-directeur précédent.

Le système admet :
• Exactement une solution lorsque toutes les colonnes de la matrice échelonnée
réduite contiennent un 1-directeur et que les lignes qui ne contiennent pas
de 1-directeur ne contiennent que des zéros,
• Aucune solution lorsque le membre de gauche d’une ligne de la matrice éche-
lonnée réduite ne contient que des zéros et que son membre de droite est
non-nul,
• Plusieurs (une infinité) de solutions dans les autres cas. Les solutions sont
alors écrites en fonction de paramètres qui peuvent prendre n’importe quelle
valeur réelle. Ces paramètres sont les variables des colonnes dans lesquelles
la matrice échelonnée réduite n’a pas de 1-directeur.

24
Voici le pseudo-code de l’algorithme que nous avons esquissé. Dans ce pseudo-
code la matrice augmentée est notée A. Rappelons que cette matrice possède n
lignes et m + 1 colonne, puisqu’on doit compter une colonne pour chacune des m
inconnues et une colonne supplémentaire pour le membre de droite. Le coefficient
de ligne i et de colonne j de A est noté A(i, j). La variable l contient la ligne du
pivot (elle est initialisée à 1) et la variable c contient la colonne du pivot.
Poser l ← 1
Pour c qui va de 1 à m faire
Poser i ← l
Tant que i ≤ n et A(i, c) = 0 faire
Poser i ← i + 1
Fin tant que
Si i ≤ n alors
Échanger les lignes i et l de la matrice A
Diviser la ligne l de A par A(l, c)
Pour i qui va de 1 à n faire
Si i 6= l alors
Soustraire A(i, c) fois la ligne l à la ligne i dans la matrice A
Fin si
Fin pour
Poser l ← l + 1
Fin si
Fin pour

III Équations différentielles ordinaires

Une équation différentielle est une équation dont l’inconnue est une fonction. Ici il
s’agira de fonctions du type y : R → R. Ce type d’équation affirme que la dérivé
i-ème de y est une expression qui dépend de y et des (i − 1)-ièmes premières
dérivées de y. Par exemple :
y 0 = y.
Cette équation différentielle affirme que y et sa dérivée sont égales. Supposons
que la variable de la fonction y est notée t. Les solutions de cette équations sont

25
les fonctions de la forme y = Cet où C est une constrante réelle.

Objectif. Nous allons considérer des points t0 , t1 , ..., tN régulièrement disposés

sur l’axe horizontal, et pour chacun d’entre eux, nous allons obtenir une valeur
approchée ui de y(ti ).

III.1 Quelques exemples types

Saut en parachute.
Prenons l’exemple d’une personne qui fait un saut en parachute. Supposons que
la trajectoire de ce parachute est verticale. La personne ouvre le parachute au
temps t = 0 alors qu’il a une vitesse v(0). On suppose que la force qu’exerce l’air
sur le parachute est proportionnelle à la vitesse :

F (t) = −a v

Le bilan de forces sur le parachute produit l’équation suivante :

dv a
v0 = = − v + g.
dt m
Il s’agit d’une équation différentielle en v. Sa solution est la vitesse comme fonction
du temps.

Dynamique des populations.

• On considère un population de bactéries dans un environnement confiné,

• Pas plus de B individus ne peuvent coexister,
• Au temps initial, il y a y0 < B individus,
• Le taux de croissance de la population est une constante positive C,
• La vitesse de croissance de la population est proportionnelle au nombre de
bactéries sous contrainte que ce nombre ne peut dépasser B,
• On peut traduire cela par l’équation suivante :
dy  y
= Cy 1 − ,
dt B
dont la solution y(t) est le nombre de bactéries au temps t.
• Si on a deux populations antagonistes, l’une étant composée de prédateurs et

26
 l’autre de proies, on obtient un système de deux équations à deux inconnues :

dy1

 = C1 y1 (1 − b1 y1 − d2 y2 ) ,
dt


 dy2 = C2 y2 (1 − b2 y2 − d1 y1 ) .


dt
Ici y1 est la population de proies et y2 est la population de prédateurs, C1 et
C2 sont les taux de croissances de ces populations (indépendants de l’autre
population).

III.2 Le problème de Cauchy

Concentrons-nous sur les équations différentielles dites du premier ordre :
(
y 0 (t) = f (t, y(t)),
y(t0 ) = y0 .

On dit que y0 est la donnée initiale. Notons que le problème de Cauchy est un pro-
blème un peu plus fort que simplement l’équation différentielle, puisqu’il précise
la valeur de y(t0 ).

Exemples. Dans ce cas de l’équation y 0 = y, f (t, y) = y. Dans le cas du para-

chute, f (t, v) = −(a/m)v + g et v(0) est la donnée initiale.

Sous certaines conditions, la solution du problème de Cauchy existe et est unique :

Théorème. Supposons que la fonction f (t, y) est :

i. Continue par rapport à ses deux variables,
ii. Lipschitzienne de rapport L quelconque par rapport à sa deuxième variable,
c’est-à-dire :
|f (t, y1 ) − f (t, y2 )| ≤ L|y1 − y2 |

pour tout t, y1 , y2 .
Dans ce cas le problème de Cauchy admet une unique solution y et cette solution
est de classe C 2 .

• Malheureusement, il n’existe pas toujours de solution explicite à une telle

équation. Par exemple, la solution de y 0 = (y − t)/(y + t) vérifie la relation

27
 implicite
1
ln(t2 + y 2 ) + arctan(y/t) = C,
2
où C est une constante. Cependant, il n’existe pas de formule simple du type
y(t) = ... Il existe même des équations différentielles pour lesquelles on ne
2
peut pas trouver de relation implicite. Par exemple y 0 = e−t .
• Pour ces raisons, nous chercherons des solutions approchées.

III.3 Les méthodes d’Euler

Considérons un intervalle [t0 , T ] et découpons-le en N intervalles de même taille.
Ces sous-intervalles auront une longueur h = (T − t0 )/N appelée pas de discréti-
sation. Appelons ti la borne supérieure du i-ème de ces intervalles :

ti = t0 + i h.

Nous allons trouver une approximation ui de y(ti ). La suite des valeurs u0 , ...,
uN est appelée la solution humérique. Pour cela, une première méthode consiste
à construire u1 , ..., un itérativement avec u0 = y0 et :

ui+1 = ui + h f (ti , ui ). (7)

Ceci revient à supposer que y est confondue avec sa tangente en ti dans l’intervalle
[ti , ti+1 ]. Cette méthode s’appelle la méthode d’Euler explicite. Elle consiste à
calculer ui+1 à partir de ui par l’équation (7).

Le pseudo-code pour la méthode d’Euler explicite s’écrit :

Poser u0 = y0 , i = 0.
Tant que i < n, faire
Poser ui+1 = ui + h f (ti , ui ).
Incrémenter i de 1.

Exemple. Considérons l’équation différentielle ordinaire :

√
y 0 = 2 y,
√
c’est-à-dire f (t, y) = 2 y. Supposons qu’on veuille résoudre le problème de Cau-

28
chy avec t0 = 0 et y0 = 1. La solution de ce problème est :

y(t) = (t − 1)2 .

En effet, pour cette fonction,

(
y(t0 ) = y(0) = 1 = y0 ,
p p
y 0 (t) = 2(t − 1) = 2 (t − 1)2 = 2 y(t).

Appliquons la méthode d’Euler explicite avec h = 1 :

• On pose u0 = y0 = 1.
√
• i = 0 ; u1 = u0 + 2 u0 = 3.
√ √
• i = 1 ; u2 = u1 + 2 u1 = 3 + 2 3 ≈ 6.4641.
√ √ p √
• i = 2 ; u3 = u2 + 2 u1 = 3 + 2 3 + 2 3 + 2 3 ≈ 11.549.
Une autre méthode consiste à construire u1 , ..., un itérativement avec u0 = y0 et :

ui+1 = ui + h f (ti+1 , ui+1 ). (8)

Ceci revient à supposer que y est confondue avec sa tangente en ti+1 dans l’in-
tervalle [ti , ti+1 ].Cette méthode s’appelle la méthode d’Euler implicite. En effet,
l’équation (8) exprime ui+1 en fonction de ui mais aussi de lui-même ! En d’autres
termes, il faut trouver ui+1 comme une solution de cette équation.

Le pseudo-code pour la méthode d’Euler implicite s’écrit :

Poser u0 = y0 , i = 0.
Tant que i < n, faire
Résoudre ui+1 = ui + h f (ti+1 , ui+1 ).
Incrémenter i de 1.

Exemple. Considérons l’équation différentielle ordinaire :

y 0 = −y 2 ,

c’est-à-dire f (t, y) = −y 2 . Supposons qu’on veuille résoudre le problème de Cau-

chy avec t0 = 1 et y0 = 1. La solution de ce problème est :

y(t) = 1/t.

29
En effet, pour cette fonction,
(
y(t0 ) = y(1) = y0 = 1,
y 0 (t) = −1/t2 = −(1/t)2 = −(y(t))2 .

Appliquons la méthode d’Euler implicite avec h = 1 :

• On pose u0 = y0 = 1.
• i = 0 ; L’équation à résoudre est u1 = u0 − u21 , c’est-à-dire X 2 + X − 1 = 0.
Cette équation admet deux solutions :
√ √
−1 − 5 −1 + 5
X= ≈ −1.6180 et X = ≈ 0.6180.
2 2
On garde la solution la plus proche de u0 . C’est-à-dire u1 = 0.6180.
• i = 1 ; L’équation à résoudre est u2 = u1 −u22 , c’est-à-dire X 2 +X −0.6180 =
0. Cette équation admet deux solutions :
√ √
−1 − 3.4721 −1 + 3.4721
X= ≈ −1.4317 et X = ≈ 0.4317.
2 2
On garde la solution la plus proche de u1 . C’est-à-dire u1 = 0.4317.
• i = 2 ; L’équation à résoudre est u3 = u2 −u23 , c’est-à-dire X 2 +X −0.4317 =
0. Cette équation admet deux solutions :
√ √
−1 − 3.4721 −1 + 3.4721
X= ≈ −1.3256 et X = ≈ 0.3256.
2 2
On garde la solution la plus proche de u2 . C’est-à-dire u2 = 0.3256.

Remarque. La méthode d’Euler implicite est plus compliquée à mettre en oeuvre,

mais elle a de meilleures propriétés de stabilité que la méthode d’Euler explicite.

Une troisième méthode, dite méthode d’Euler modifiée consiste à supposer que y 0
est constante dans l’intervalle ]ti , ti+1 ], que cette constante est y 0 (ti + h/2) :

ui+1 = ui + hy 0 (ti + h/2).

Or, d’après notre équation différentielle, y 0 (ti + h/2) = f (ti + h/2, y(ti + h/2)),
ce qui donne :
ui+1 = ui + h f (ti + h/2, y(ti + h/2)).

30
Finalement, on estime y(ti + h/2) en confondant y avec sa tangente en ti dans
l’intervalle ]ti , ti + h/2] :
h
y(ti + h/2) = ui + f (ti , ui ).
2
Cette méthode consiste donc à construire u1 , ..., un itérativement avec u0 = y0
et :
h
ui+1 = ui + h f (ti + h/2, ui + f (ti , ui )). (9)
2
Cette méthode est plus précise que la méthode d’Euler explicite.

Le pseudo-code pour la méthode d’Euler modifiée s’écrit :

Poser u0 = y0 , i = 0.
Tant que i < n, faire
Poser k = ui + h f (ti , ui )/2.
Poser ui+1 = ui + h f (ti + h/2, k).
Incrémenter i de 1.

Exemple. Reprenons l’exemple précédent :

y 0 = −y 2 ,

c’est-à-dire f (t, y) = −y 2 . Supposons qu’on veuille résoudre le problème de Cau-

chy avec t0 = 1 et y0 = 1. Appliquons la méthode d’Euler modifiée avec h = 1 :
• On pose u0 = y0 = 1.
• i = 0 ; k = u0 + f (t0 , u0 )/2 = 1 − 1/2 = 1/2.
u1 = u0 + f (t0 + h/2, k) = u0 − 1/4 = 3/4 = 0.75.
• i = 1 ; k = u1 + f (t1 , u1 )/2 = 3/4 − 9/16 = 3/16.
u2 = u1 + f (t1 + h/2, k) = u1 − 9/256 = 3/4 − 9/256 ≈ 0.7148.

Exemple. Reprenons l’exemple utilisé pour la méthode d’Euler explicite :

√
y 0 = 2 y,
√
c’est-à-dire f (t, y) = 2 y. Supposons qu’on veuille résoudre le problème de Cau-
chy avec t0 = 0 et y0 = 1. Appliquons la méthode d’Euler modifiée avec h = 1 :
• On pose u0 = y0 = 1.

31
 • i = 0 ; k = u0 + f (t0 , u0 )/2 = 2.
√
u1 = u0 + f (t0 + h/2, k) = 1 + 2 2 ≈ 3.8284.
√ p √
• i = 1 ; k = u1 + f (t1 , u1 )/2 = 1 + 2 2 + 1 + 2 2 ≈ 5.7851.
√ √ √
q p
u2 = u1 + f (t1 + h/2, k) = 1 + 2 2 + 2 1 + 2 2 + 1 + 2 2 ≈ 8.6389.

III.4 Exemple d’application

Considérons le problème de Cauchy
(
y 0 (t) = sin(t) + y(t), t ∈]0, 1]
y(0) = 0.

• Dans ce cas, t0 = 0, y0 = 0 et f (t, y) = sin(t) + y.

• Commençons par étudier l’existence et l’unicité d’une solution à ce problème :
i. f est continue par rapport à ses deux variables parce que c’est une somme
de fonctions continues (t → sin(t) et y → y),
ii. Lipschitzienne de rapport par rapport L = 1 à sa deuxième variable. En
effet, pour tout réels y1 et y2 , et tout t ∈ [0, 1],

|f (t, y1 ) − f (t, y2 )| = |y1 − y2 | ≤ L|y1 − y2 |.

Par conséquent, d’après le théorème, notre problème de Cauchy admet

une unique solution y, et cette solution est de classe C 2 .
• La solution de ce problème de Cauchy est
1 t

y(t) = e − sin(t) − cos(t) .
2
En effet, dans ce cas, on a bien y(0) = 0 et
1 t
y 0 (t) =


e − cos(t) + sin(t)
2
1 t

= sin(t) + e − sin(t) − cos(t)
2
= sin(t) + y(t).

• Le graphe de y sur [0, 1] est le suivant :

32
 (1.0 ; 0.6683)
0.6

(0.8 ; 0.4057)
0.4

(0.6 ; 0.2161)
0.2
(0.4 ; 0.0907) (0.2 ; 0.0213)

0
0 0.2 0.4 0.6 0.8 1

Application de la méthode d’Euler explicite.

Ici t0 = 0 et T = 1. Prenons N = 5.
Dans ce cas, h = 0.2, t0 = 0, t1 = 0.2, t2 = 0.4, t3 = 0.6, t4 = 0.8, t5 = 1.0.

• Initialisation : u0 = y0 = 0 et i = 0
• i = 0 : u1 = u0 + 0.2(sin(t0 ) + u0 ) = 0.2 sin(0).
Donc u1 = 0.
• i = 1 : u2 = u1 + 0.2(sin(t1 ) + u1 ) = 0.2 sin(0.2)
Donc u2 ≈ 0.0397.
• i = 2 : u3 = u2 + 0.2(sin(t2 ) + u2 ) ≈ 0.0397 + 0.2(sin(0.4) + 0.0397)
Donc u3 ≈ 0.1256.
• i = 3 : u4 = u3 + 0.2(sin(t3 ) + u3 ) ≈ 0.1256 + 0.2(sin(0.6) + 0.1256)
Donc u4 ≈ 0.2636.
• i = 4 : u5 = u4 + 0.2(sin(t4 ) + u4 ) ≈ 0.2636 + 0.2(sin(0.8) + 0.2636)
Donc u5 ≈ 0.4598.

Application de la méthode d’Euler modifiée.

Utilisons les mêmes paramètres que pour la méthode d’Euler explicite :
h = 0.2, t0 = 0, t1 = 0.2, t2 = 0.4, t3 = 0.6, t4 = 0.8, t5 = 1.0

• Initialisation : u0 = y0 = 0 et i = 0
• i = 0 : k = u0 + 0.1(sin(t0 ) + u0 ) = 0.1 sin(0) = 0.
u1 = u0 + 0.2(sin(t0 + 0.1) + k) = 0.2 sin(0.1) ≈ 0.0200.

33
 • i = 1 : k = u1 + 0.1(sin(t1 ) + u1 ) ≈ 0.0200 + 0.1(sin(0.2) + 0.0200).
Donc k ≈ 0.0418
u2 = u1 + 0.2(sin(t1 + 0.1) + k) ≈ 0.0200 + 0.2(sin(0.3) + 0.0418).
Donc u2 ≈ 0.0874.
• i = 2 : k = u2 + 0.1(sin(t2 ) + u2 ) ≈ 0.0874 + 0.1(sin(0.4) + 0.0874).
Donc k ≈ 0.1351
u3 = u2 + 0.2(sin(t2 + 0.1) + k) ≈ 0.0874 + 0.2(sin(0.5) + 0.1351).
Donc u3 ≈ 0.2103.
• i = 3 : k = u3 + 0.1(sin(t3 ) + u3 ) ≈ 0.2103 + 0.1(sin(0.6) + 0.2103).
Donc k ≈ 0.2878
u4 = u3 + 0.2(sin(t3 + 0.1) + k) ≈ 0.2103 + 0.2(sin(0.7) + 0.2878).
Donc u4 ≈ 0.3968.
• i = 4 : k = u4 + 0.1(sin(t4 ) + u4 ) ≈ 0.3968 + 0.1(sin(0.8) + 0.3968).
Donc k ≈ 0.5082
u5 = u4 + 0.2(sin(t4 + 0.1) + k) ≈ 0.3968 + 0.2(sin(0.9) + 0.5082).
Donc u5 ≈ 0.6551.

Application de la méthode d’Euler implicite.

Dans le cadre de cette méthode, une étape préliminaire nécessaire consiste à écrire
l’équation à résoudre à l’étape i, donc la solution est ui+1 . Pour notre problème
de Cauchy :

ui+1 = ui + hf (ti+1 , ui+1 ) = ui + h(sin(ti+1 ) + ui+1 ).

Nous pouvons écrire une solution exacte de cette équation :

1
ui+1 = (ui + h sin(ti+1 )).
1−h
Il suffira donc de calculer ui+1 en fonction de ui , de h et de ti+1 en utilisant
l’expression ci-dessus. Attention, il est souvent nécessaire d’utiliser une méthode
numérique approchée pour résoudre l’équation de la méthode implicite !

Utilisons de nouveau les paramètres suivants :

h = 0.2, t0 = 0, t1 = 0.2, t2 = 0.4, t3 = 0.6, t4 = 0.8, t5 = 1.0

34
 • Initialisation : u0 = y0 = 0 et i = 0
• i = 0 : Calculons :
1 1
u1 = (u0 + 0.2 sin(0.2)) = sin(0.2) ≈ 0.0497.
0.8 4
• i = 1 : Calculons :
1 1
u2 = (u1 + 0.2 sin(0.4)) = (0.0497 + 0.2 sin(0.4)) ≈ 0.1594.
0.8 0.8
• i = 2 : Calculons :
1 1
u3 = (u2 + 0.2 sin(0.6)) = (0.1594 + 0.2 sin(0.6)) ≈ 0.3405.
0.8 0.8
• i = 3 : Calculons :
1 1
u4 = (u3 + 0.2 sin(0.8)) = (0.3405 + 0.2 sin(0.8)) ≈ 0.6049.
0.8 0.8
• i = 4 : Calculons :
1 1
u5 = (u4 + 0.2 sin(1)) = (0.6049 + 0.2 sin(1)) ≈ 0.9665.
0.8 0.8

Plaçons les approximations obtenues avec chacune de ces méthodes, ansi que le
graphe de la solution exacte sur la même figure :

(1.0 ; 0.6683)
0.6

(0.8 ; 0.4057)
0.4

(0.6 ; 0.2161)
0.2
(0.4 ; 0.0907) (0.2 ; 0.0213)

0
0 0.2 0.4 0.6 0.8 1

Sur cette figure, les approximations par les méthode d’Euler explicite, modifiée et
implicite sont respectivement en bleu, rouge et vert.
On observe que :
• La méthode modifiée est la plus précise,
• Les méthodes explicite et implicite donnent des erreurs comparables,
• La méthode explicite sous-estime y et la méthode explicite sur-estime y.

35
III.5 La méthode de Crank-Nicolson
Puisque les méthodes explicites et implicites donnent des erreurs opposées et
comparables en valeur absolue, il est naturel d’envisager une méthode mixte dans
laquelle l’itéré suivant ui+1 est la moyenne de ce que donneraient les méthodes
explicites et implicites, c’est-à-dire :
f (ti , ui ) + f (ti+1 , ui+1 )
ui+1 = ui + h .
2
Cette méthode s’appelle la méthode de Crank-Nicolson. Comme la méthode d’Eu-
ler implicite, elle nécessite d’obtenir ui+1 à chaque itération comme solution d’une
équation. Appliquons-là avec le problème de Cauchy précédent.

Commençons par écrire l’équation à résoudre à l’étape i, donc la solution est ui+1 .
Pour notre problème de Cauchy :
f (ti , ui ) + f (ti+1 , ui+1 ) h
ui+1 = ui + h = ui + (sin(ti ) + sin(ti+1 ) + ui + ui+1 ).
2 2
Nous pouvons écrire une solution exacte de cette équation :
  
2 h h h
ui+1 = 1+ ui + sin(ti ) + sin(ti+1 ) .
2−h 2 2 2
Il suffira donc de calculer ui+1 en fonction de ui , de h, de ti et de ti+1 en utilisant
l’expression ci-dessus. Attention, là aussi, il est généralement nécessaire d’utili-
ser une méthode numérique approchée pour résoudre l’équation de la méthode
implicite pour des fonctions f (t, y) plus compliquées !

Nous utiliserons de nouveau les paramètres suivants :

h = 0.2, t0 = 0, t1 = 0.2, t2 = 0.4, t3 = 0.6, t4 = 0.8, t5 = 1.0
• Initialisation : u0 = y0 = 0 et i = 0
• i = 0 : Calculons :
2
u1 = (1.1u0 + 0.1 sin(0) + 0.1 sin(0.2)).
1.8
0.2
Donc u1 = sin(0.2) ≈ 0.0221.
1.8
• i = 1 : Calculons :
2
u2 = (1.1u1 + 0.1 sin(0.2) + 0.1 sin(0.4)).
1.8
36
 2
Donc u2 = (1.1 · 0.0221 + 0.1 sin(0.2) + 0.1 sin(0.4)) ≈ 0.0923.
1.8
• i = 2 : Calculons :
2
u3 = (1.1u2 + 0.1 sin(0.4) + 0.1 sin(0.6)).
1.8
2
Donc u3 = (1.1 · 0.0923 + 0.1 sin(0.4) + 0.1 sin(0.6)) ≈ 0.2188.
1.8
• i = 3 : Calculons :
2
u4 = (1.1u3 + 0.1 sin(0.6) + 0.1 sin(0.8)).
1.8
2
Donc u4 = (1.1 · 0.2188 + 0.1 sin(0.6) + 0.1 sin(0.8)) ≈ 0.4099.
1.8
• i = 5 : Calculons :
2
u5 = (1.1u4 + 0.1 sin(0.8) + 0.1 sin(1.0)).
1.8
2
Donc u5 = (1.1 · 0.4099 + 0.1 sin(0.8) + 0.1 sin(1.0)) ≈ 0.6742.
1.8

Plaçons les approximations obtenues sur la figure précédente :

(1.0 ; 0.6683)
0.6

(0.8 ; 0.4057)
0.4

(0.6 ; 0.2161)
0.2
(0.4 ; 0.0907) (0.2 ; 0.0213)

0
0 0.2 0.4 0.6 0.8 1

Sur cette figure, les approximations par les méthode d’Euler explicite, modifiée et
implicite sont respectivement en bleu, rouge et vert, et les approximations par la
méthode de Crank-Nicolson sont en violet.

Références
[1] A. Quarteroni, F. Salieri, P. Gervasio, Calcul Scientifique. Cours, exercices
corrigés et illustrations en MATLAB et Octave, Springer, 2010.

37