Analyse Numérique

2ème Année Licence Mathématique

Probabilités et Statistique1

2021 − 2022

1 Department of Probabilities & Statistics.

Faculty of Mathematics.
University of Science and Technology Houari Boumediene.
BP 32 El-Alia, U.S.T.H.B, Algeria.
Table des matières

1 Introduction 3
1.1 Qu’est-ce que l’analyse numérique ? . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Qu’est-ce qu’un bon algorithme ? . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Pourquoi un ordinateur fait-il des calculs faux ? . . . . . . . . . . . . . . . 4
1.3.1 Ecriture en base b . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.2 Exemple 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.3 Exemple 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2 Résolution f (x) = 0 5
2.1 Résolution d’équations non linéaires . . . . . . . . . . . . . . . . . . . . . . 5
2.1.1 Racines de l’équation f (x) = 0 . . . . . . . . . . . . . . . . . . . . . 5
2.1.2 Séparation des racines . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Approximation des racines : Méthodes itératives . . . . . . . . . . . . . . . 6
2.2.1 Méthode de bisection (dichotomie) . . . . . . . . . . . . . . . . . . 7
2.2.2 Méthode de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . 8
2.2.3 Méthode du point fixe . . . . . . . . . . . . . . . . . . . . . . . . . 9

3 Interpolation Polynômial 13
3.1 Méthode de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.1.1 Détermination pratique des polynômes de Lagrange : . . . . . . . . 14
3.2 Méthode de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2.1 Calcul des différences divisées . . . . . . . . . . . . . . . . . . . . . 16

4 Intégration Numérique 19
4.1 Méthode des trapèzes (n = 1) . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.2 Méthode de Simpson (n = 2) . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.3 Méthode de Newton (n = 3) . . . . . . . . . . . . . . . . . . . . . . . . . . 23

5 Résolution d’Équations Différentielles Ordinaires (ÉDO) 25

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
5.2 Méthode d’Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
5.3 Méthodes de Runge-Kutta d’ordre 2 . . . . . . . . . . . . . . . . . . . . . . 28
5.4 Comparaison entre les deux méthodes . . . . . . . . . . . . . . . . . . . . . 29

i
ii TABLE DES MATIÈRES

6 TP Matlab 31
6.1 Résolution d’équations non linéaires . . . . . . . . . . . . . . . . . . . . . . 31
6.1.1 Méthode de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . 31
6.1.2 Méthode du point fixe . . . . . . . . . . . . . . . . . . . . . . . . . 32
6.2 Interpolation polynômial . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
6.2.1 Méthode de Lagrange et de Newton . . . . . . . . . . . . . . . . . . 33
6.2.2 Calcul des différences divisées . . . . . . . . . . . . . . . . . . . . . 34
6.3 Intégration numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
6.3.1 Méthode des trapèzes . . . . . . . . . . . . . . . . . . . . . . . . . . 35
6.3.2 Méthode de Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . 36
6.3.3 Méthode de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . 36
6.4 Résolution d’équations différentielles ordinaires (ÉDO) . . . . . . . . . . . 37
6.4.1 Méthode d’Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
6.4.2 Méthodes de Runge-Kutta d’ordre 2 . . . . . . . . . . . . . . . . . 40
1 Introduction

1.1 Qu’est-ce que l’analyse numérique ?

Les calculatrices et les ordinateurs nous permettent de faire beaucoup d’opérations et ce
très rapidement. Mais pour que les machines soient capables de faire ces calculs, il faut les
programmer !...et donc écrire les programmes. C’est l’objet essentiel de l’analyse numérique
qui s’est développée avec l’apparition des ordinateurs (et qui est alors une spécialité récente
des mathématiques).
Mais faire beaucoup d’opérations ne veut pas dire faire n’importe quoi : les méthodes
ont un "coût", lié d’une part au temps de calcul, i.e au nombre d’opérations élémentaires
(additions, soustractions, multiplications et divisions) à faire, et d’autre part à l’espace
mémoire nécessaire pour stocker les données et les résultats.
Par ailleurs,
√ l’ordinateur ne fait pas de calculs exacts (à cause du mode de représentation
des réels 2 ≈ 1.414214. Ceci constitue un inconvénient car cela provoque des erreurs
d’arrondis et de troncature.
L’objectif de l’analyse numérique est de développer des algorithmes, de les comparer
entre eux (c’est-à-dire d’estimer leurs performances a priori) et de les étudier afin de sélec-
tionner les bons algorithmes.
décrire le comportement de systèmes physiques pour prédire, optimiser, contrôler le
comportement de ces systèmes, mise en équation des phénomènes physiques observés

1.2 Qu’est-ce qu’un bon algorithme ?

Afin de choisir le meilleur algorithme possible, il faut pouvoir les comparer. Un bon
algorithme est alors un algorithme :
1. le moins coûteux possible en place mémoire,
2. le moins coûteux possible en temps de calcul : c’est-à-dire qui minimise le nombre
d’opérations nécessaires. C’est ce qu’on appelle un problème de complexité.
3. le plus stable possible, ie le moins sensible aux erreurs d’arrondi.
4. le plus précis possible : la solution que j’obtiens est une solution approchée, je veux
savoir si je suis près ou loin de la solution exacte. C’est ce qu’on appelle l’estimation
d’erreur.

3
4 CHAPITRE 1. INTRODUCTION

1.3 Pourquoi un ordinateur fait-il des calculs faux ?

Tout simplement parce qu’il ne connaît qu’un nombre fini de nombres ! Par exemple
ceux qui possèdent√un nombre fini donné de chiffres non nuls après la virgule, or ce n’est
pas le cas de 13 et 2 qui ont un nombre infini de chiffres non nuls après la virgule.

1.3.1 Ecriture en base b

La représentation des nombres que nous utilisons quotidiennement est la représentation
en base 10 mais les ordinateurs travaillent en base 2 (représentation binaire) ou en base 8
ou encore en base 16.

Définition 1.3.1. Soit b un entier strictement supérieur à 1. Pour tout entier n supérieur
ou égal à 1, il existe un unique entier p et des entiers di (0 ≤ i ≤ p) compris entre 0 et
b − 1 avec dp 6= 0 tels que :
p
b
di bi = dp dp−1 . . . d0
X
n=
i=0

1.3.2 Exemple 1
10
En base 10 : 234 = 2 × 102 + 3 × 10 + 4 × 100 ; donc 234 = 234
2
En base 2 : 234 = 27 + 26 + 25 + 23 + 2 ; donc 234 = 11101010
8
En base 8 : 234 = 3 × 82 + 5 × 8 + 2 ; donc 234 = 352
On généralise cette écriture aux nombres réels. Mais attention le nombre de chiffres
après la virgule étant infini la somme va de −∞ à p.

Définition 1.3.2. Soit b un entier strictement supérieur à 1. Pour tout réel x non nul, il
existe un unique entier p et des entiers di (i ≤ p) compris entre 0 et b − 1 avec dp 6= 0 tels
que :
p
b
di bi = dp dp−1 . . . d0 , d−1 . . . d−q . . .
X
x=
i=−∞

1.3.3 Exemple 2
10
1
En base 10 : 0.625 = 0 × 1010 + 6 × 10 + 2 × 1012 + 5 × 1013 ; donc 0.625 = 0.625
2
En base 2 : 0.625 = 0.500 + 0.125 = 0 × 210 + 1 × 21 + 0 × 212 + 1 × 213 ; donc 0.625 = 0.101
2 Résolution f (x) = 0

2.1 Résolution d’équations non linéaires

Le numéricien est souvent confronté à la résolution d’équations algébriques de la forme :

f (x) = 0 (2.1)

et ce dans toutes sortes de contextes. Introduisons dès maintenant la terminologie qui nous
sera utile pour traiter ce problème.

2.1.1 Racines de l’équation f (x) = 0

Définition 2.1.1. Soit f une fonction de R dans R dont le domaine de définition est une
partie Df de R. On dit que α ∈ Df est une racine de l’équation (2.1) si :

f (α) = 0 (2.2)

Résoudre l’équation (2.1) c’est trouver tous les nombres réels α tels que (2.2) soit vérifiée.
Théorème 2.1 (Valeurs Intermédiaires (TVI)). Soit I un intervalle. Soient a et b dans I
avec a < b. Soit f une application continue et dérivable dans I. Si : f (a)f (b) < 0 alors il
existe α ∈ [a, b] telque f (α) = 0.
Soit f (x) = x(1 + ex ) − ex et I = [0, 1]. f (x) est continue et dérivable sur I.

f (0) = −1 et f (1) = 1 ⇒ f (0)f (1) < 0

d’autre part f 0 (x) = 1 + xex ≥ 0 sur I (monotone). Donc d’après le théorème des valeurs
intermédiaires il existe une seule solution α ∈ [0, 1] telque f (α) = 0.

2.1.2 Séparation des racines

Définition 2.1.2. On dit qu’une racine α d’une fonction f (x) est séparable si on peut
trouver un intervalle [a, b] telque α soit la seule racine de cette équation dans [a, b], la
racine α est alors dite séparée (isolée).
Méthode graphique :

5
6 CHAPITRE 2. RÉSOLUTION F (X) = 0

1. Soit on trace le graphe de la fonction f (x) et on cherche son intersection avec l’axe
Ox.
2. Soit on décompse la fonction f (x) en deux fonctions f1 (x) et f2 (x) simples à étudier,
telle que :
f (x) = f1 (x) − f2 (x)
et on cherche les points d’intersection des graphes de f1 (x) et f2 (x), dont les abscisse
sont exactement les racines de l’équation f (x) = 0.
Soit la fonction : f (x) = ex + x. Alors on a :

f = ex
1 (x)
f (x) = f1 (x) − f2 (x) = ex − (−x), avec f2 (x) = −x
20
15

f1(x) = exp(x)
10
5

α●
0

f2(x) = − x
−5

−3 −2 −1 0 1 2 3

Figure 2.1 – Méthode graphique : f1 (x) ∩ f2 (x) = {α}, avec f (α)=0.

d’après la figure 2.1 : f1 (x) ∩ f2 (x) = {α} ⇒ α ∈ [−1, 0]. f (x) est continue et dérivable
sur [−1, 0]
f (−1) = −0, 63 et f (0) = 1 ⇒ f (0)f (1) < 0
d’autre part f 0 (x) = 1 + ex ≥ 0 sur [−1, 0] (monotone). Donc d’après le théorème des
valeurs intermédiaires il existe une seule solution α ∈ [−1, 0] telque f (α) = 0.

2.2 Approximation des racines : Méthodes itératives

Définition 2.2.1. On appelle méthode itérative un procédé de calcul de la forme :
xk+1 = F (xk ), k = 0, 1, 2, . . . (2.3)
2.2. APPROXIMATION DES RACINES : MÉTHODES ITÉRATIVES 7

dans lequel on part d’une valeur x0 (approchée) pour calculer x1 , puis à l’aide de x1 on
calculer x2 , ect. La formule (2.3) est dite formule de récurrence.

2.2.1 Méthode de bisection (dichotomie)

Nous nous intéressons à une fonction f continue sur l’intervalle I = [a, b] et telle
que f (a)f (b) < 0 ; en d’autres termes, nous avons encadré la racine. Il découle de ces
hypothèses que f s’annule au moins une fois dans I et c’est ce zéro que nous souhaitons
localiser précisément.

Définition 2.2.2. Le principe de la méthode est facile à comprendre. Soit [a0 , b0 ] un

intervalle dans lequel f (a0 )f (b0 ) < 0, on note c0 = (a0 + b0 )/2 le centre de l’intervalle. Si
f (a0 )f (c0 ) < 0, alors la racine α appartient à l’intervalle [a0 , c0 ]. On reprend le procédé
avec a1 = a0 et b1 = c0 . Sinon, c’est-à-dire si f (a0 )f (c0 ) > 0, on pose a1 = c0 et b1 = b0 . On
construit ainsi une suite d’intervalles emboîtés [ak , bk ] de longueur (a0 + b0 )/2k . Les suites
ak et bk convergent vers la solution α.

L’algorithme de bissection s’écrit :

Algorithme 2.2.1 (Méthode de bissection).

1. Étant donné l’intervalle I = [a, b], et  la valeur maximale de l’erreur (un critère
d’arrêt).
2. Effectuer :
x(k) = (a + b)/2
– Si f (a)f (x(k) ) < 0 alors b = x(k)
– sinon a = x(k)
3. – Si |f (x(k) )| <  écrire la solution α = x(k) .
– arrêté.
4. Si non, retour à l’étape 2.

L’erreur absolue de la méthode de dichotomie est :

b−a
2k+1
après k étapes car l’erreur est diminuée de moitié à chaque étape. Ainsi, la méthode
converge linéairement, ce qui est très lent par comparaison avec la méthode de Newton.
Si l’on se donne la tolérance relative (valeur maximale de l’erreur) , on sait majorer le
nombre d’itérations nécessaires pour satisfaire cette tolérance :
 
b−a
ln 
k≥
ln 2
Soit : f (x) = ex + x, I = [−1, 0] et  = 10−2 .
8 CHAPITRE 2. RÉSOLUTION F (X) = 0

k ≥ ln (1/10−2 ) / ln 2 = 6.643856 ≈ 7

k a b x(k) f (a) f (b) f (x(k) ) |f (x(k) )| < 10−2

1 −1 0 −0.5 −0.6321 1.0000 0.1065 FALSE
2 −1 −0.5 −0.75 −0.6321 0.1065 −0.2776 FALSE
3 −0.75 −0.5 −0.625 −0.2776 0.1065 −0.0897 FALSE
4 −0.625 −0.5 −0.5625 −0.0897 0.1065 0.0072 FALSE
5 −0.625 −0.5625 −0.5938 −0.0897 0.0072 −0.0415 FALSE
6 −0.5938 −0.5625 −0.5781 −0.0415 0.0072 −0.0171 FALSE
7 −0.5781 −0.5625 −0.5703 −0.0171 0.0072 −0.0049 TRUE

La solution de l’équation f (x) = 0 est α = x(7) ≈ −0, 5703.

2.2.2 Méthode de Newton-Raphson

La méthode de Newton-Raphson est l’une des méthodes les plus utilisées pour la réso-
lution des équations non linéaires (2.1).
Soit f : R → R une fonction dérivable. Soit x(0) un point donnée. On définit x(k+1)
comme étant le point ou cette droite intersecte l’axe Ox. On en déduit que :
 
f x(k)
x(k+1) = x(k) − , ∀ k = 0, 1, 2, . . . (2.4)
f 0 (x(k) )

A la k ème itération, la méthode de Newton-Raphson nécessite l’évaluation des deux fonc-

tions f et f 0 au point x(k) . Cet effort de calcul supplémentaire est plus que compensé par
une accélération de la convergence, la méthode de Newton-Raphson étant d’ordre 2.

Remarque 2.1. Choix de l’initialisation x0 , il doit vérifier la condition f (x0 )f 00 (x0 ) > 0.

Algorithme 2.2.2 (Méthode de Newton-Raphson).

1. Étant donné x(0) une valeur initiale de la solution, et  la valeur maximale de l’erreur
(un critère d’arrêt).
2. Effectuer :  
f x(k)
x(k+1) = x(k) − , ∀ k = 0, 1, 2, . . .
f 0 (x(k) )
f (x(k) )
3. Si x(k+1) − x(k) = f 0 (x(k) )
≤:
(k+1)
– écrire la solution α = x .
– arrêté.
4. Si non, retour à l’étape 2.
2.2. APPROXIMATION DES RACINES : MÉTHODES ITÉRATIVES 9

Soit : f (x) = ex + x, x(0) = 0 et  = 10−3 . On a f 0 (x) = 1 + ex , l’itération :

(k)
(k+1) (k) ex + x(k)
x =x − x(k) , ∀ k = 0, 1, 2, . . .
e +1

(k)
ex +x(k)
k x(k+1) (k) ≤ 10−3
ex +1
− x(0) = 0 −
0 x(1) = −0, 5000 0.5000
1 x(2) = −0, 5664 0.0664
2 x(3) = −0, 5672 0.0008 ≤ 10−3

La solution de l’équation f (x) est α = x(3) = −0, 5672.

2.2.3 Méthode du point fixe

Avant d’aborder les méthodes des points fixes, il importe de définir ce qu’est un point
fixe d’une fonction.

Définition 2.2.3. Un point fixe d’une fonction g(x) est une valeur de x qui reste invariante
pour cette fonction, c’est-à-dire toute solution de :

x = g(x)

est un point fixe de la fonction g(x).

Théorème 2.2 (théorème du point fixe). Supposons que g(x) est définie sur l’intervalle
[a, b] et satisfait aux conditions suivantes :
1. g([a, b]) ⊂ [a, b], i.e. ∀x ∈ [a, b], a ≤ g(x) ≤ b.
2. g(x) est contractante, i.e. ∃L ∈ R, 0 ≤ L ≤ 1 tel que :

∀x, y ∈ [a, b] : |g(x) − g(y)| ≤ L|x − y|

Alors g(x) admet un point fixe unique x∗ ∈ [a, b]. De plus, pour tout point x0 ∈ [a, b] la
suite xk définie par xk+1 = g(xk ) converge vers x∗ .

Remarque 2.2. La condition 2 du théorème du point fixe entraine :

1. La continuité de g dans [a, b].
2. En divisant l’inégalité par |x − y| et en passant à la limite quand y tend vers x, on
obtient : |g 0 (x)| ≤ L < 1, et ceci en supposant que g est dérivable au point x.
10 CHAPITRE 2. RÉSOLUTION F (X) = 0

Il existe un algorithme très simple permettant de déterminer des points fixes. Il suffit
en effet d’effectuer les itérations de la façon suivante :

x
0 donné
xn+1 = g (xn )

à partir d’une valeur estimée initiale x0 . L’intérêt de cet algorithme réside dans sa généralité
et dans la relative facilité avec laquelle on peut en faire l’analyse de convergence. Il en
résulte l’algorithme plus complet suivant.
Algorithme 2.2.3 (Méthode de point fixe).
1. Étant donné x0 une valeur initiale de la solution, et  la valeur maximale de l’erreur
(un critère d’arrêt).
2. Effectuer :
xk+1 = g(xk ), ∀ k = 0, 1, 2, . . .
3. Si |xk+1 −xk |
|xk+1 |
≤:
– écrire la solution α = xk+1 .
– arrêté.
4. Si non, retour à l’étape 2.
On peut résoudre des équations non linéaires de la forme f (x) = 0 en utilisant l’algo-
rithme des points fixes. Il suffit pour ce faire de transformer l’équation f (x) = 0 en un
problème équivalent de la forme x = g(x). L’ennui, c’est qu’il y a une infinité de façons
différentes de le faire.
Soit la fonction f (x) = ex + x avec x0 = −1 et  = 10−2 . On a :

f (x) = 0 ⇔ ex + x = 0 ⇔ x = −ex

donc on trouver : g(x) = −ex .

– La fonction g(x) est décroissant dans [−1, 0], on a :

−1 ≤ g(0) = −1 ≤ g(−1) = −0.3678 ≤ 0 ⇒ g([−1, 0]) ⊂ [−1, 0]

alors : ∃L ∈ [0, 1] tel que :

|g(b) − g(a)| ≤ L|b − a| ⇔ |g(0) − g(−1)| ≤ L

donc la fonction g(x) (voir la Figure 2.2) est contractante dans l’intervalle [−1, 0].
– On g 0 (x) = −ex , d’ou |g 0 (−1)| = 0.3678 < 1.
D’après le théoème de convergence du point fixe, notre itéation proposée x = − exp(x)
converge vers la solution unique x∗ de l’équation f (x) = ex + x = 0.
L’itération du point fixé est :

xk+1 = −exk , ∀ k = 0, 1, 2, . . .
2.2. APPROXIMATION DES RACINES : MÉTHODES ITÉRATIVES 11

g([− 1, 0]) ⊂ [− 1, 0]
La fonction g(x) est contractante
0.0

−0.2
g(x) = − exp(x)

−0.4

−0.6

−0.8

−1.0

−1.0 −0.8 −0.6 −0.4 −0.2 0.0

Figure 2.2 – g(x) = −ex .

|xk+1 −xk |
k xk+1 = −exk |xk+1 |
≤ 10−2
0 x1 = −0, 36787 1.71828
1 x2 = −0.69220 0.46853
2 x3 = −0.50047 0.38309
3 x4 = −0.60624 0.17446
4 x5 = −0.54540 0.11156
5 x6 = −0.57961 0.05903
6 x7 = −0.56012 0.03480
7 x8 = −0.57114 0.01930
8 x9 = −0.56488 0.01108
9 x10 = −0.56843 0.00624 ≤ 10−2

La solution de l’équation f (x) est α = x10 = −0.56843.

12 CHAPITRE 2. RÉSOLUTION F (X) = 0
3 Interpolation Polynômial

Soit y = f (x) une fonction dont on ne connait que les valeurs yi qu’elle prend aux
(n + 1) points distincts xi ; i = 0, 1, . . . , n. On a donc :

f (xi ) = yi ; i = 0, 1, . . . , n.

Problème : Déterminer un polynôme Pn (x) de degré inférieur ou égal à n tel que :

Pn (xi ) = yi = f (xi ); i = 0, 1, . . . , n.

de manière à pouvoir estimer les valeurs f (x) au moyen de Pn (x) pour x (xmin ≤ x ≤ xmax ).
C’est ce qu’on appelle : Interpoler la fonction f (x) par le polynôme Pn (x) au points
x0 , x1 , . . . , xn .

3.1 Méthode de Lagrange

Résolvons d’abord le problème partiel suivant : Construire un polynôme Li (x), de degré
n, tel que : 
1 si i = j
Li (xj ) = δij = 
0 si i 6= j
Le polynôme Li (x) à obtenir s’annule en n points x0 , x1 , . . . , xi−1 , xi+1 , . . . , xn . Il s’écrit
alors :

Li (x) = Ki (x − x0 ) . . . (x − xi−1 )(x − xi+1 ) . . . (x − xn ) ou Ki = cte (3.1)

Pour x = xi , on a :

Li (x) = Ki (xi − x0 ) . . . (xi − xi−1 )(xi − xi+1 ) . . . (xi − xn ) = 1

et donc :
1
Ki =
(xi − x0 ) . . . (xi − xi−1 )(xi − xi+1 ) . . . (xi − xn )
ou encore :
1
Ki = Qn
j=0 (xi − xj )
j6=i

13
14 CHAPITRE 3. INTERPOLATION POLYNÔMIAL

et dela : Qn
j=0 (x − xj )
j6=i
Li (x) = Qn avec degré de Li (x) égale à n
j=0 (xi − xj )
j6=i

ou :
(x − x0 ) . . . (x − xi−1 )(x − xi+1 ) . . . (x − xn )
Li (x) =
(xi − x0 ) . . . (xi − xi−1 )(xi − xi+1 ) . . . (xi − xn )
Pour chaque i = 0, 1, . . . , n, Li (x) est appelé : Polynôme de Lagrange au points xi .
Passons à présent à la résolution du problème général qui consiste à former Pn (x)
vérifiant les conditions indiquées plus haut ie : Pn (xi ) = yi . Ce polynôme est de la forme :
n
X
Pn (x) = yi Li (x)
i=0
Qn
n j=0 (x − xj )
yi Qnj6=i
X
=
i=0 j=0 (xi − xj )
j6=i
n
X (x − x0 ) . . . (x − xi−1 )(x − xi+1 ) . . . (x − xn )
= yi
i=0 (xi − x0 ) . . . (xi − xi−1 )(xi − xi+1 ) . . . (xi − xn )

Pn (x) s’appelle : Le polynôme d’interpolation de Lagrange.

Constuire le polynôme d’interpolaltion de Lagrange de la fonction y = f (x) = sin(πx)
pour x0 = 0, x1 = 1/6 et x2 = 1/2.
Solution : P2 (x) = 72 x − 3x2 (degré n ≤ 2)

3.1.1 Détermination pratique des polynômes de Lagrange :

La détermination des (n + 1) polynômes Li (x) peut être conduite de la façon suivante :

x − x0 x0 − x 1 x0 − x2 x0 − x3 · · · x0 − xn
x1 − x0 x − x1 x 1 − x2 x1 − x3 · · · x1 − xn
x2 − x0 x2 − x1 x − x 2 x2 − x3 · · · x2 − xn
.. .. .. .. ..
. . . . .
xn − x0 xn − x1 xn − x2 xn − x3 · · · x − xn

et on voit que :

Le produit des termes diagonaux du tableau

Li (x) =
Le produit des termes de la i ème linge du tableau

Déterminer les 5 polynômes Li (x) attachés aux valeurs : x0 = 0, x1 = 4, x2 = −4,

x3 = 2 et x4 = 1. Considérons le tableau ci-après :
3.1. MÉTHODE DE LAGRANGE 15

1.0
●

0.8
function(x) sin(pi * x)

0.6

●
0.4
0.2
0.0

● ● ●

0.0 0.2 0.4 0.6 0.8 1.0

Figure 3.1 – f (x) = sin(πx) en noir. P2 (x) = 27 x − 3x2 en rouge.

x −4 4 −2 −1
4 x−4 8 2 3
−4 −8 x + 4 −6 −5
2 −2 6 x−2 1
1 −3 5 −1 x − 1

x(x − 4)(x + 4)(x − 2)(x − 1) 3 7 3 1

L0 (x) = = −1 − x + x2 + x3 − x4
x(−4)4(−2)(−1) 2 16 32 32
x(x − 4)(x + 4)(x − 2)(x − 1) 1 5 1 3 1 4
L1 (x) = = x − x2 + x + x
4(x − 4)(8)(2)(3) 24 96 192 192
x(x − 4)(x + 4)(x − 2)(x − 1) 1 7 2 7 3 1 4
L2 (x) = =− x+ x − x + x
−4(−8)(x + 4)(−6)(−5) 120 480 960 960
x(x − 4)(x + 4)(x − 2)(x − 1) 2 2 1 1
L3 (x) = = − x + x2 + x3 − x4
2(−2)6(x − 2)1 3 3 24 24
x(x − 4)(x + 4)(x − 2)(x − 1) 32 16 2 1
L4 (x) = = − x − x2 − x3 + x4
1(−3)5(−1)(x − 1) 15 15 15 15
16 CHAPITRE 3. INTERPOLATION POLYNÔMIAL

3.2 Méthode de Newton

Définition 3.2.1. Soit f une fonction dont on connait les valeurs f (x0 ), f (x1 ),. . . ,f (xn )
qu’elle prend aux abscisses distinctes x0 , x1 , . . . , xn . On définit les différences divisées de f
aux points x0 , x1 , . . . , xn par les relations de récurrence :
δ(xi ) = f (xi )
δ(xi ) − δ(xi+1 )
δ(xi , xi+1 ) =
xi − xi+1
δ(xi , xi+1 ) − δ(xi+1 , xi+2 )
δ(xi , xi+1 , xi+2 ) = (3.2)
xi − xi+2
.. ..
. = .
δ(xi , xi+1 , . . . , xi+p−1 ) − δ(xi+1 , xi+2 , . . . , xi+p )
δ(xi , xi+1 , . . . , xi+p ) =
xi − xi+p
La dernière relation du ce système est appelée la différence divisée d’ordre p de la
fonction f aux points xi , xi+1 , . . . , xi+p
Propriété 3.2.1. La différence divisée δ(xi , xi+1 , . . . , xi+p ) est invariante par toute per-
mutation des abscisses xi , xi+1 , . . . , xi+p .

3.2.1 Calcul des différences divisées

Pour claculer la différence divisée d’ordre n de la fonction f aux points x0 , x1 , . . . , xn ,
on forme le tableau suivant en appliquant les formules (3.2) colonne après colonne :
xi f (xi ) δ1 δ2 ··· δn−1 δn
x0 f (x0 )
δ(x0 , x1 )
x1 f (x1 ) δ(x0 , x1 , x2 )
δ(x1 , x2 )
x2 f (x2 ) δ(x1 , x2 , x3 )
δ(x2 , x3 )
x3 f (x3 )
.. ..
. . ··· δ(x0 , . . . , xn−1 )
··· ··· δ(x0 , . . . , xn )
.. ..
. . ··· δ(x1 , . . . , xn )
xn−1 f (xn−1 ) δ(xn−2 , xn−1 , xn )
δ(xn−1 , xn )
xn f (xn )

Théorème 3.1. Le polynôme Pn (x) qui prend les valeurs f (xi ), i = 0, 1, . . . , n aux points
distincts xi , peut s’écrire :
Pn (x) = δ(x0 ) + δ(x0 , x1 )(x − x0 ) + δ(x0 , x1 , x2 )(x − x0 )(x − x1 ) + · · ·
· · · + δ(x0 , . . . , xn )(x − x0 )(x − x1 ) · · · (x − xn−1 ) (3.3)
Pn (x) s’appelle : Le polynôme d’interpolation de Newton.
3.2. MÉTHODE DE NEWTON 17

Lemme 3.1. Si Pn (x) est un polynôme de degré n, sa différence divisée d’ordre (n + 1)

est nulle, c’est à dire δ(x0 , . . . , xn ) = 0.

Remarque 3.1. L’équation (3.3) c’est la forme la plus utilisée pour calculer le polynôme
d’interpolation, car elle nécessite le moins de calculs pour obtenir numériquement ses co-
efficients.

Constuire le polynôme d’interpolaltion de Newton de la fonction y = f (x) = sin(πx)

pour x0 = 0, x1 = 1/6 et x2 = 1/2.
Solution : P2 (x) = 72 x − 3x2
Déterminer le polynôme d’interpolation de Newton satisfaisant au tableau ci-dessous :
Solution : P3 (x) = 53
30
x3 − 7x2 + 253
30
x−1

x 0 2 3 5
f (x) -1 2 9 87
18 CHAPITRE 3. INTERPOLATION POLYNÔMIAL
4 Intégration Numérique

On s’intéresse dans cette section à l’approximation de la fonctionnelle linéaire faisant

correspondre à une fonction f son intégrale définie sur un intervalle fixé [a, b] :
Z b
f → I(f ) = f (x)dx (4.1)
a

Si l’intégrale définie est connue, on utilise alors la formule de Newton-Leibnitz : ab f (x)dx =

R

F (b) − F (a) et tout se ramène au calcul d’une fonction F (primitive de f ) pour des valeurs
données de la variable.
Pourtant on aura intérèt à utiliser une méthode approchée dans le cas ou, bien que la
primitive soit connue, sa complexité la rend impropre au calcul numérique. On se propose
de chercher des approximations de I(f ) au moyen d’expression de la forme :
n
X
Ai f (xi )
i=0

4.1 Méthode des trapèzes (n = 1)

Posant x0 = a et x1 = b. L’approximation d’ordre 1 de I(f ) :
Z b 1
X
I(f ) = f (x)dx ' Ai f (xi )
a i=0
' A0 f (x0 ) + A1 f (x1 )

h
avec : A0 = A1 = 2
(h = b − a). Donc l’approximation d’ordre 1 est :
Z b
b−a
I(f ) = f (x)dx ' [f (a) + f (b)] (4.2)
a 2
La formule (4.2) est connue sous le nom de "Formule des Trapèzes".

Remarque 4.1. En général, la précision fournie par un seul trapèze n’est pas suffisante.
On divise alors le segment [a, b] en segment partiel égaux, et à chacun on applique la

19
20 CHAPITRE 4. INTÉGRATION NUMÉRIQUE

formule des trapèzes. Soient : x0 = a, x1 = x0 + h, x2 = x1 + 2h,. . . ,xn = xn−1 + nh. Avec

h = (b − a)/n, on trouve :
Z xn Z x1 Z x2 Z xn
f (x)dx ' f (x)dx + f (x)dx + · · · + f (x)dx
x0 x0 x1 xn−1
" #
f (x0 ) + f (xn )
'h + f (x1 ) + f (x2 ) + · · · + f (xn−1 )
2
ou encore :
" #
Z b
b − a f (a) + f (b)
f (x)dx ' + f (x1 ) + f (x2 ) + · · · + f (xn−1 ) (4.3)
a n 2
La formule (4.3) est connue sous le nom de "Formule des Trapèzes Générale (Com-
posée)".
L’erreur globale de l’approximation (4.3) étant donnée par :
(b − a) 2
|E(h)| ≤ h max |f 00 (x)|
12 x∈[a,b]
3
(b − a)
≤ max |f 00 (x)|
12n2 x∈[a,b]
La méthode des trapèzes générale (composée) est d’ordre 2.
Remarque 4.2. Pour appliquer la méthode des trapèzes au calcul d’une intégrale avec
une erreur inférieur à , on choisit le nombre n de trapèzes grâce à la majoration :
(b − a)3
2
max |f 00 (x)| ≤ 
12n x∈[a,b]
Calculer l’intégrale : Z 8 √
I1 = x x2 + 1dx
0
1. Par la méthode des trapèzes.
2. Par la méthode des trapèzes, forme générale, avec n = 4.
Solution :
√
1. I1 = 08 x x2 + 1dx ' 257, 9922
R
√
2. I1 = 08 x x2 + 1dx ' 179, 4203
R

3. Solution exacte à 4 chiffres 174, 3489.

Calculons une valeur approchée de 01 exp(−x2 ) avec une erreur inférieur à 10−5 , en
R

utilisant la formule des trapèzes ?

Pour cette fonction, on peut montrer que maxx∈[0,1] |f 00 (x)| = 2, si on utilise la formule
des trapèzes, on doit donc choisir n tel que :
2
2
≤ 10−5 ⇔ n ≥ 130
12n
4.2. MÉTHODE DE SIMPSON (N = 2) 21

4.2 Méthode de Simpson (n = 2)

a+b
Posons x0 = a, x1 = 2
et x2 = b. L’approximation d’ordre 2 de I(f ) :
Z x2 2
X
I(f ) = f (x)dx ' Ai f (xi )
x0 i=0
' A0 f (x0 ) + A1 f (x1 ) + A2 f (x2 )

avec : Ai = I(Li ) = xx02 Li (x)dx, ou Li est le polynôme de Lagrange au point xi (∀ i =

R

0, 1, 2).
Posons h = x1 − x0 = x2 − x1 = (b − a)/2. Afin de calculer les polynômes de Lagrange
on considère le tableau suivant :

x − x0 x0 − x1 x0 − x2
x1 − x0 x − x1 x1 − x2
x2 − x0 x2 − x1 x − x2
on a alors :

(x − x1 )(x − x2 ) (x − x1 )(x − x2 )
L0 (x) = =
(x0 − x1 )(x0 − x2 ) 2h2
(x − x0 )(x − x2 ) (x − x0 )(x − x2 )
L1 (x) = =−
(x1 − x0 )(x1 − x2 ) h2
(x − x0 )(x − x1 ) (x − x0 )(x − x1 )
L2 (x) = =
(x2 − x0 )(x2 − x1 ) 2h2

Il vient, après calcul :

Z x2
1 Z x2 h
A0 = I(L0 ) = L0 (x)dx = 2 (x − x1 )(x − x2 )dx =
x0 2h x0 3
Z x2
1 Z x2 4h
A1 = I(L1 ) = L1 (x)dx = − 2 (x − x0 )(x − x2 )dx =
x0 h x0 3
Z x2
1 Z x2
h
A2 = I(L2 ) = L2 (x)dx = 2 (x − x0 )(x − x1 )dx =
x0 2h x0 3

Donc l’approximation d’ordre 2 est :

I(f ) ' A0 f (x0 ) + A1 f (x1 ) + A2 f (x2 )

h 4h h
' f (x0 ) + f (x1 ) + f (x2 )
3 3 3
h
' [f (x0 ) + 4f (x1 ) + f (x2 )]
3
22 CHAPITRE 4. INTÉGRATION NUMÉRIQUE

avec h = (b − a)/2 et c = (a + b)/2, on a lors :

Z b
b−a
I(f ) = f (x)dx ' [f (a) + 4f (c) + f (b)] (4.4)
a 6
La formule (4.4) est connue sous le nom de "Formule de Simpson".
Remarque 4.3. Comme dans le cas de la méthode des trapèzes, il y a souvent intérêt
à décomposer au préalable l’intervalle total e, intervalles partiels égaux, et à appliquer
séparement la méthode de Simpson à chacun d’eux. Soient, en posant n = 2k, x0 = a,
x1 = x0 + h, x2 = x1 + h,. . . ,xn = xn−1 + h. Avec h = (b − a)/n on trouve :
Z xn Z x1 Z x2 Z xn
f (x)dx = f (x)dx + f (x)dx + · · · + f (x)dx
x0 x0 x1 xn−1
hh
' f (x0 ) + f (xn ) + 4 (f (x1 ) + f (x3 ) + f (x5 ) + · · · + f (xn−1 ))
3 i
+ 2 (f (x2 ) + f (x4 ) + f (x6 ) + · · · + f (xn−2 ))

ou encore :
n n
−1
 
Z b
b−a X2 2
X
f (x)dx ' f (a) + f (b) + 4 f (x2i−1 ) + 2 f (x2i ) (4.5)
a 3n i=1 i=1

Tous les termes de rang impair sont multipliés par 4 tandis que ceux de rang pair sont
multipliés par 2, sauf le premier f (a) et le dernier f (b). La formule (4.5) est connue sous
le nom de "Formule de Simpson Générale".
L’erreur globale de l’approximation (4.5) étant donnée par :
(b − a) 4
|E(h)| ≤ h max |f (4) (x)|
180 x∈[a,b]
5
(b − a)
≤ max |f (4) (x)|
180n4 x∈[a,b]
La méthode des Simpson générale est d’ordre 4.
Remarque 4.4. Pour déterminer le nombre (suffisant) de sous intervalles partiels de l’in-
tervalle d’intégration tel que l’erreur soit inférieur à , il suffit de faire la majoration :
(b − a)5
max |f (4) (x)| ≤ 
180n4 x∈[a,b]
Calculer l’intégrale : Z 8 √
I1 = x x2 + 1dx
0
1. Par la méthode de Simpson.
2. Par la méthode de Simpson, forme générale, pour n = 4(k = 2) et n = 6(k = 3).
4.3. MÉTHODE DE NEWTON (N = 3) 23

Solution :
√
1. I1 = 08 x x2 + 1dx ' 173, 9570
R
√
2. ◦ I1 = R08 x√x2 + 1dx ' 174, 2385 pour n = 4.
R

◦ I1 = 08 x x2 + 1dx ' 174, 3079 pour n = 6.

3. Solution exacte à 4 chiffres 174, 3489.
Calculons une valeur approchée de 01 exp(−x2 ) avec une erreur inférieur à 10−5 , en
R

utilisant la formule de Simpson ?

Pour cette fonction, on peut montrer que maxx∈[0,1] |f (4) (x)| = 12, si on utilise la formule
de Simpson, on doit donc choisir n tel que :
12
≤ 10−5 ⇔ n ≥ 10
180n4

4.3 Méthode de Newton (n = 3)

Soient x0 = a, x1 = x0 + h, x2 = x0 + 2h, x3 = x0 + 3h = b, avec h = (b − a)/3.
L’approximation d’order 3 de la fonctionnelle I(f ) est :
Z x3 3
X
I(f ) = f (x)dx ' Ai f (xi )
x0 i=0
' A0 f (x0 ) + A1 f (x1 ) + A2 f (x2 ) + A3 f (x3 )
R x3
avec : Ai = I(Li ) = x0 Li (x)dx, ou Li est le polynôme de Lagrange au point xi (∀ i =
0, 1, 2, 3).
Après calcul, il vient :
3 9 9 3
A0 = h, A1 = h, A2 = h, A3 = h
8 8 8 8
et alors :
Z x3
3
I(f ) = f (x)dx ' h [f (x0 ) + 3f (x1 ) + 3f (x2 ) + f (x3 )]
x0 8
avec h = (b − a)/3, c = (2a + b)/3 et d = (a + 2b)/3, on a alors :
Z b
b−a
I(f ) = f (x)dx '
[f (a) + 3f (c) + 3f (d) + f (b)] (4.6)
a 8
La formule (4.6) est connue sous le nom de "Formule de Newton".
Calculer l’intégrale : Z 8 √
I1 = x x2 + 1dx
0
1. Par la méthode de Newton.
Solution :
√
1. I1 = 08 x x2 + 1dx ' 174, 1024
R

2. Solution exacte à 4 chiffres 174, 3489.

24 CHAPITRE 4. INTÉGRATION NUMÉRIQUE
5 Résolution d’Équations Différen-
tielles Ordinaires (ÉDO)

5.1 Introduction
Nous prenons comme point de départ la formulation générale d’une équation différen-
tielle d’ordre 1 avec condition initiale. La tâche consiste à déterminer une fonction y(t)
solution de problème de Cauchy :

y 0 = f (t, y)
y (t0 )
(5.1)
= y0
La variable indépendante t représente très souvent (mais pas toujours) le temps. La
variable dépendante est notée y et dépend bien sûr de t. La fonction f est pour le moment
une fonction quelconque de deux variables que nous supposons suffisamment différentiable.
La condition y (t0 ) = y0 est la condition initiale et en quelque sorte l’état de la solution
au moment où on commence à s’y intéresser. Il s’agit d’obtenir y(t) pour t ≥ t0 , si on
cherche une solution analytique, ou une approximation de y(t), si on utilise une méthode
numérique.
Le théorème ci-après nous donne des conditions qui assurent l’existence et l’unicité de
la solution y(t) théorique (explicite) de problème (5.1).
Théorème 5.1. Si f est continue et s’il existe une constante L strictement positive telle
que pour tout t ∈ [t0 , t0 + T ], et pour tout y1 , y2 ∈ R, on ait :

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

alors, le problème de Cauchy admet une solution unique, quelque soit y0 ∈ R. On dit alors
que f est L− lipschitzienne, ou L est la constante de Lipschitz.
Remarque 5.1. Nous nous plaçerons toujours dans les conditions du théorème 5.1.
Soit l’équation différentielle du premier ordre :

y 0 (t) =t
(5.2)
y (0) =1

25
26CHAPITRE 5. RÉSOLUTION D’ÉQUATIONS DIFFÉRENTIELLES ORDINAIRES (ÉDO)

Voilà certainement l’un des exemples les plus simples que l’on puisse imaginer. La solution
particulière est alors :
t2
y(t) = + 1 (5.3)
2
qui vérifie à la fois l’équation différentielle et la condition initiale.
Soit l’équation différentielle :

y 0 (t)
= ty(t)
y (0) = 2
(5.4)

ne suffit pas dans ce cas d’intégrer les deux côtés de l’équation pour obtenir la solution.
On doit d’abord séparer les variables en écrivant par exemple :

dy
= ty(t)
dt
qui devient :
dy
= tdt
y
la solution particulière est :
2 /2
y(t) = 2et (5.5)
qui vérifie à la fois l’équation différentielle et la condition initiale.
L’objectif de ce chapitre est de décrire un certain nombre de méthodes permettant de
résoudre numériquement le problème de Cauchy (5.1).
Étant donné une subdivision t0 < t1 < · · · < tN = t0 + T de l’intervalle [t0 , t0 + T ],
on cherche à déterminer des valeurs approchées y0 , y1 , . . . , yN des valeurs ytn prises par la
solution exacte y(t). On notera les pas succéssives :

h = tn+1 − tn , ∀ 0≤n≤N −1

Dans toutes les méthodes numériques développées par la suite, on subdivise l’intervalle
[t0 , t0 + T ] en N intervalles de longueur h donnée par :

T − t0
h=
N

5.2 Méthode d’Euler

Algorithme 5.2.1 (Méthode d’Euler).
1. Étant donné un pas de temps h, une condition initiale (t0 , y0 ) et un nombre maximal
d’itérations N .
2. Pour 0 ≤ n ≤ N − 1 :
5.2. MÉTHODE D’EULER 27

n+1 = yn + hf (tn , yn )
y
– 
tn+1 = tn + h
– Écrire tn+1 et yn+1
3. Arrêt.

Soit le problème de Cauchy :


y 0 (t)
= ty(t)
(5.6)
y (0) = 2

On veut approcher la solution de (5.6) en t = 1 à l’aide de la méthode d’Euler, en subdi-

visant l’intervalle [0, 1] en 10 parties égales.
Selon l’agorithme d’Euler :

y 0

 = 2; t0 = 0; T = 1; N = 10; h = (T − t0 )/N = 0, 1

yn+1 = yn + htn yn , ∀ 0 ≤ n ≤ 9 (5.7)


tn+1 = tn + h

On calculer les valeurs :

n 0 1 2 3 4 5 6 7 8 9 10
tn 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
yn 2,0000 2,0000 2,0200 2,0604 2,1222 2,2071 2,3175 2,4565 2,6285 2,8387 3,0942
y(tn ) 2,0000 2,0100 2,0404 2,0921 2,1666 2,2663 2,3944 2,5552 2,7543 2,9986 3,2974

On trouve y(1) ≈ 3, 0942. La solution exacte de l’équation (5.6) est donnée par l’équa-
tion (6.2), ce qui donne y(1) = 3, 2974. L’approximation calculée est donc très grossière.

Figure 5.1 – h = 0, 05. Figure 5.2 – h = 0, 025.

28CHAPITRE 5. RÉSOLUTION D’ÉQUATIONS DIFFÉRENTIELLES ORDINAIRES (ÉDO)

Figure 5.3 – h = 0, 01. Figure 5.4 – h = 0, 001.

5.3 Méthodes de Runge-Kutta d’ordre 2

Carl Runge (1856-1927) et Martin Kutta (1867-1944) ont proposé en 1895 de résoudre
le problème de Cauchy (5.1). en introduisant un schéma numérique de la forme :

y
n+1 = yn + hΦ(tn , yn )
tn+1 = tn + h

où la fonction Φ est une approximation de f (t, y) sur l’intervalle [tn , tn+1 ]. La méthode
Runge-Kutta est plus symétrique puisqu’elle revient à évaluer numériquement la dérivée
au centre de l’intervalle [tn , tn+1 ]. Ceci conduit à l’algorithme suivant :

Algorithme 5.3.1 (Méthode de Runge-Kutta d’ordre 2).

1. Étant donné un pas de temps h, une condition initiale (t0 , y0 ) et un nombre maximal
d’itérations N .
2. Pour 0 ≤ n ≤ N − 1 :




k1 = f (t ,y )
n n 
k = f t + h , y + h k


2 n 2 n 2 1
–




yn+1 = yn + hk2


tn+1 = tn + h
– Écrire tn+1 et yn+1
3. Arrêt.

On veut approcher la solution de (5.6) en t = 1 à l’aide de la méthode de Runge-Kutta

d’ordre 2, en subdivisant l’intervalle [0, 1] en 10 parties égales.
5.4. COMPARAISON ENTRE LES DEUX MÉTHODES 29

Selon l’agorithme de Runge-Kutta d’ordre 2 :




 y0 = 2; t0 = 0; T = 1; N = 10; h = (T − t0 )/N = 0, 1

k1 = tn yn






k2 = 0.5yn (h + 2 (1 + 0.25h2 ) tn + ht2n ) (5.8)
2 2

yn+1 = yn [1 + 0.5h (h + 2 (1 + 0.25h ) tn + htn )] , ∀ 0 ≤ n ≤ 9






tn+1 = tn + h
On calculer les valeurs :
n 0 1 2 3 4 5 6 7 8 9 10
tn 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
yn 2,0000 2,0100 2,0403 2,0918 2,1661 2,2656 2,3933 2,5535 2,7517 2,9950 3,2923
y(tn ) 2,0000 2,0100 2,0404 2,0921 2,1666 2,2663 2,3944 2,5552 2,7543 2,9986 3,2974

On trouve y(1) ≈ 3, 2923. La solution exacte de l’équation (5.6) est donnée par l’équa-
tion (6.2), ce qui donne y(1) = 3, 2974. L’approximation calculée est donc très bonne.

Figure 5.5 – Runge-Kutta d’ordre 2 : y 0 (t) = ty(t) avec y(0) = 2 et h = 0, 1.

5.4 Comparaison entre les deux méthodes

n 0 1 2 3 4 5 6 7 8 9 10
tn 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
Euler 2,0000 2,0000 2,0200 2,0604 2,1222 2,2071 2,3175 2,4565 2,6285 2,8387 3,0942
RK2 2,0000 2,0100 2,0403 2,0918 2,1661 2,2656 2,3933 2,5535 2,7517 2,9950 3,2923
Exacte 2,0000 2,0100 2,0404 2,0921 2,1666 2,2663 2,3944 2,5552 2,7543 2,9986 3,2974
30CHAPITRE 5. RÉSOLUTION D’ÉQUATIONS DIFFÉRENTIELLES ORDINAIRES (ÉDO)
 6 TP Matlab

6.1 Résolution d’équations non linéaires

6.1.1 Méthode de Newton-Raphson

1 f u n c t i o n [ xvect , x d i f , fx , n i t ]= newton ( x0 , t o l , fun , dfun )

2 e r r=t o l +1; n i t =0; x v e c t=x0 ; x=x0 ; f x=e v a l ( fun ) ; x d i f = [ ] ;
3 w h i l e e r r >t o l
4 n i t=n i t +1;
5 x=x v e c t ( n i t ) ;
6 d f x=e v a l ( dfun ) ;
7 i f d f x==0
8 e r r=t o l ∗ 1 . e −10;
9 f p r i n t f ( ’ a r r e t c a r dfun e s t n u l l e ’ ) ;
10 else
11 xn=x−f x ( n i t ) / d fx ; e r r=abs ( xn−x ) ; x d i f =[ x d i f ; e r r ] ;
12 x=xn ; x v e c t =[ x v e c t ; x ] ; f x =[ f x ; e v a l ( fun ) ] ;
13 end
14 end
15 return
Soit : f (x) = ex + x et f 0 (x) = 1 + ex , avec : x(0) = 0 et  = 10−3 .
1 >> f = ’ exp ( x )+x ’ ;
2 >> d f = ’1+exp ( x ) ’ ;
3 >> [ xvect , x d i f , fx , n i t ]= newton (0 ,10^ −3 , f , d f )
4

5 xvect =
6 0
7 −0.5000
8 −0.5663
9 −0.5671
10

11 xdif =
12 0.5000

31
 32 CHAPITRE 6. TP MATLAB

13 0.0663
14 0.0008
15

16 fx =
17 1.0000
18 0.1065
19 0.0013
20 0.0000
21

22 nit =
23 3
La solution de l’équation f (x) est α = x(3) = −0, 5671.
6.1.2 Méthode du point fixe

1 f u n c t i o n [ xvect , x d i f , fx , n i t ]= f i x p o i n t ( x0 , t o l , fun , p h i )
2 e r r=t o l +1; n i t =0;
3 x v e c t=x0 ; x=x0 ; f x=e v a l ( fun ) ; x d i f = [ ] ;
4 while err > t o l
5 n i t=n i t +1;
6 x=x v e c t ( n i t ) ;
7 xn=e v a l ( ph i ) ;
8 e r r=abs ( xn−x ) / abs ( xn ) ;
9 x d i f =[ x d i f ; e r r ] ;
10 x=xn ; x v e c t =[ x v e c t ; x ] ; f x =[ f x ; e v a l ( fun ) ] ;
11 end
12 return
Soit la fonction f (x) = ex + x avec x0 = −1 et  = 10−2 . On a :

f (x) = 0 ⇔ ex + x = 0 ⇔ x = −ex

donc on trouver : g(x) = −ex

1 >> f = ’ exp ( x )+x ’ ;
2 >> p h i = ’−exp ( x ) ’ ;
3 >> [ xvect , x d i f , fx , n i t ]= f i x p o i n t ( −1 ,10^ −2 , f , ph i )
4

5 xvect =
6 −1.0000
7 −0.3679
8 −0.6922
9 −0.5005
10 −0.6062
11 −0.5454
 6.2. INTERPOLATION POLYNÔMIAL 33

12 −0.5796
13 −0.5601
14 −0.5711
15 −0.5649
16 −0.5684
17

18 xdif =
19 1.7183
20 0.4685
21 0.3831
22 0.1745
23 0.1116
24 0.0590
25 0.0348
26 0.0193
27 0.0111
28 0.0062
29

30 fx =
31 −0.6321
32 0.3243
33 −0.1917
34 0.1058
35 −0.0608
36 0.0342
37 −0.0195
38 0.0110
39 −0.0063
40 0.0035
41 −0.0020
42

43 nit =
44 10
La solution de l’équation f (x) est α = x(10) = −0.5684.
6.2 Interpolation polynômial
6.2.1 Méthode de Lagrange et de Newton
Utiliser la formule d’interpolation (Lagrange ou celle de Newton) pour trouver le poly-
nôme passant par : x0 = 0, 4 ; x1 = 0, 5 ; x2 = 0, 7 ; x3 = 0, 8 pour la fonction f (x) = sin(x).
1 >> x = [ 0 . 4 0 . 5 0 . 7 0 . 8 ] ;
2 >> y = s i n ( x ) ;
 34 CHAPITRE 6. TP MATLAB

3 >> p = p o l y f i t ( x , y , l e n g t h ( x ) −1)
4

5 p =
6 −0.1372 −0.0342 1.0145 −0.0021
7

8 >> poly2sym ( p , ’ x ’ )
9

10 ans =
11 3 2
12 −0.0021 − 0 . 1 3 7 2 x − 0.0342 x + 1.0145 x
13 >>
On trouve : P3 (x) = −0.1372x3 − 0.0342x2 + 1.0145x − 0.0021.

6.2.2 Calcul des différences divisées

1 f u n c t i o n [ d]= d i v i d i f ( x , y )
2 [ n ,m]= s i z e ( y ) ;
3 i f n == 1 , n = m; end
4 n = n−1;
5 d = z e r o s ( n+1,n+1) ;
6 d(: ,1) = y ’;
7 f o r j = 2 : n+1
8 f o r i = j : n+1
9 d ( i , j ) = ( d ( i −1, j −1)−d ( i , j −1) ) / ( x ( i −j +1)−x ( i ) ) ;
10 end
11 end
12 d=[x ’ , d ] ;
13 return
Claculer les différences divisées de la fonction f (x) aux points xi :

x 0 2 3 5
f (x) -1 2 9 87

1 >> x = [ 0 2 3 5 ] ;
2 >> y = [−1 2 9 8 7 ] ;
3 >> Tab=d i v i d i f ( x , y )
4

5 Tab =
6 0 −1.0000 0 0 0
7 2.0000 2.0000 1.5000 0 0
8 3.0000 9.0000 7.0000 1.8333 0
9 5.0000 87.0000 39.0000 10.6667 1.7667
 6.3. INTÉGRATION NUMÉRIQUE 35

6.3 Intégration numérique

6.3.1 Méthode des trapèzes

1 f u n c t i o n [ i n t ] = t r a p e z ( a , b ,m, fun )
2 h=(b−a ) /m;
3 x=[a : h : b ] ;
4 dim=l e n g t h ( x ) ;
5 y=e v a l ( fun ) ;
6 i f s i z e ( y )==1
7 y=d i a g ( ones ( dim ) ) ∗y ;
8 end
9 i n t=h ∗ ( 0 . 5 ∗ y ( 1 )+sum ( y ( 2 :m) ) +0.5∗ y (m+1) ) ;
10 return
Calculer l’intégrale : Z 8 √
I(f ) = x x2 + 1dx
0
1. Par la méthode des trapèzes n = 1.
2. Par la méthode des trapèzes, forme générale, avec n = 4.

1 >> f = ’ x . ∗ s q r t (1+x . ^ 2 ) ’ ;
2 >> I 1 = t r a p e z ( 0 , 8 , 1 , f )
3

4 I1 =
5 257.9922
6

7 >> I 2 = t r a p e z ( 0 , 8 , 4 , f )
8

9 I2 =
10 179.4203
R8 √
1. I1 = 0 x x2 + 1dx ' 257, 9922
R8 √
2. I1 = 0 x x2 + 1dx ' 179, 4203

1 >> syms x ;
2 >> do ub le ( i n t ( x∗ s q r t (1+x ^2) , 0 , 8 ) )
3

4 ans =
5

6 174.3489

1. Solution exacte à 4 chiffres 174, 3489.

 36 CHAPITRE 6. TP MATLAB

6.3.2 Méthode de Simpson

1 f u n c t i o n [ i n t ] = simpson ( a , b ,m, fun )

2 h=(b−a ) /m;
3 x=[a : h / 2 : b ] ;
4 dim= l e n g t h ( x ) ;
5 y=e v a l ( fun ) ;
6 i f s i z e ( y )==1
7 y=d i a g ( ones ( dim ) ) ∗y ;
8 end
9 i n t =(h / 6) ∗( y ( 1 ) +2∗sum ( y ( 3 : 2 : 2 ∗m−1) ) +4∗sum ( y ( 2 : 2 : 2 ∗m) )+y
(2∗m+1) ) ;
10 return
Calculer l’intégrale :
Z 8 √
I(f ) = x x2 + 1dx
0

1. Par la méthode de Simpson.

2. Par la méthode de Simpson, forme générale, pour n = 4(k = 2) et n = 6(k = 3).

1 >> f = ’ x . ∗ s q r t (1+x . ^ 2 ) ’ ;
2 >> I 1=simpson ( 0 , 8 , 1 , f )
3

4 I1 =
5 173.9570
6

7 >> I 2=simpson ( 0 , 8 , 2 , f )
8

9 I2 =
10 174.2385
11

12 >> I 3=simpson ( 0 , 8 , 3 , f )
13

14 I3 =
15 174.3079
R8√
1. I(f ) = x x2 + 1dx ' 173, 9570
0
√
2. ◦ I(f ) = R08 x√x2 + 1dx ' 174, 2385 pour n = 4(k = 2).
R

◦ I(f ) = 08 x x2 + 1dx ' 174, 3079 pour n = 6(k = 3).

6.3.3 Méthode de Newton

 6.4. RÉSOLUTION D’ÉQUATIONS DIFFÉRENTIELLES ORDINAIRES (ÉDO) 37

1 f u n c t i o n [ i n t ] = newtonc ( a , b ,m, fun )

2 h=(b−a ) /m;
3 n2=f i x (m/ 2 ) ;
4 i f m > 6 , e r r o r ( ’m vaut <= 6 ’ ) ; end
5 a03 =1/3; a08 =1/8; a45 =1/45; a288 =1/288; a140 =1/140;
6 a lp h a = [ 0 . 5 0 0 0 ; . . .
7 a03 4∗ a03 0 0 ; . . .
8 3∗ a08 9∗ a08 0 0 ; . . .
9 14∗ a45 64∗ a45 24∗ a45 0 ; . . .
10 95∗ a288 375∗ a288 250∗ a288 0 ; . . .
11 41∗ a140 216∗ a140 27∗ a140 272∗ a140 ] ;
12 x=a ; y ( 1 )=e v a l ( fun ) ;
13 f o r j =2:m+1
14 x=x+h ;
15 y ( j )=e v a l ( fun ) ;
16 end
17 i n t =0;
18 j = [ 1 : n2 + 1 ] ; i n t=sum ( y ( j ) . ∗ alpha (m, j ) ) ;
19 j =[n2 +2:m+ 1 ] ; i n t=i n t+sum ( y ( j ) . ∗ alpha (m,m−j +2) ) ;
20 i n t=i n t ∗h ;
21 return
Calculer l’intégrale : Z 8 √
I(f ) = x x2 + 1dx
0
1. Par la méthode de Newton (n = 3).
1 >> f = ’ x . ∗ s q r t (1+x . ^ 2 ) ’ ;
2 >> I =newtonc ( 0 , 8 , 3 , f )
3

4 I =
5 174.1024
R8 √
On trouve : I(f ) = 0 x x2 + 1dx ' 174, 1024.

6.4 Résolution d’équations différentielles ordinaires (ÉDO)

Soit l’équation différentielle du premier ordre :

y 0 (t)
= ty(t)
(6.1)
y (0) = 2

Voilà certainement l’un des exemples les plus simples que l’on puisse imaginer. la solution
particulière est :
2
y(t) = 2et /2 (6.2)
 38 CHAPITRE 6. TP MATLAB

qui vérifie à la fois l’équation différentielle et la condition initiale.

1 >> S o l = d s o l v e ( ’Dy = t ∗y ’ , ’ y ( 0 ) = 2 ’ )
2

3 Sol
4 2
5 2 exp ( 1 /2 t )
6 >> e z p l o t ( S o l , [ 0 2 ] )
7 >> l e g e n d ( ’ S o l u t i o n Exacte ’ )

2 /2
Figure 6.1 – y(t) = 2et avec y(0) = 2.

6.4.1 Méthode d’Euler

1 f u n c t i o n [ t , y ] = ER( f , t0 , T, N, y0 )
2 h = (T−t 0 ) /N ;
3 t = l i n s p a c e ( t0 , T,N+1) ;
4 y ( 1 ) = y0 ;
5 f o r n = 2 : N+1
6 y ( n ) =y ( n−1) + h∗ f e v a l ( f , t ( n−1) , y ( n−1) ) ;
7 end
8 end
On veut approcher la solution de (6.1) en t = 1 à l’aide de la méthode d’Euler, en
subdivisant l’intervalle [0, 1] en 10 parties égales.
1 >> f=i n l i n e ( ’ t ∗y ’ , ’ t ’ , ’ y ’ ) ;
2 >> [ tn , ytn ]=ER( f , 0 , 1 , 1 0 , 2 )
 6.4. RÉSOLUTION D’ÉQUATIONS DIFFÉRENTIELLES ORDINAIRES (ÉDO) 39

4 tn =
5 Columns 1 through 6
6

7 0 0.1000 0.2000 0.3000 0.4000

0.5000
8

9 Columns 7 through 11
10

11 0.6000 0.7000 0.8000 0.9000 1.0000

12

13 ytn =
14 Columns 1 through 6
15

16 2.0000 2.0000 2.0200 2.0604 2.1222

2.2071
17

18 Columns 7 through 11
19

20 2.3175 2.4565 2.6285 2.8387 3.0942

On trouver les valeurs :

tn 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

yn 2,0000 2,0000 2,0200 2,0604 2,1222 2,2071 2,3175 2,4565 2,6285 2,8387 3,0942

1 >> e z p l o t ( fun , [ 0 1 ] )
2 >> h o l d on
3 >> p l o t ( tn , ytn , ’ r−o ’ )
4 >> l e g e n d ( ’ Exacte ’ , ’ E u l e r ’ )
5 >> hold o f f
 40 CHAPITRE 6. TP MATLAB

Figure 6.2 – Méthode d’Euler : y 0 (t) = ty(t) avec y(0) = 2 et h = 0, 1.

6.4.2 Méthodes de Runge-Kutta d’ordre 2

1 f u n c t i o n [ t , y ] = RK2( f , t0 , T, N, y0 )
2 h = (T−t 0 ) /N ;
3 t = l i n s p a c e ( t0 , T,N+1) ;
4 y ( 1 ) = y0 ;
5 f o r n = 2 :N+1
6 k1 = f e v a l ( f , t ( n−1) , y ( n−1) ) ;
7 k2 = f e v a l ( f , t ( n−1)+0.5∗h , y ( n−1)+0.5∗h∗k1 ) ;
8 y ( n ) = y ( n−1) + h∗k2 ;
9 end
10 end
On veut approcher la solution de (6.1) en t = 1 à l’aide de la méthode de Runge-Kutta
d’ordre 2, en subdivisant l’intervalle [0, 1] en 10 parties égales.
1 >> f=i n l i n e ( ’ t ∗y ’ , ’ t ’ , ’ y ’ ) ;
2 >> [ tn , ytn ]=RK2( f , 0 , 1 , 1 0 , 2 )
3

4 tn =
5 Columns 1 through 6
6

7 0 0.1000 0.2000 0.3000 0.4000

0.5000
8
 6.4. RÉSOLUTION D’ÉQUATIONS DIFFÉRENTIELLES ORDINAIRES (ÉDO) 41

9 Columns 7 through 11
10

11 0.6000 0.7000 0.8000 0.9000 1.0000

12

13 ytn =
14 Columns 1 through 6
15

16 2.0000 2.0100 2.0403 2.0918 2.1661

2.2656
17

18 Columns 7 through 11
19

20 2.3933 2.5535 2.7517 2.9950 3.2923

On trouver les valeurs :
tn 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
yn 2,0000 2,0100 2,0403 2,0918 2,1661 2,2656 2,3933 2,5535 2,7517 2,9950 3,2923

1 >> e z p l o t ( fun , [ 0 1 ] )
2 >> h o l d on
3 >> p l o t ( tn , ytn , ’ r−o ’ )
4 >> l e g e n d ( ’ Exacte ’ , ’RK2 ’ )
5 >> hold o f f

Figure 6.3 – Méthodes de Runge-Kutta d’ordre 2 : y 0 (t) = ty(t) avec y(0) = 2 et h = 0, 1.