Ministère de l’Enseignement Supérieur, de la Recherche Scientifique et de la Technologie

Programmation Mathématique
Simulations sur SageMath/Python et Projet sur les RN

1ère année du cycle d’Ingénieur à l’Ecole Supérieure

de la Statistique et de l’Analyse de l’Information

Ines ABDELJAOUED TEJ

Année Universitaire 2017/2018

ESSAI, 6 rue des métiers - BP 675, Tunis 1080 ; Tél : 71 703 717 / 71 704 309 - Fax : 71 704 329
Mél : inestej@[Link] - Web : [Link]
Cours d’Analyse Numérique Ines Abdeljaoued Tej

Attention, malgré une re-lecture soigneuse, il est possible que ce cours contienne des fautes de
frappe ou des erreurs. Si vous en trouvez, merci de nous en avertir de préférence par e-mail :
inestej@[Link].

Avant-propos
Ce document regroupe les notes du cours enseigné en première année de l’Ecole Supérieure de
la Statistique et de l’Analyse de l’Information à l’Université de Carthage. Cet enseignement se
compose à la fois de cours magistraux et de séances de travaux pratiques.

Il s’agit de pésenter plusieurs méthodes numériques de base utilisées dans l’approximation des
fonctions par interpoation polynomiale. Nous introduisons également la résolution des équations
non linéaires et celle des systèmes linéaires. Le but de cette matière est de familiariser les élèves-
ingénieurs avec les techniques de simulations numériques, en abordant notamment la notion de
convergence, de précision et de stabilité. L’analyse des méthodes sont complétées par un travail
d’implémentation et d’application réalisé par les étudiants sur le logiciel SAGEmath/Python.

ESSAI - 2017/2018 1 inestej@[Link]

Table des matières

1 SAGEmath : une plateforme mathématique complète et libre 7

1.1 SAGEmath par l’exemple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.1.1 Commandes de base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.1.2 Exemples divers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.2 Application au codage/décodage RSA . . . . . . . . . . . . . . . . . . . . . . . . 13
1.2.1 Algorithme d’Al Kindi . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.2.2 Algorithme de Miller-Rabin . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.2.3 Algorithme RSA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.2.4 Correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2 Interpolation et Approximation 21
2.1 Formule d’Interpolation de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2 Algorithme d’Aitken . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.2.1 Formule de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.2.2 La méthode des différences divisées . . . . . . . . . . . . . . . . . . . . . 24
2.2.3 Algorithme de Neville . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.3 Complexité d’algorithmes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.4 Estimation de l’erreur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.5 Interpolation Spline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.6 Introduction à l’Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.6.1 Résultats fondamentaux . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.6.2 Approximation polynomiale discrète . . . . . . . . . . . . . . . . . . . . . 29
2.7 Exercices d’Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3 Intégration numérique 33
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.2 Formules de Newton-Côtes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.2.1 Méthode du trapèze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.2.2 Méthode de Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.2.3 Remarque sur la Formule Composite . . . . . . . . . . . . . . . . . . . . 35
3.3 Exemples d’Implémentations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

4 Racine d’équation non linéaire 36

4.1 Quelques algorithmes classiques . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.1.1 Méthode de la bissection . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.1.2 Méthode de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.1.3 Méthode de la sécante . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

2
Cours d’Analyse Numérique Ines Abdeljaoued Tej

4.1.4 Méthode du point fixe . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

4.2 Méthodes d’accélération de la convergence . . . . . . . . . . . . . . . . . . . . . 41
4.2.1 Procédé ∆2 d’Aitken . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.2.2 Méthode de Steffensen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.3 Exercices d’Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

5 Résolution de systèmes d’équations linéaires 43

5.1 Rappels et compléments sur les matrices . . . . . . . . . . . . . . . . . . . . . . 44
5.2 Normes matricielles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.2.1 Conditionnement d’une matrice . . . . . . . . . . . . . . . . . . . . . . . 45
5.2.2 La méthode de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.2.3 La factorisation LU d’une matrice . . . . . . . . . . . . . . . . . . . . . . 48
5.2.4 La factorisation de Cholesky . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.2.5 Valeurs propres et vecteurs propres . . . . . . . . . . . . . . . . . . . . . 50
5.2.6 La méthode de la puissance . . . . . . . . . . . . . . . . . . . . . . . . . 50

6 Equations différentielles ordinaires 54

6.1 Rappel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
6.2 Problème de Cauchy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
6.3 Méthode d’Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
6.4 Méthode de Taylor d’ordre 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
6.5 Méthode de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

7 Application aux réseaux de neurones 57

7.1 Historique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
7.2 Définitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
7.2.1 Neurone formel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
7.2.2 Fonctions d’activation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
7.2.3 Perceptron multi-couches . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
7.3 Fonctionnement du réseau de neurones . . . . . . . . . . . . . . . . . . . . . . . 59
7.4 Statistique et PMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
7.5 Logiciels utilisés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
7.5.1 nnet de R pour le perceptron multi-couches . . . . . . . . . . . . . . . . 64
7.5.2 SageMath, Python . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
7.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

A Travaux Pratiques sur SAGEmath 66

A.1 Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
A.2 Résolution d’équations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
A.3 Méthode d’Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
A.4 Analyse Numérique Matricielle : SVD . . . . . . . . . . . . . . . . . . . . . . . . 70
A.4.1 A faire sur Python : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
A.4.2 Correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
A.4.3 Un code plus intéracif . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
A.5 Matrice bien ou mal conditionnée . . . . . . . . . . . . . . . . . . . . . . . . . . 74
A.6 Factorisation LU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
A.7 Factorisation de Cholesky . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
A.8 Méthodes itératives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

ESSAI - 2017/2018 3 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

B Exerices de Révision et Applications 78

C Tests, DS et Examens corrigés 82

C.1 2008-2009 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
C.2 2009-2010 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
C.3 2010-2011 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

ESSAI - 2017/2018 4 inestej@[Link]

Introduction

Scientifiques et ingénieurs approchent les mathématiques de haut niveau de deux manières - en

utilisant le Calcul Numérique ainsi que le Calcul Formel. L’analyse numérique qui fait partie
du Calcul Scientifique est basée sur les ordinateurs.

L’avantage de l’analyse numérique est de pouvoir donner une réponse numérique à un pro-
blème qui n’a pas de solution analytique. Par conséquent, les résultats ne sont qu’approximatifs
contrairement aux résultats du Calcul Formel. Un cours d’analyse numérique va vous apprendre
à écrire les problèmes mathématiques en code, ainsi que vous apprendre quelles sortes de mo-
dèles mathématiques fonctionnent le mieux avec différentes sortes d’équations.

La démarche du numéricien commence par la modélisation du problème ou sa mise en équation.

Une analyse mathématique permet ensuite d’étudier le problème, de déterminer sous quelles
conditions il existe une solution unique et sinon, choisir celle qui correspond au mieux à la réalité.

Lorsque l’ingénieur ne peut pas calculer toutes les valeurs, le mieux est d’approcher le problème
en calculant des valeurs de certains points. C’est la discrétisation, qui doit être accompagnée
d’une mesure de l’erreur commise, c’est-à-dire mesurer la distance entre la solution approchée et
la véritable solution du problème. Une étape de simplification, comme la linéarisation pour des
problèmes non linéaires, permet d’obtenir des problèmes pour lesquels on dispose de méthodes
de résolution numérique performantes : résolution de systèmes linéaires, optimisation par une
méthode de gradient, recherche de valeurs propres.

Comme suite aux résultats du calcul par ordinateur, il convient d’être circonspect et d’étudier
les sources d’erreurs : les données elles-mêmes, l’arrondi, la méthode utilisée. La démarche qui
vient d’être présentée est celle que suit le numéricien en charge de la simulation numérique d’un
problème.

Le cours est structuré en trois grands chapitres :

— Interpolation et Approximation
— Racine d’équations non linéaires
— Résolution de systèmes d’équations linéaires.

Etant donné la valeur connue d’une certaine fonction en un certain nombre de points, quelle
valeur prend cette fonction en un autre point quelconque situé entre deux points donnés ? Une
méthode très simple est d’utiliser l’interpolation linéaire, qui suppose que la fonction inconnue
évolue linéairement entre chaque paire de points successifs connus. Nous introduisons dans le

5
Cours d’Analyse Numérique Ines Abdeljaoued Tej

premier chapitre l’interpolation polynomiale ainsi que d’autres méthodes d’interpolation utili-
sant des fonctions localisées telles que les splines.

Un autre problème fondamental est le calcul des solutions d’une équation donnée. Les algo-
rithmes de recherche de racines d’une fonction sont utilisés pour rśoudre les équations non
linéaires. Si la fonction est différentiable et que sa dérivée est connue, alors nous choisirons la
méthode de Newton sinon, d’autres méthodes seront détaillés dans le deuxième chapitre.

Le troisième chapitre porte sur l’étude des méthodes de calcul du vecteur x solution du système
d’équations linéaires Ax = b. La première partie de ce chapitre sera consacrée aux méthodes
directes et la deuxième aux méthodes itératives. Nous aborderons ensuite en troisième partie
les méthodes de résolution de problèmes aux valeurs propres.

Le champ de l’analyse numérique est divisé en différentes disciplines suivant le type de pro-
blème à résoudre, et chaque discipline étudie diverses méthodes de résolution des problèmes
correspondants. Il s’agit dans un premier temps d’étudier les principales méthodes en Analyse
Numérique (Méthode de Lagrange, algorithmes du point fixe, de Newton, du gradient conjugué,
la factorisation LU, etc...). Nous nous intéressons ensuite aux principes qui sous-tendent ces
méthodes ainsi qu’à leur mise en oeuvre pratique en langage Scilab et Maxima. Il s’agira enfin
d’implémenter et de tester plusieurs algorithmes d’Analyse Numérique et plusieurs approches
de résolution.

Nous donnons à la fin de chaque chapitre une liste d’exercices de difficultés variables. Ces exer-
cices ne sont profitables que si l’élève-ingénieur les travaille. Qu’il garde à l’esprit ce proverbe
chinois :

J’entends et j’oublie, (cours oral)

je vois et je retiens, (étude du cours)
je fais et je comprends (exercices).

Pour conclure, ce cours est une introduction à des approches plus complexes des mathématiques
de l’ingénieur : d’autres thématiques seront développés dans les cours d’Optimisation Convexe,
de Calcul de Variations, de Recherche Opérationnelle, etc...

Enfin, rappelez-vous la devise suivante :

That which is learned without desire is soon forgotten.

ESSAI - 2017/2018 6 inestej@[Link]

Chapitre 1

SAGEmath : une plateforme mathématique

complète et libre

Python est un langage de programmation adapté au calcul scientifique. Il est interprété, inter-
actif et orienté objet. Il possède de nombreuses extensions pour le calcul, le traçage, le stockage
de données. Combiné avec des fonctionnalités comme Tkinter, une interface graphique pour les
utilisateurs, Python permet de développer de bonnes interfaces graphiques pour les codes. Par
exemple, les principaux mots cles sous Python sont :

and else import assert except in break

exec is class finally lambda continue for
not def from or del global pass elif if print
raise return try while yield

L’aspect le plus intéressant est qu’il existe un ensemble de modules pour le calcul scientifique
sous Python. Des modules pour pratiquement tout (vecteurs, tenseurs, transformations, déri-
vées, algèbre linéaire, transformées de Fourier, statistiques, groupes, etc) sont disponibles. Il
est aussi possible d’utiliser des bibliothèques C et Fortran à partir de Python. Finalement,
écrire un programme numérique est chose aisée en Python. Il est également possible de réaliser
des algorithmes et des calculs parallèles : des interfaces existent vers netCDF (fichiers binaires
portables), MPI et BSPlib.
Comme exemple, existe Scipy, la bibliothèque d’outils scientifiques libres pour Python. Elle
inclut des modules pour créer des graphiques et faire du traçage, des modules d’optimisation,
d’intégration, de fonctions spéciales, de manipulation de signaux et d’images, d’algorithmes
génétiques, de solveurs d’équations différentielles ordinaires, etc.
Le responsable, K. Hinsen, propose un bon didacticiel intitulé Scientific Computing in Python
(Calcul scientifique avec Python). Consuler aussi I. Lima dans [?].

1.1 SAGEmath par l’exemple

La plateforme SAGEmath est un ensemble de logiciels libres et payants (sous réserve d’avoir la
licence) qui cohabitent ensemble grâce au langage Python.

7
Cours d’Analyse Numérique Ines Abdeljaoued Tej

1.1.1 Commandes de base

Afin de mieux illustrer les utilisations de base de SAGEmath, nous nous sommes largement
inspirés de [?]. Mais avant tout, télécharger [Link] depuis [Link] (280MB) et
taper les commandes suivantes :
$ tar xf [Link]
$ cd sage-4.3.1
$ make
Aller prendre un café (ou plusieurs). . .
1. Les opérations arithmétiques sont données par :
2+3
15 + 3 * 2 / 8 - 2^3
type(21)
--2
30/3
n(30/3, 4)
2^3
8%5
boite=7
7
2. Les chaînes de caractères :
a="Bonjour, je suis une chaine de caractere"
print(a)
type(a)
<type ’str’>.
3. Opérateurs conditionnels :
x = 2
y = 3
print x == y
print x <> y
print x, "!=", y, ":", x != y
print x < y
print x > y
print x >= y
False.
4. Utiliser ce code avec la variable s et obtenir la solution :
var(’a,b,c’)
eqn = [a+b*c==1, b-a*c==0, a+b==5]
s = solve(eqn, a,b,c)
La solution de eqn = [bc + a = 1, −ac + b = 0, a + b = 5] est :
" √ √ #
25i 79 + 25 5i 79 + 5 1 √ 1
a= √ ,b = √ , c = i 79 +
6i 79 − 34 i 79 + 11 10 10

ESSAI - 2017/2018 8 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

" √ √ #
25i 79 − 25 5i 79 − 5 1 √ 1
a= √ ,b = √ , c = − i 79 +
6i 79 + 34 i 79 − 11 10 10

5. var(’x’)
f=-x^3+3*x^2+7*x^2-4
On construit l’objet courbe, stocké dans la variable p :
p = plot(f, x, -5, 5)
350

300

250

200

150

100

50

-4 -2 2 4

6. La commande if :
x = 3
print x > 5
if x > 5:
print x
print "plus grand"
print "Fin du programme"
else:
print "Introduire une autre valeur"
7. L’opérateur booléen and :
a = 7
b = 9
print a > 5 and b > 10
if a > 5 and b < 10:
print "Les deux expressions sont vraies"
8. L’opérateur booléen or :
print a > 5 or b > 10
if a > 5 or b > 10:
print "Au moins une des deux expressions est vraie"
9. L’opérateur booléen not :
a = 7
print a > 5
print not a > 5
10. La boucle while :
x = 1 #initialisation
while x <= 40:
print x,
x = x + 2 #incrementation de x par 2
11. Les ensembles :

ESSAI - 2017/2018 9 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

x = Set([2,3,51,21])
print type(x)
print x[0]
print x
y = Set([51])
[Link](y)
a = Set([1,3,2,5,5])
[Link]()
3 in a
[Link](x)
[Link](x)
12. Les listes :
y = [1,1/2,0.75,’Bonjour’,x]
type(y)
y[3]
y[4]
add([1,2,3,4])
13. Les n-uplets :
z = (1,3,5,4)
type(z)
z[0]
a,b,c = (1,2,3)
a; b; c
14. Les opérateurs in et not in :
54 in [2,54,2]
30 not in [2,54,2]
15. La boucle for :
for i in [1,2,3,4,5]:
i,
16. Les fonctions :
def salut(s):
"""
Retourne Bonjour ou Au revoir.
"""
if s>3:
return "Bonjour"
else:
return "Au revoir"

html("<H1>Exemple de fonctions interactive</H1>")

html("<I> Juste pour voir</I>")
@interact
def _(p=1):
salut(p)

ESSAI - 2017/2018 10 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

17. Quelle est l’utilité des commandes suivantes :

denominator(12/23)
divisors(20)
gcd(40,116)
mul([2,3,4])
len([1,2,3,’hello’])
srange(11)
a=srange(4,step=1.3)
b=srange(-5,-1)
zip(a,b)
[(0.000000000000000, −5) , (1.30000000000000, −4) , (2.60000000000000, −3) , (3.90000000000000, −2)] .
18. La boucle for :
for t in srange(6):
t,

d=zip(a,b)
for (x,y) in d:
x,y

[2*t for t in [0,1,2,3]]

[t^2 for t in range(20) if t%2==0]

19. Obtenir des résultats numériques :
a = pi
a
n(a)
a.n(digits=30)
a.n(prec=12)
20. Calcul symbolique :
var(’a’)
type(2*a)
m = 5 + a
n = 8 + a
y = m * n
z = [Link]()
y(5)
130.
21. Equations symboliques : [a = (−5)] .
22. Divers :
A=matrix(3,3,[1,2,3,2,3,4,1,2,3])
latex(A)
23. Courbes en 2 ou 3 dimensions :

ESSAI - 2017/2018 11 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

24. Calcul de 0 + 1 + · · · + N − 1 :
def sum1(N):
s = 0
for k in range(N):
s += k
return s
SAGEmath: time sum1(10^6)
Cette fonction peut être sauvegardée dans un fichier ([Link] par exemple). Les com-
mandes suivantes permettent de lire ou d’exécuter le fichier :
SAGEmath: cat [Link]
SAGEmath: load [Link]
SAGEmath: time sum2(10^6)

1.1.2 Exemples divers

1. On considère les vecteurs et les matrices suivants
     
1 1 2 0 0 0 1 2 3 4 5



 2 


 0 0 3 2 1 


 2 3 4 5 6 

l= 1 2 3 4 5 ,v = 
 3 ,A = 
  0 0 0 1 1 ,B = 
  3 4 5 6 7 

 4   0 0 0 0 2   4 5 6 7 8 
5 1 1 1 0 0 5 6 7 8 9

Créer et commenter le résultat des instructions l ∗ v, A + B, A ∗ B, B/A, t v.

2. On appelle méthode de Monte-Carlo toute méthode visant à calculer une approximation
numérique par utilisation de procédés aléatoires. Nous allons calculer une valeur appro-
chée de π par la méthode de Monte-Carlo 1 .
On tire au hasard des points de coordonnées (x, y) dans un carré de longueur égale à 2.
On vérifie s’ils appartiennent ou non à un disque de rayon 1 et de centre, le centre du
carré. Notons que ces points peuvent être tirés avec la même probabilité dans l’ensemble
du carré. D’après la loi des grands nombre, lorsque le nombre de tirage tend vers l’infini,
le rapport entre le nombre de points tirés dans le disque et le nombre de points tirés au
total tend vers le rapport des surfaces du cercle et du carré, soit π4 .
3. Etant donné n > 0, écrire en utilisant le moins de lignes possible, un script SAGEmath
1
qui déclare la matrice de Cauchy de coefficients ai,j = i+j .
4. On peut toujours transformer un problème du type f (x) = 0 en un problème de la forme
g(x) = x et ce d’une infinité de façons. Le théorème suivant caractérise les fonctions g
qui donnent des suites convergentes : Si dans l’intervalle [a, b], g vérifie les conditions
suivantes :
(a) x ∈ [a, b] alors g(x) ∈ [a, b],
0
(b) l’application g est strictement contractante, c’est-à-dire maxx∈[a,b] | g (x) |< 1 .
Alors pour tout x0 ∈ [a, b], la suite définie par xn+1 = g(xn ) converge vers l’unique point
fixe x∗ de g sur [a, b].

1. N. Metropolis et S. Ulam ont inventé en 1947 ce type de méthode. Le nom fait allusion aux jeux de hasard
pratiqués dans les casinos de Monte-Carlo

ESSAI - 2017/2018 12 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Ecrire une fonction prenant en entrée g, x0 ,  et un nombre d’itérations maximal N et

qui rend une approximation du point fixe de g à  près.

1.2 Application au codage/décodage RSA

Ce TP est utilisé pour une initialisation à Python. C’est un bon exemple qui permet la manipu-
lation des grands nombres pseudo-aléatoires, il permet de voir la puissance d’un logiciel comme
Python sur des problèmes de décryptage de données. On utilise principalement les résultats
d’arithmétique (petit théorème de Fermat) pour répondre à la problématique de chiffrement.
Nous encourageons les étudiants à suivre les étapes de ce TP pour se familiariser avec Python.
Le système de chiffrement RSA est une méthode qui repose sur certains concepts de la théorie
des nombres. C’est une méthode à clé publique, nommée à partir des noms de famille de ses
inventeurs Ronald Rivest, Adi Shamir et Leonard Adleman. Ce TP a pour but d’étudier le
principe de fonctionnement de l’algorithme RSA. Pour cela, nous allons d’abord étudier un
code permettant de décrypter et de décrypter un message. A la suite, nous allons voir comment
générer des clés publiques et privées pour crypter un message.
Mais avant cela, la première section présentera un petit exemple ludique de décryptage d’un
message, permettant d’introduire un certain nombre de vocabulaire lié à la cryptographie. Cet
exemple comprend une brève explication du cryptage de Jules César ainsi que de l’algorithme
de décriptage d’Al Kindi.

1.2.1 Algorithme d’Al Kindi

Pour décoder un message texte, on va utiliser l’algorithme de Abu Yusuf Yaqub Ibn Ishaq Al
Kindi qui a vécu à Bagdad de 801 à 873. Cet algorithme consiste à comparer la fréquence
d’apparition des lettres codées avec la fréquence des lettres dans la langue française. En effet,
une lettre qui apparaît souvent dans la langue française, va apparaître souvent (via son code)
dans un message codé. Ci-dessous les fréquences des lettres dans la langue française :

Lettres E S A I T N R U L O D C P
Fréquences (%) 14.7 7.9 7.6 7.5 7.2 7.1 6.6 6.3 5.5 5.4 3.7 3.3 3
Soit le message codé suivant 2 :
code = "YR PNEER QR Y ULCBGRAHFR "
code += "RFG RTNY FV WR AR Z NOHFR "
code += "N YN FBZZR QRF PNEERF "
code += "QRF QRHK NHGERF PBGRF"
On crée l’objet message composé des lettres de code et l’ensemble lettres des différentes
lettres apparaissant dans le message codé :
message = [i for i in code]
lettres = set(message)
1. Si une lettre apparaît souvent dans un message, on s’attend à ce qu’elle apparaisse autant
dans un message codé. On calcule donc le tableau d’effectifs en comptant les lettres du
message :
2. TP inspiré des supports du Laboratoire d’informatique et de mathématiques de l’Institut de recherche sur
l’enseignement des mathématiques de la Réunion : http ://[Link]/IMG/pdf/[Link]

ESSAI - 2017/2018 13 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

for j in lettres:
print(j+’-->’+str([Link](j)))
Comparer ce résultat avec le tableau ci-dessus, et donnez le décodage de la lettre E, S
et A.
2. Pour continuer le décodage, et vu que les fréquences d’apparition des lettres sont trop
proches, on peut utiliser les propriétés lingustiques, comme le fait qu’un mot long com-
mence souvent par une consonne ou qu’il y a peu de mots de deux lettres. On peut
également tester si le procédé de chiffrement utilisé est bien celui proposé par Jules Cé-
sar.

Ce procédé date de 120 avant J.C., il utilise la technique suivante : il suffit de procéder
à une permutation circulaire de toutes les lettres de l’alphabet, en remplaçant chaque
lettre par celle qui est située à s rangs plus loin (ici s correspond à la clé de chiffrement).
Dans cet exemple, la lettre A correspond à N, la lettre E correspond à R et la lettre F
correspond à S. Vérifiez, grâce au code suivant, que la clé de chiffrement est bien s = 13
et en déduire le message déchiffré.
s = 13 # s est la cle de chiffrement selon la méthode de Jules César
n=len(message)
MessageDechiffre=""
for i in range(n):
caractere = message[i]
p = ord(caractere) #Varie entre 65 et 91
q = p + s
if q > 91:
q = q - 26
if q == 91:
q = 65
r = chr(q)
print caractere, p, r, q
MessageDechiffre = MessageDechiffre + r
print MessageDechiffre
3. Pourquoi est-ce que le codage de César est facile à casser ? Indication : trouver le nombre
de clés différentes qu’on peut utiliser pour le codage de César.

1.2.2 Algorithme de Miller-Rabin

Le but de cette section est de construire des nombres premiers, le plus grand possible. Pour
cela, on va utiliser l’algorithme de Miller-Rabin qui permet de tester si un nombre donné est
premier ou pas, en un temps raisonnable.

L’algorithme se base sur le théorème suivant : Théorème. Soit p un nombre premier. On pose
s la puissance maximale de 2 divisant p − 1 :

p − 1 = 2s .d avec d pair.

ESSAI - 2017/2018 14 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Pour tout 1 ≤ n ≤ p − 1, on a l’une des possibilités suivantes :

jd
nd = 1 mod p ou alors n2 = −1 mod p pour un certain 0 ≤ j < s.
Il s’agit donc de choisir au hasard un nombre n jusqu’à ce que l’une des deux conditions précé-
dente n’est pas satisfaite. Cet n est appelé "témoin" du fait que p est un nombre composé. Si
p n’est pas premier, alors il existe 3/4 de chance que n soit "témoin". L’idée est donc d’utiliser
le petit théorème de Fermat i.e. le fait que ap−1 = 1 mod p.

1. On désigne par t la taille en bits de la clé RSA que l’on souhaite utiliser. Il nous faut donc
trouver deux nombres premiers de taille t/2. Pour cela, on utilise la fonction random pour
trouver un nombre de t/2 bits au hasard puis netx_prime le nombre premier suivant.
2. Déterminez le plus petit nombre premier de 10 chiffres (pensez à la conversion en écriture
binaire et à la commande random_prime(2ˆ1023, 2ˆ1024)).
import random
randint(10\^9, 10\^10)
[Link](2048)
3. Implémentez l’algorithme de Miller-Rabin. Le pseudo-code obtenu de Wikipédia est
donné ci-dessous. La fonction Témoin_de_Miller(a,n) renvoie Vrai si a est un témoin
de Miller que n est composé, Faux si n est fortement pseudo-premier en base a. Elle
prend en entrée un nombre n entier, impair ≥ 3 et a un entier > 1. On calcule s et d
tels que n − 1 = 2s d avec d impair. Ici s > 0 car n est impair. On note x le reste de la
division de ad par n, puis dans la boucle Tant que, le reste de la division de x2 par n.
Témoin_de_Miller(a, n):
x := a^d % n #x entier reste de la division de ad par n
si x = 1 ou x = n - 1
renvoyer(Faux)
Tant que s > 1
x := x^2 % n
si x = n - 1
renvoyer(Faux)
s := s - 1
Fin de boucle tant que
renvoyer(Vrai)
Le test de Miller-Rabin peut alors être décrit comme suit, Miller-Rabin(n, k) renvoie
Vrai si n est fortement pseudo-premier en base a pour k entiers a, Faux s’il est composé.
C’est un algorithme qui a une complexité polynomiale en terme de temps de calcul :
O([Link](n)3 ) obtenue en remarquant que la décomposition n − 1 = 2s d se calcule en
O(log(n)).
Miller-Rabin(n,k):
répéter k fois :
choisir a aléatoirement dans l’intervalle [2, n-2]
si Témoin_de_Miller(a,n)
renvoyer(Faux)
Fin de boucle répéter
renvoyer(Vrai)

ESSAI - 2017/2018 15 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

4. La fonction is_prime est un test de primalité qui prouve si un nombre donné est premier
ou pas. La fonction is_pseudoprime est un test probabiliste de primalité : l’option > 0
utilise le test de primalité probabiliste de Miller-Rabin. Ce test s’appuie sur le petit
théorème de Fermat. Testez ces deux fonctions et expliquez à quoi servent les fonctions
netx_prime et random_prime ou encore next_probable_prime 3 .

1.2.3 Algorithme RSA

Le principe du système de chiffrement RSA est le suivant 4 :

Trouver deux nombres premiers p et q p et q doivent rester secrets

Calculer n = pq n est clé publique
Choisir e tel que son PGCD avec (p − 1)(q − 1) soit 1 e est clé publique
Trouver d tel que ed ≡ 1 (modulo (p − 1)(q − 1)) d est la clé privée.
Le cryptage fonctionne de la manière suivante : il y a deux personnes Ali et Besma qui sou-
haitent communiquer en toute sécurité. Soit M un nombre qui correspond au message à envoyer.
B génère n ainsi que les exposants de déchiffrements e et d. B envoie à A les nombres n et e
qui constituent la clé publique 5 .

Ici Ali est l’émetteur et il ne connaît que la clé publique, c’est-à-dire les nombres e et n. Il
calcule M e ≡ N (modulo n) et envoie au destinataire le message N . Le récepteur qui est Besma
calcule N d ≡ W (modulo n). On constate que W est exactement le message M .

Exercice à faire : Les messages qu’on peut chiffrer sont les entiers M compris entre 0 et n-1.
Le codage du message M s’obtient par le calcul de N et le déchiffrement par le calcul de W .
1. Calculez N = M e (modulo n). Pour cela utilisez la commande help(mod) ou encore la
commande %.
2. Le décodage ou le déchiffrement d’un message codé s’obtient par le calcul de W = Nd
modulo n. Donnez deux fonctions de cryptage et de décryptage, par clés RSA : on
appellera les fonctions N = encrypt(e, n, M) et W = decrypt(d, n, N).
3. Les attaques actuelles du RSA se font essentiellement en factorisant l’entier n de la
clé publique. La sécurité du RSA repose donc sur la difficulté à factoriser de grands
entiers. Pour garantir une certaine fiabilité des transactions, il faut choisir des clés
plus grandes : des clés de 21024 voire 22048 bits, i.e. des clés à 1024*ln(2)/ln(10) ou
2048*ln(2)/ln(10) chiffres.
La comande factor factorise les entiers qu’on lui passe en paramètre. Générez des clés
RSA de plus en plus grande et constatez l’augmentation du temps de calcul de la facto-
risation de ces données.
Le RSA est une méthode de chiffrement assez sûre, mais il n’est pas exclu que des services
secrets comme la NSA ou le FSB aient réussi le décryptage du RSA.
3. http : //[Link]/les − nombres − en − python/, voir aussi
http ://[Link]/staff/[Link]/[Link] pour une bonne utilisation des outils de
génération de nombre pseudo-aléatoires
4. Le nombre e est appelé la graine et on pose φ = (p − 1)(q − 1).
5. Le code python pour l’algorithme RSA a été inspiré de https ://[Link]/JonCooperWorks/5314103
alors que la structure du TP a suivi les travaux pratiques du master Informatique de l’Université de Lille 1 :
https ://[Link]/prive/sgambs/TP7_intro_securite.pdf

ESSAI - 2017/2018 16 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

1.2.4 Correction
Algorithme d’Al Kindi
Soit le message chiffré suivant :

YR PNEER QR Y ULCBGRAHFR RFG RTNY FV WR AR Z NOHFR N YN FBZZR QRF PNEERF QRF

QRHK NHGERF PBGRF

Il y a 20 lettres distinctes : qui apparaissent avec une fréquence détaillée ci-dessous :

Lettres A C B E G F H K L O N Q P R U T W V Y Z
Fréquences (%) 2 1 3 5 4 9 4 1 1 1 7 4 3 17 1 1 1 1 4 3

Les trois premières fréquences indiquent les lettres E, S et N. Ainsi, la lettre R correspond très
probablement à la lettre E, la lettre F correspond à la lettre S et la lettre N correspond à la
lettre A. En tenant compte de la clé de chiffrement s = 13, on peut en déduire le message
déchiffré :

LE CARRE DE L HYPOTENUSE EST EGAL SI JE NE M ABUSE A LA SOMME DES CARRES DES

DEUX AUTRES COTES

Les procédures sur Python sont les suivantes. On peut s’amuser à coder et à décoder, grâce à
la méthode de déchiffrement de Cesar, n’importe quel texte contenant les lettres majuscules de
l’alphabet.

def codage_cesar(cle, message):

#message est sous format string
n = len(message)
MessageChiffre = ""
for i in range(n):
caractere = message[i]
p = ord(caractere) #Varie entre 65 et 91
q = p + s
if q > 91:
q = q - 26
if q == 91:
q = 65
r = chr(q)
#print caractere, p, r, q
MessageChiffre=MessageChiffre+r
return(MessageChiffre)

message = "LE CARRE DE L HYPOTENUSE"

codage_cesar(cle, message)

def decodage_cesar(cle,code):

ESSAI - 2017/2018 17 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

n = len(code)
MessageDechiffre = ""
for i in range(n):
caractere = code[i]
p = ord(caractere) #Varie entre 65 et 91
q = p + s
if q > 91:
q = q - 26
if q == 91:
q = 65
r = chr(q)
#print caractere, p, r, q
MessageDechiffre=MessageDechiffre+r
return(MessageDechiffre)

decodage_cesar(cle,code)

Algorithme de Miller-Rabin
sage: import random
sage: def decompose(n):
... s = 0
... d = n
... while d % 2 == 0:
... d = d/2
... s += 1
... return s, d
sage: decompose(33)
(0, 33)
sage: def Temoin_de_Miller(a,n):
... s,d = decompose(n-1)
... #print n-1, 2^s*d
... x = a^d % n
... if x == 1 or x == n - 1:
... return(False)
... while s > 1:
... x = x^2 % n
... if x == n - 1:
... return(False)
... s = s - 1
... return(True)
sage: def Miller_Rabin(n,k):
... i = 0
... while i<k:
... a = randint(2,n-2)
... if Temoin_de_Miller(a,n)==True:
... return(False)

ESSAI - 2017/2018 18 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

... i += 1
... return(True)
sage: n = 122201
sage: a = 10
sage: k = 10
sage: print Temoin_de_Miller(a, n), Miller_Rabin(n, k)
False True
sage: next_prime(122187)
122201

Algorithme RSA
sage: import random
sage: """
sage: Ce que fait Besma :
... D’abord B génère n ainsi que e et d.
... Ensuite B envoie les clés publiques n et e à Ali.
...
sage: """
sage: p = 100
sage: p = next_prime(p)
sage: q = next_prime(p)
sage: print p,q
sage: n = p*q
sage: #la graine est choisie par hasard ;
sage: #il faut vérifier que e et phi sont premiers entre-eux :
sage: phi = (p-1)*(q-1)
sage: e = [Link](1,phi)
sage: print e, n, "les exposants de déchiffrements"
101 103
2522 10403 les exposants de déchiffrements
sage: phi = (p-1)*(q-1)
sage: #On cherche d tel que d*e+a*phi=1
sage: #print phi, e, n, d
sage: """
sage: Ali reçoit les clés de déchiffrement n et e.
sage: Pour envoyer le message M codé, A calcule N :
sage: """
sage: # Soit un message M :
sage: M = randint(2,n-1)
sage: print M
sage: def encrypt(e, n, M):
... return pow(M,e,n)
sage: N = encrypt(e, n, M)
sage: print N
sage: """
sage: Ce que fait B : En recevant N de Ali, Besma calcule W :

ESSAI - 2017/2018 19 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

... B calcule d qui est la clé privée

... Ensuite, elle détermine W qui est le M déchiffré.
sage: """
sage: def decrypt(d, n, N):
... return(pow(N, d, n))
sage: g, a, d = xgcd(phi, e) #phi*a+d*e=g
sage: while g<>1 and d<=0:
... g, a, d = xgcd(phi, e)
... e = [Link](1,phi)
...
...
sage: W = decrypt(d, n, N)
sage: print p, q, n, e, d, M, N, W, M==W
101 103 10403 9019 19 3805 9308 4727 False
sage: """
sage: Pour casser le codage, il suffit de factoriser n et de récupérer p et q
sage: """
sage: factor(n)
101 * 103

ESSAI - 2017/2018 20 inestej@[Link]

Chapitre 2

Interpolation et Approximation

Soient n + 1 couples : (x0 , f0 ), (x1 , f1 ), . . ., (xn , fn ) ∈ R × R tels que les xi soient tous distincts.
On cherche un polynôme P : R −→ R tel que P (xi ) = fi pour i = 0..n.

Nous introduisons dans un premier temps la formule d’interpolation de Lagrange qui prouve
l’existence et l’unicité du polynôme P de degré au plus n. Ensuite nous présenterons l’algorithme
d’Aitken et la méthode des différence divisées qui utilise la formule d’interpolation de Newton.
Puis, nous introduirons l’algorithme de Neville qui présente une méthode très pratique pour le
calcul sur machine du polynôme P en un point. Nous introduirons enfin l’interpolation Spline
qui est un outil fondamental du graphisme sur ordinateur.

2.1 Formule d’Interpolation de Lagrange

Généralement les réels fi représentent le résultat d’expériences et sont l’expression d’une fonc-
tion f continue dans R en les points xi . La formule d’interpolation de Lagrange est un moyen
qui permet d’approcher la fonction f .

Théorème 2.1.1. Soient (xi , fi ) pour i = 0..n des points réels ou complexes, avec xi 6= xj pour
i 6= j. Il existe un unique polynôme P de degré inférieur ou égal à n vérifiant P (xi ) = fi pour
i = 0..n.

La preuve du théorème 2.1.1 se base sur les polynômes d’interpolation de Lagrange

Yj=n x − xj
Li (x) =
j=0,j6=i xi − xj

qui
Pi=nvérifient Li (xi ) = 1 et pour j 6= i : Li (xj ) = 0. Le polynôme P est alors égal à P (x) =
i=0 fi Li (x) et il est appelé polynôme ou formule de Lagrange.

Démonstration : Soient P et Q deux polynômes de degré au plus égal à n et vérifiant P (xi ) =

Q(xi ) = fi pour i = 0..n. Alors le polynôme r = P − Q (qui est de degré au plus égal à n)
s’annule en x0 , x1 , . . . , xn . C’est donc un polynôme qui a n + 1 racines : il est nécessairement
nul. Donc P = Q, d’où l’unicité.

21
Cours d’Analyse Numérique Ines Abdeljaoued Tej

Pour montrer l’existence d’un tel polynôme, il suffit de remarquer que la famille (L0 , L1 , . . . , Ln )
forme une base de l’espace vectoriel des polynômes
Pi=n de degré inférieur ou égal à n (qui est de
dimension n + 1). En effet, pour tout x, i=0 αi Li (x) = 0 implique en particulier que αk = 0
pour k = 0..n (onPremplace x par xk ). On retrouve la formule de Lagrange en écrivant P dans
cette base : P = i=ni=0 fi Li (x) et P (xi ) = fi pour tout i = ..n.

Exemple 2.1.1. Le polynôme P de degré inférieur ou égal à 3 vérifiant

P (0) = 1, P (1) = 2, P (−1) = −1 et P (2) = 4
est P (x) = 1.L0 (x) + 2.L1 (x) + (−1).L−1 (x) + 4.L2 (x) = 13 x3 − 21 x2 + 67 x + 1 où :
(x−1)(x+1)(x−2) (x−0)(x+1)(x−2)
L0 (x) = (0−1)(0+1)(0−2)
, L1 (x) = (1−0)(1+1)(1−2)
,
(x−0)(x−1)(x−2) (x−0)(x−1)(x+1)
L−1 (x) = (−1−0)(−1−1)(−1−2)
, L2 (x) = (2−0)(2−1)(2+1)
.

2
Exercice 2.1.1. Construire les polynômes d’interpolation de Lagrange de f (x) = e−x sur les
entiers de [−2, 2].
Exercice 2.1.2. Ecrire un algorithme qui calcule Li (x) pour n, i et x fixés.

Remarque 2.1.1. Si on utilise la formule de Lagrange, le nombre d’opérations nécessaires pour

calculer P (x) est de l’ordre de 4n2 opérations pour chaque valeur de x.

Sur ordinateur, on préfère utiliser la formule barycentrique :

Pk=n ak fk
k=0 x−xk
P (x) = Pk=n ak
k=0 x−xk

avec ak = v0 (x1 k ) et v(x) = (x−x0 )(x−x1 ) . . . (x−xn ). Cette formule nécessite pour une première
valeur de x l’ordre de 2n2 opérations afin de déterminer P (x). Pour d’autres valeurs de x, 2n
opérations sont nécessaire au calcul de P a condition de conserver les valeurs ak .

2.2 Algorithme d’Aitken

Soit la relation de récurrence suivante :
(
Fk,0 (x) = fk , k = 0..n
Fj,j (x)(xk −x)−Fk,j (x)(xj −x)
Fk,j+1 (x) = xk −xj
, k = j + 1..n

Théorème 2.2.1. Le polynôme de Lagrange P (x) sur les points x0 , x1 , . . . , xn est égal à Fn,n (x).

La preuve du théorème 2.2.1 se base sur le fait que Fk,j (x) avec k ≥ j est le polynôme de
Lagrange sur les points x0 , x1 , . . . , xj−1 , xk .

Algorithme. On calcule P (x) = Fn,n (x) à l’aide du tableau triangulaire suivant :

ESSAI - 2017/2018 22 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

k/j 0 1 2 ... j−1 j ... n−1 xk − x

0 F0,0 x0 − x
1 F1,0 F1,1 x1 − x
2 F2,0 F2,1 F2,2 x2 − x
3 F3,0 F3,1 F3,2 F3,3 x3 − x
.. ..
. ... .
k Fk,0 Fk,1 Fk,2 Fk,3 Fk,j Fk,j+1 xk − x
.. ..
. .
n Fn,0 Fn,1 Fn,2 Fn,n xn − x

Par exemple, pour calculer Fk,2 (x), on considère le produit croisé suivant : F1,1 (x)(xk − x) −
Fk,1 (x)(x1 − x). Nous obtenons :
F1,1 (x)(xk − x) − Fk,1 (x)(x1 − x)
Fk,2 (x) = .
xk − x1
Nous calcuons les éléments du tableau progressivement afin de déterminer Fn,n (x) = P (x).

Démonstration du théorème 2.2.1 : Par récurrence sur j : Fk,0 (x) = fk pour k = 0..n.
Supposons que Fk,j (x) est le polynôme de Lagrange sur les points x0 , x1 , . . . , xj−1 , xk pour tout
k = j..n.
En particulier, et d’après l’hypothèse de récurrence, Fj,j (x) est le polynôme de Lagrange sur
x0 , x1 , . . . , xj−1 , xj .
Montrons que Fk,j+1 (x) est le polynôme de Lagrange sur x0 , x1 , . . . , xj , xk pour k = j + 1..n.
Pour i = 0..j − 1 :
Fj,j (xi )(xk − xi ) − Fk,j (xi )(xj − xi )
Fk,j+1 (xi ) = = fi
xk − xj
Nous utilisons pour cela l’hypothèse de récurrence : Fj,j (xi ) = fi et Fk,j (xi ) = fi pour i =
0..j−1). D’autre part, en utilisant la définition, on trouve que Fk,j+1 (xj ) = fj et Fk,j+1 (xk ) = fk .
En prenant j = n, on obtient le résultat du théorème.

Remarque 2.2.1. Le polynôme Fk,j (x) avec k ≥ j est le polynôme de Lagrange sur les points
x0 , x1 , . . . , xj−1 , xk de degré inférieur ou égal à j.

Exemple 2.2.1. L’algorithme d’Aitken appliqué à l’exemple 2.1.1 donne :

k/j 0 1 2 xk − x
0 1 −x
1 2 F1,1 = 1 + x 1−x
2
2 -1 F2,1 = 2x + 1 F2,2 = − x2 + 3x
2
+1 −1 − x
x2 x3 x2
3 4 F3,1 = 1 + 3x
2
x
F3,2 = 2 + 2 + 1 F3,3 = 3
− 2
+ 7x
6
+1 2−x

2.2.1 Formule de Newton

Soient n ∈ N et 0 ≤ i < j ≤ n deux entiers distincts. Notons Di,j (x) le polynôme de Lagrange
sur xi , xi+1 , . . . , xj . C’est un polynôme de degré inférieur ou égal à j − i vérifiant Di,j (xk ) = fk
pour i ≤ k ≤ j.

ESSAI - 2017/2018 23 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

2.2.2 La méthode des différences divisées

Proposition 2.2.2 (relation des différences divisées). Notons ci,j le coefficient de xj−i dans
Di,j (x). Alors, on a la relation suivante :

 ci,i = fi pour i = 0..n,
ci+1,j −ci,j−1
 c
i,j = xj −xi
pour 0 ≤ i < j ≤ n.

Démonstration : Soit n ∈ N. Pour tout 0 ≤ i < j ≤ n, on a

Di+1,j (x)(x − xi ) − Di,j−1 (x)(x − xj )
Di,j (x) = . (2.2.1)
xj − xi
En effet, le degré de Di+1,j (x) est ≤ j − i − 1 et le degré de Di,j−1 (x) est ≤ j − i − 1
(x−xi )Di+1,j (x)−(x−xj )Di,j−1 (x)
ce qui montre que le degré du polynôme xj −xi
est ≤ j − i. De plus,
Di+1,j (xk )(xk −xi )−Di,j−1 (xk )(xk −xj )
xj −xi
= fk (il faut distinguer les trois cas : k = i, k = j et k = i+1..j−1
où Di+1,j (xk ) = Di,j−1 (xk ) = fk ).
D (x)(x−x )−D (x)(x−x )
Le polynôme i+1,j i
xj −xi
i,j−1 j
est alors le polynôme de Lagrange sur les points xi , . . . , xj
et puisqu’il est unique (voir théorème 2.1.1) on a l’égalité (2.2.1) qui est équivalente à :

Di+1,j (x)(x − xi ) − Di,j−1 (x)(x − xi )

Di,j (x) = + Di,j−1 (x) .
xj − xi
D (x)(x−xi )−Di,j−1 (x)(x−xi )
On remarque que i+1,j xj −xi
est un polynôme de degré ≤ j − i qui s’annule en
xi , xi+1 , ..xj−1 alors :

Di,j (x) = ci,j (x − xi )(x − xi+1 ) . . . (x − xj−1 ) + Di,j−1 (x) . (2.2.2)

En examinant le coefficient de xj−i de part et d’autre des deux égalités (2.2.1) et (2.2.2), nous
c −ci,j−1
obtenons la relation ci,j = i+1,j
xj −xi
:

Di+1,j (x)(x − xi ) − Di,j−1 (x)(x − xj )

ci,j (x − xi )(x − xi+1 ) . . . (x − xj−1 ) + Di,j−1 (x) = .
xj − xi

En itérant le résultat de l’équation (2.2.2), on obtient :

Di,j (x) = ci,j (x − xi )(x − xi+1 ) . . . (x − xj−1 ) + ci,j−1 (x − xi ) . . . (x − xj−2 ) + · · · + ci,i

et pour i = 0 et j = n, on a :
Définition 2.2.1 (Formule de Newton). Le polynôme de Lagrange P (x) sur les points x0 , x1 , . . . , xn
est donné par la formule de Newton :
k=n
X
P (x) = D0,n (x) = f0 + c0,k (x − x0 )(x − x1 ) . . . (x − xk−1 ) .
k=1

c0,k est le coefficient de xk appelée la différence divisée d’ordre k aux points x0 , x1 , . . . , xk .

ESSAI - 2017/2018 24 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

La formule de Newton permet de déterminer le polynôme de Lagrange. Pour cela, nous utilisons
un tableau analogue à la méthode d’Aitken : les coefficients de P (x) sont sur la diagonale du
tableau triangulaire.
Exemple 2.2.2. En utilisant la méthode des différences divisées et plus précisément la proposi-
tion 2.2.2 dans l’exemple 2.1.1 on obtient :
xj i=j i=j−1 i=j−2 i=j−3
j =0 0 c0,0 = 1
j =1 1 c1,1 = 2 c0,1 = 1
j =2 −1 c2,2 = −1 c1,2 = 32 c0,2 = − 21
j =3 2 c3,3 = 4 c2,3 = 53 c1,3 = 61 c0,3 = 1
3
c −c
Nous remplissons le tableau ligne par ligne avec ci,i = fi pour i = 0..n et ci,j = i+1,j i,j−1
xj −xi
pour 0 ≤ i < j ≤ 3. Ainsi, D0,3 (x) = 31 (x − 0)(x − 1)(x + 1) + (− 21 )(x − 0)(x − 1) + x + 1 =
1 3
3
x − 12 x2 + 67 x + 1.

2.2.3 Algorithme de Neville

L’algorithme de calcul du polynôme de Lagrange P (x) = D0,n (x) pour x fixé est constitué des
n étapes suivantes :

Etape 1 Pour i allant de 0 à n − 1, calculer

fi+1 (x − xi ) − fi (x − xi+1 )
Di,i+1 (x) =
xi+1 − xi
Etape 2 Pour i allant de 0 à n − 2, calculer
Di+1,i+2 (x)(x − xi ) − Di,i+1 (x)(x − xi+2 )
Di,i+2 (x) =
xi+2 − xi
..
.
Etape n − 1 Pour i allant de 0 à 1, calculer
Di+1,i+n−1 (x)(x − xi ) − Di,i+n−2 (x)(x − xi+n−1 )
Di,i+n−1 (x) =
xi+n−1 − xi
Etape n Pour i = 0 Calculer Di,i+n (x) = P (x).

2.3 Complexité d’algorithmes

Que devient le temps de calcul si on multiplie la taille des données par 2 ? Pour répondre à cette
question, on considère le nombre d’opérations élémentaires (affectations, comparaisons, opéra-
tions arithmétiques) effectuées par l’algorithme. Ce nombre est couramment appelé complexité
algorithmique et il s’exprime en fonction de la taille n des données. On dit que la complexité de
l’algorithme est O(h(n)) où h est généralement une combinaison de polynômes, de logarithmes
ou d’exponentielles. Ceci reprend la notation mathématique classique, et signifie que le nombre
d’opérations effectuées est borné par ch(n), où c est une constante, lorsque n tend vers l’infini.

La croissance de cette complexité est en fonction de la taille des données. Les algorithmes usuels
peuvent être classés en un certain nombre de grandes classes de complexité (source wikipédia) :

ESSAI - 2017/2018 25 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

— Les algorithmes sub-linéaires, dont la complexité est en général en O(logn). C’est le cas
de la recherche d’un élément dans un ensemble ordonné fini de cardinal n.
— Les algorithmes linéaires en complexité O(n) ou en O(nlogn) sont considérés comme
rapides, comme l’évaluation de la valeur d’une expression composée de n symboles ou
les algorithmes optimaux de tri.
— Plus lents sont les algorithmes de complexité située entre O(n2 ) et O(n3 ), c’est le cas de
la multiplication des matrices et du parcours dans les graphes.
— Au delà, les algorithmes polynomiaux en O(nk ) pour k > 3 sont considérés comme
lents, sans parler des algorithmes exponentiels (dont la complexité est supérieure à tout
polynôme en n) que l’on s’accorde à dire impraticables dès que la taille des données est
supérieure à quelques dizaines d’unités.

2.4 Estimation de l’erreur

Pour évaluer l’erreur d’interpolation de Lagrange, on a le théorème suivant :
Théorème 2.4.1. Soient f une fonction de classe C n+1 sur un intervalle [a, b] et P le polynôme
de Lagrange de f en les points x0 < x1 < · · · < xn ∈ [a, b]. Alors
(x − x0 )(x − x1 ) . . . (x − xn ) (n+1)
f (x) − P (x) = f (ζ) (2.4.1)
(n + 1)!
où a ≤ min(x, x0 ) < ζ < max(x, xn ) ≤ b.
Démonstration : La relation est triviale pour x = xi ∀i = 0..n. Soit x 6= xi , pour tout i = 0..n.
f (x)−P (x)
On pose k(x) = (x−x0 )(x−x 1 )...(x−xn )
et w(t) un polynôme qui s’annule en x0 , x1 , . . ., xn , et x
définit par : w(t) = f (t) − P (t) − (t − x0 )(t − x1 ) . . . (t − xn )k(x). Il s’annule en (n + 2) points.
En appliquant le thórème de Rolle, on montre que wn+1 (t) s’annule en au moins un point ζ :
f n+1 (ζ)
wn+1 (ζ) = f n+1 (ζ) − (n + 1)! k(x) = 0 et k(x) = .
(n + 1)!

On remarque que l’erreur dépend de la fonction f n+1 (ζ) et des points d’interpolation par (x −
x0 )(x − x1 ) . . . (x − xn ). Considérons le cas où les points d’interpolation sont équidistants, c’est-
à-dire, xi = a + hi pour i = 0..n avec h = b−a n
. Si f est continue et que toutes ses dérivées sont
continues jusqu’à l’ordre n + 1, alors
1 b − a n+1
| f (x) − P (x) |≤ ( ) maxa≤x≤b | f (n+1) (x) | (2.4.2)
2(n + 1) n
Cette quantité ne tend pas nécessairement vers 0 car les dérivées f (n+1) de f peuvent grandir
très vite lorsque n croît. Dans l’exemple suivant, nous présentons deux cas sur la convergence
de l’interpolation de Lagrange.
Exemple 2.4.1. Soit f (x) = sin(x), | f k (x) | ≤ 1 : l’interpolé P converge vers f quelque soit
1
le nombre de points d’interpolation et leur emplacement. Par contre, pour f (x) = 1+25x 2 sur

[−1, 1] et bien que f soit indéfiniment continûment dérivable sur [−1, 1], les grandeurs
max−1≤x≤1 | f k (x) | , k = 1, 2, 3, . . .
explosent très rapidement et l’inégalité (2.4.2) ne nous assure plus la convergence de l’interpo-
lation.

ESSAI - 2017/2018 26 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Le choix des points xi apporte une amélioration sensible de l’interpolation : Lorsque a ≤ x ≤ b,

on choisit les abscisses d’interpolation
a+b b−a 2i + 1
xi = + cos( π)
2 2 2(n + 1)
pour i = 0..n, on obtient :
(b − a)n+1
| f (x) − P (x) | ≤ maxx∈[a,b] | f (n+1) (x) | .
(n + 1)!22n+1
Les points xi = a+b
2
+ b−a
2
2i+1
cos( 2(n+1) π) pour i = 0..n sont calculés à partir des zéros des poly-
nômes de Tchebycheff (Tn+1 (x) = 2xTn (x) − Tn−1 (x)).

Nous venons de voir que l’interpolant de Lagrange peut diverger lorsque le nombre de points
d’interpolation n devient grand et lorsque les points d’interpolation sont équidistants. Une
première solution consiste à découper l’intervalle [a, b] en plusieurs sous-intervalles, c’est l’in-
terpolation par intervalles. On peut également utiliser des points judicieusement placés (zéros
des polynômes de Tchebychev par exemple).
Exercice :
Les polynômes de Tchebychev constituent un outil important dans le domaine de l’interpola-
tion. Soit n ∈ N. Noter que pour minimiser l’erreur d’interpolation de Lagrange, les racines des
polynômes de Tchebychev sont les candidats pour être les points d’interpolations. Le polynôme
de Tchebychev de degré n est noté par Tn (x) = cos(n arccos(x)) pour tout x ∈ [−1, 1] et il
vérifie Tn+2 (x) = 2xTn+1 (x) − Tn (x) avec T0 (x) = 1 et T1 (x) = x.

1. Soit n ∈ N∗ . Montrer par récurrence que le coefficient du terme de degré n de Tn est

2n−1 .
2. Vérifier que Tn admet n racines simples définies par
2k + 1
xk = cos( π), k = 0, 1, . . . , n − 1, n ∈ N∗
2n
3. Sachant que la dérivée de Tn (x) est définie par Tn0 (x) = √ n
1−x2
sin(n arccos(x)). Montrer
que Tn admet exactement (n + 1) extrema définis par
kπ
x0k = cos( ), k = 0, 1, . . . , n.
n
En déduire que Tn (x0k ) = (−1)k , pour k = 0, 1, . . . , n.
4. Vérifier que maxx∈[−1,1] |Tn (x)| = 1 et en déduire que
1
max |(x − x0 )(x − x1 ) . . . (x − xn )| = .
x∈[−1,1] 2n

2.5 Interpolation Spline

Nous rencontrons les splines dans la majorité des programmes de CAO et il est certainement
bon pour un ingénieur d’en connaître l’existence et les principes.

ESSAI - 2017/2018 27 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

L’interpolation spline est une approche locale en faisant de l’interpolation par morceaux : sur
un sous-intervalle donné, on fait de l’interpolation polynomiale de degré faible (≤ 3). Les splines
réalisent une interpolation polynomiale par morceaux et on imposera aux noeuds un degré de
régularité suffisant. Avec les splines, on obtient une très bonne approximation, sans effets de
bord, contrairement à l’interpolation de Lagrange.

Etant donnée une fonction f définie sur [a, b] et un ensemble de noeuds a ≤ x0 < x1 < · · · <
xn ≤ b. Nous appelons spline d’ordre 3 ou spline cubique interpolant f aux noeuds, une fonction
S vérifiant les conditions suivantes :
1. S coïncide avec un polynôme Pi de degré 3 sur chacun des intervalles [xi−1 , xi ] pour
i = 1..n.
2. S(xi ) = Pi (xi ) = fi pour i = 0..n.
3. Pi (xi ) = Pi+1 (xi ) pour i = 0..n − 1, (continuité de S).
0 0 0
4. Pi (xi ) = Pi+1 (xi ) pour i = 0..n − 1, (continuité de S ).
00 00 00
5. Pi (xi ) = Pi+1 (xi ) pour i = 0..n − 1, (continuité de S ).
6. S satisfait une des conditions au bord suivantes :
00 00 0 0 0 0
S (x0 ) = S (xn ) = 0 (spline libre) ou bien S (x0 ) = f0 et S (xn ) = fn (spline forcée
lorsque f est dérivable).
Proposition 2.5.1. Les polynômes Pi (x) sont données par les formules suivantes :
(xi − x)((xi − x)2 − h2i ) (x − xi−1 )((x − xi−1 )2 − h2i ) fi−1 (xi − x) fi (x − xi−1 )
Pi (x) = Mi−1 + Mi + +
6hi 6hi hi hi

avec hi = xi −xi−1 pour i = 1..n, M0 = Mn = 0 et M1 , M2 , . . . , Mn−1 sont solutions du système

linéaire :
hi hi + hi+1 hi+1 fi+1 − fi fi − fi−1
Mi−1 + Mi + Mi+1 = − pour i = 1..n − 1.
6 3 6 hi+1 hi
00
Indications sur la preuve : Pour montrer ces égalités, on commence par exprimer Pi en
00 00 00 00 i+1 −x 00
fonction de fi et fi+1 par une interpolation linéaire : Pi (x) = fi xxi+1 −xi
x−xi
+ fi+1 xi+1 −xi
. En
−x) 3 3
“ (xi −x)
intègrant deux fois, on obtient : Pi (x) = fi“ (xi+1
6hi
− fi+1 6hi
+ ai (xi+1 − x) − bi (xi − x) où
ai et bi sont obtenues par les conditions d’interpolation :
fi 00 hi fi+1 00 hi
ai = − fi , bi = − fi+1
hi 6 hi 6
0
En tenant compte des raccords aux noeuds pour Pi et i = 1..n − 1, on obtient un système
linéaire de (n − 1) équations et (n + 1) inconnues.
hi hi + hi+1 hi+1 fi+1 − fi fi − fi−1
Mi−1 + Mi + Mi+1 = − pour i = 1..n − 1.
6 3 6 hi+1 hi

Exemple 2.5.1. On reprend l’exemple 2.1.1 de la section précédente. On cherche la spline cubique
S(x) telle que S(0) = 1, S(1) = 2, S(−1) = −1 et S(2) = 4. On pose x0 = −1, x1 = 0, x2 = 1
et x3 = 2. Ils sont bien ordonnées. Alors f0 = −1, f1 = 1, f2 = 2, f3 = 4 et hi = 1 pour i = 1..3.
On prend M0 = M3 = 0 ; M1 et M2 sont solution du système :
 2
 3 M1 + 16 M2 = −1
1 2
M1 + M2 = 1

6 3

ESSAI - 2017/2018 28 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Nous obtenons :

 P1 (x) = − 13 x3 − x2 + 34 x + 1 x ∈ [−1, 0]
S(x) = P2 (x) = 23 x3 − x2 + 43 x + 1 x ∈ [0, 1]
1 3 2 5
P3 (x) = − 3 x + 2x − 3 x + 2 x ∈ [1, 2]


2.6 Introduction à l’Approximation

Soit f une fonction réelle définie soit de façon discrète, soit de façon continue. Le problème
général de l’approximation consiste à déterminer une fonction g de R dans R de forme donnée,
qui, en un certain sens, approche le mieux possible la fonction f . Un cas particulier est l’in-
terpolation polynomiale (la fonction g est alors un polynôme). Mais l’interpolation est un outil
de peu de valeur pour le traitement des données expérimentales. On l’utilise seulement lorsque
les données sont très précises. Cette section traite particulièrement de la méthode des moindres
carrés discrets.

2.6.1 Résultats fondamentaux

Il est relativement fréquent d’avoir à ajuster sur des données une loi empirique ou encore d’avoir
à déterminer la valeur des paramètres d’une loi théorique de façon à reproduire les résultats
d’une expérience. Une façon possible pour aborder de tels problèmes est de chercher les valeurs
de paramètres minimisant un critère.

Soient u0 , u1 , . . ., up , p + 1 élémets d’un espace vectoriel normé E. Pour tout y ∈ E, le problème

de la meilleure approximation Pde y admet une solution, c’est-à-dire que la fonction norme de
l’erreur (a0 , . . . , an ) =k y − ni=0 ai ui k2 admet un minimum a dans E. Pour cela, on montre
que cette fonction est continue sur un fermé borné : elle atteint donc son minimum.

Notons V le sous-espace vectoriel engendré par les p + 1 points u0 , u1 , . . . , up de E P

deux à deux
distincts. La meilleure approximation de y sur V au sens des moindres carrés est i=p i=0 ai ui où
les ai , i = 0..p sont solutions du système linéaire :
i=p
X
ai < ui , uj > = < y, uj > pour j = 0..p. (2.6.1)
i=0

2.6.2 Approximation polynomiale discrète

Il y a une manière assez simple de trouver l’équation d’une courbe. On positionne les données
sur un graphique, on deplace la règle jusqu’à ce que la ligne se trouve équidistante de tous les
points. L’analyse numérique change cette supposition à vue d’oeil en un processus plus scien-
tifique c’est l’approximation au sens des moindres carrés.

Notons E = R[x] l’espace vectoriel des polynômes en la variable x et à coefficient dans R. Soit
V l’espace vectoriel des polynômes de degré ≤ p. On pose u0 = xp , u1 = xp−1 , . . . , up−1 = x et
up = 1.

ESSAI - 2017/2018 29 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

On considère les abscisses x0 , x1 , . . . , xn ∈ R et f une fonction de R dans R telle que f (xi ) = yi

pour i = 0..n. On veut approcher l’ensemble discret des points (xi , yi )0≤i≤n par un polynôme
de degré p à coefficient dans R :
Pp (x) = a0 xp + a1 xp−1 + · · · + ap−1 x + ap .
Pour tout g1 , g2 ∈ V , on définit le produit scalaire sur V par
i=n
X
< g1 (x), g2 (x) >= g1 (xi )g2 (xi ) .
i=0

Le système (2.6.1) est alors équivalent au système suivant :

 Pi=n 2p Pi=n 2p−1 Pi=n 2p−2 Pi=n p 
a0

p
i=0 xi i=0 xi i=0 xi ... ... i=0 xi
P
y i xi
 
  
  a1
 P Pi=n p−1     
i=n 2p−1 Pi=n 2p−2
p−1
i=0 xi i=0 xi ... ... ... i=0 xi
   P 







 yi xi 

 . 
 .
 
.

 P
i=n 2p−2 .  .
 
.

i=0 xi .
   


 .
 =
  . 

 . .  .
   . 
 . .  .
  . 
. .
P .
   
  .   
. yi xi
    
 .
.
   
 Pi=n  .   
 . x
i=0 i  .  P
Pi=n p Pi=n p−1 Pi=n yi
i=0 xi i=0 xi ... ... i=0 xi p+1 ap

Lorsque p = 1, on parle de régression et dans ce cas, on cherche une droite P1 (x) = a0 x + a1

tel que
Xi=n i=n
X
2
(yi − a0 xi − a1 ) ≤ (yi − α0 xi − α1 )2 ∀α0 , α1 ∈ R .
i=0 i=0

Le système (2.6.1) est alors équivalent à :

 Pi=n 2 Pi=n     Pi=n 
i=0 xi i=0 xi a0 i=0 xi yi
  =  . (2.6.2)
Pi=n Pi=n
i=0 xi n+1 a1 i=0 yi

En posant X = (x0 , x1 , . . . , xn ) et Y = (y0 , y1 , . . . , yn ), on obtient à l’aide des formules de la

covariance et de la moyenne :
Pi=n
a0 = σ(X,Y
σ(X 2 )
)
avec σ(X, Y ) = n+1 1
i=0 (xi − X̄)(yi − Ȳ )
1
P i=n
a1 = Ȳ − a0 X̄ et X̄ = n+1 i=0 xi

Le point moyen (X̄, Ȳ ) appartient à la droite d’équation P1 (x) = a0 x + a1 .

Exemple 2.6.1. On considère la fonction f (x) = x2 + 2x sur [−1, 1]. Calculer les valeurs de f
aux points d’abscisses −1, − 21 , 0, 12 et 1. Donner la droite de meilleure approximation discrète de
f pour ces 5 points au sens des moindres carrés. On note x0 = −1, x1 = − 12 , x2 = 0, x3 = 12 et
x4 = 1. Alors y0 = −1, y1 = − 34 , y2 = 0, y3 = 45 et y4 = 3. On calcule i=4 x2i = 52 , i=4
P P
Pi=4 Pi=4 i=0 i=0 xi = 0,
5
i=0 xi yx = 5 et i=0 yi = 2 . On obtient le système suivant à partir du système (2.6.2) :
 Pi=4 2 5 Pi=n    Pi=n 
i=0 xi = 2 i=0 xi = 0 a0 i=0 xi yi = 5
  =  .
Pi=n Pi=n 5
i=0 xi = 0 n+1=5 a1 i=0 yi = 2

et a0 = 2, a1 = 12 . Donc P1 (x) = 2x + 12 .

ESSAI - 2017/2018 30 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

2.7 Exercices d’Applications

1. Implémenter l’algorithme de Neville dans le logiciel de calcul SAGEmath pour x fixé.
2
2. Construire le polynôme de Lagrange de f (x) = e−x sur les entiers de [−2, 2] par la
méthode des différences divisées.
1
3. Soit la fonction f de R dans R définie par f (x) = 1+x2
.
(a) Construire le polynôme de Lagrange P de f par la méthode des différences divisées
sur 0, 1, 3 et 5.
(b) Comparer P (4) et f (4).
(c) Que peut-on conclure ?
4. On recueille les données suivantes concernant la météo dans la région de Bizerte le 13
janvier 2003 :
temps (heures) 0h00 5h00 7H00 12h00
pluviométrie (millimètres) 10 3 5 2
(a) Utiliser la méthode d’Aitken pour calculer le polynôme interpolant la fonction “plu-
viométrie”.
(b) Approximer la pluviométrie dans la région de Bizerte au temps t=10h00.
5. Soit f (x) = x + exp(−x2 ). Donner la formule du polynôme d’interpolation de Lagrange
Li (x) sur les points (i, f (i)) pour i fixé.
6. Poser f (x) = exp (−x2 ), a = −1, b = 2 et
a+b b−a 2k + 1
xk = + cos( π)
2 2 2(n + 1)
pour k allant de 1 à 3.
(a) Tracer l’ensemble des zéros du polynômes de Tchebycheff sur un même graphe que
la fonction f .
(b) Calculer le polynôme de Lagrange en ces points et comparer avec la fonction f .
7. Déterminer la spline cubique des points (−1, 0), (0, 1), (1, 3) et (2, 2). Comparer cette
spline avec le polynôme de Lagrange en ces même points.
8. L’évolution de la concentration c d’une espèce dans un mélange est donnée par une loi
linéaire en fonction du temps : c(t) = a0 t + a1 . Un expérimentateur fait une série de
mesures pour déterminer les paramètres inconnus :
0 1 2 3 4 5 6 7 8 9 10 11
1,54 3,09 4,97 7,43 10,65 14,92 20,6 28,2 38,42 52,15 70,65 96,6
Pi=n 2 Pi=n Pi=n Pi=n
(a) On pose n = 11. Calculer i=0 ti , i=0 ti , i=0 ti ci et i=0 ci et écrire le système
pour l’approximation au sens des moindres carrées.
(b) Utiliser la méthode de régression des moindres carrés pour obtenir les paramètres a0
et a1 . Rappelons que
σ(X, Y )
a0 = et a1 = Ȳ − a0 X̄ .
σ(X 2 )
9. (a) Trouver le polyôme P du second degré passant par (0, 0), (1, 1), (2, 8).

ESSAI - 2017/2018 31 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

(b) La fonction x3 prend les valeurs de P . Utiliser la formule vue en cours pour exprimer
l’erreur d’interpolation commise en remplaçant x3 par le polynôme d’interpolation
P.
1
10. (a) Interpoler ex aux points 0, 2
et 1 par un polynôme du second degré.
(b) Trouver une expression pour l’erreur d’interpolation et une borne de l’erreur indé-
pendante de ζ.
11. (a) Sachant que la température initiale d’un objet est de 0o C, puis de 5o C après une
minute et 3o C après deux minutes, quelle a été la température maximale et quant
a-t-elle été atteinte ?
(b) Quelle est la température moyenne entre 1 et 2 minutes ?
12. Considérons l’ensemble des points
x0 = −1, x1 = 0, x2 = 1, x3 = 2,
f0 = 7, f1 = 2, f2 = −1, f3 = −14.
(a) Calculer les polynômes d’interpolation de Lagrange L0 et L3 .
(b) Vérifier que Li (xi ) = 1 et que Li (xj ) 6= 0 pour tout i et j 6= i de 0 à 3.
(c) Citer deux méthodes permettant de calculer le polynôme de Lagrange.
(d) Vérifier que le polynôme de Lagrange sur les points (x0 , f0 ), (x1 , f1 ), (x2 , f2 ) et (x3 , f3 )
est égal à :
P (x) = −2x3 + x2 − 2x + 2 .
(e) Afin de déterminer une approximation du nuage de points (xi , fi ) pour i de 0 à 3,
calculer la droite de régression linéaire P0 par la méthode des moindres carrés discrets.
(f) Quelle est la différence entre P et P0 .

Méthode d’Hermite. Soit x0 , x1 , . . . , xn , (n + 1) points de a, b et f une fonction de classe C 1

sur a, b. Au lieu de chercher à faire coïncider f et Pn aux points xi , nous pouvons aussi chercher
à faire coïncider f et Pn ainsi que leurs dérivées jusqu’à un ordre donné aux points xi .
De la même manière que dans l’interpolation de Lagrange, nous cherchons un polynôme Pn tel
que Pn (xi ) = f (xi ) et Pn0 (xi ) = f 0 (xi ).

ESSAI - 2017/2018 32 inestej@[Link]

Chapitre 3

Intégration numérique

3.1 Introduction
Soit f une fonction continue d’un intervalle [a, b] dans R. On désire approcher numériquement
la quantité Z b
f (x) dx .
a

Pour cela, on subdivise l’intervalle [a, b] en N intervalles : xi = a + ih, i = 0..N − 1 et h = b−aN

.
Rb xi+1 −xi
On peut alors approcher a f (x) dx par la somme des airs des rectangles de longueur f ( 2 )
et de largeur xi+1 − xi . C’est la méthode des rectangles.
Exercice 3.1.1. Calculer l’aire comprise entre l’axe des abscisses, l’axe des ordonnées, la droite
2
x = 1 et la courbe f (x) = ex avec la méthode du rectangle.

La méthode de Newton-Côtes utilise les polynômes de Lagrange P (x) pour approcher la fonction
f (x) sur [a, b] ; on considère comme valeur approchée de l’intégrale de f entre a et b :
Z b
P (x) dx .
a

3.2 Formules de Newton-Côtes

Soit n ∈ N∗ . On considère la subdivision uniforme xi = a + ih avec i = 0..n et h = b−a n
. Soit
P le polynôme d’interpolation de Lagrange de degré ≤ n vérifiant P (xi ) = f (xi ), i = 0..n. On
Rb Rb RbP
obtient alors a f (x) dx ' a P (x) dx = a i=n i=0 f (xi ) Li (x) dx. En utilisant le changement de
variables x = a + ht dans les Li (x), on obtient :
Z b i=n
X Z n k=n
Y t−k
P (x) dx = h f (xi ) αi et αi = dt .
a i=0 0 i=0,k6=i i − k

Les nombres rationnels αi sont indépendants de f et de [a, b].

Remarque 3.2.1. Si f = 1, a = 0 et b = 1, on obtient que le polynôme de Lagrange associé à f
soit égal à 1, d’où :
Z 1 i=n
X i=n
X
1= dx = h f (xi ) αi = h αi .
0 i=0 i=0

33
Cours d’Analyse Numérique Ines Abdeljaoued Tej

Pi=n
Ce qui donne i=0 αi = n.
Remarque 3.2.2. Pour n = 1, on obtient α0 = α1 = 12 , d’où
b
b−a
Z
f (x) dx = (f (a) + f (b)) .
a 2

Pour n = 2, on obtient α0 = 31 , α1 = 4
3
et α2 = 13 . D’où :
b
b−a
Z
f (x) dx = (f (a) + 4f (a + h) + f (b)) .
a 6

3.2.1 Méthode du trapèze

Soit N ∈ N∗ . On considère la subdivision uniforme xi = a + ih, i = 0..N avec h = b−a
N
. La
méthode du trapèze consiste à appliquer la méthode de Newton-Côtes avec n = 1 sur chaque
intervalle [xi , xi+1 ], i = 0..N − 1. On obtient :
Z b N −1 Z xi+1 N −1
X X h
f (x) dx = f (x) dx ' T (h) et T (h) = (f (xi ) + f (xi+1 ))
a i=0 xi i=0
2

Théorème 3.2.1. La méthode du trapèze est d’ordre deux. Plus précisément, si f est suffisam-
ment dérivable et h = b−a
N
, alors il existe c ∈]a, b[ tel que
b
b − a 2 00
Z
f (x) dx − T (h) = − h f (c) .
a 12
Exercice 3.2.1. Reprendre l’exemple 3.1.1 avec la méthode du trapèze.

3.2.2 Méthode de Simpson

Soit N ∈ N∗ . On considère la subdivision uniforme xi = a + ih, i = 0..2N avec h = b−a
2N
. La
méthode de Simpson consiste à appliquer la méthode de Newton-Côtes avec n = 2 sur chaque
intervalle [x2i , x2i+2 ], i = 0..N − 1. On obtient :
Z b N
X −1 Z x2i+2
f (x) dx = f (x) dx ' S(h)
a i=0 x2i

et
N −1
X h
S(h) = (f (x2i ) + 4f (x2i+1 ) + f (x2i+2 ))
i=0
3

Théorème 3.2.2. La méthode de Simpson est d’ordre quatre. Plus précisément, si f est suffi-
samment dérivable et h = b−a
2N
, alors il existe c ∈]a, b[ tel que
b
b − a 4 (4)
Z
f (x) dx − S(h) = − h f (c) .
a 180
Exercice 3.2.2. Reprendre l’exemple 3.1.1 avec la méthode de Simpson.

ESSAI - 2017/2018 34 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

3.2.3 Remarque sur la Formule Composite

Les méthodes du rectangle, du trapèze et de simpson sont des cas particuliers de la formule
composite définie par :
N −1 j=n
X xi+1 − xi X tj + 1
αj f (xi + (xi+1 − xi ) ) .
i=0
2 j=0
2

Les n + 1 points −1 ≤ t0 ≤ · · · ≤ tn ≤ 1 sont appelés points d’intégration et les n + 1 points

α0 , . . . , αn sont appeés poids de la formule de quadrature.
Remarque 3.2.3. Pour la méthode du rectangle, on prends n = 0, t0 = 0 et α0 = 2. La formule
composite devient
N −1
X xi+1 − xi xi + xi+1
f( ) .
i=0
2 2
Remarque 3.2.4. Pour la méthode du trapèze, on prends n = 1, t0 = −1, t1 = 1, α0 = 1 et
α1 = 1. La formule composite devient
N −1
X xi+1 − xi
(f (xi ) + f (xi+1 )) .
i=0
2

Remarque 3.2.5. Pour la méthode de Simpson, on prends n = 2, t0 = −1, t1 = 0, t2 = 1,

α0 = 13 , α1 = 43 et α2 = 13 . La formule composite devient
N −1
X xi+1 − xi xi + xi+1
(f (xi ) + 4f ( ) + f (xi+1 )) .
i=0
6 2

3.3 Exemples d’Implémentations

Voir en annexe, les simulations numériques sur ces méthodes.

ESSAI - 2017/2018 35 inestej@[Link]

Chapitre 4

Racine d’équation non linéaire

Les méthodes analytiques de résolution d’équations algébriques polynomiales sont limités à cer-
taines formes et à de faible degré (équations quadratiques, cubiques, quartiques, ax2n +bxn +c =
0, etc...). La résolution des équations polynomiales de degré 4 est déjà plus pénible et à partir
du degré 5, E. Galois a montré qu’il n’est pas toujours possible de donner une forme exacte des
racines au moyen des coefficients (on utilise alors le groupe de l’équation, l’ideal des relations
entre les racines ou encore le corps de décomposition du polynôme : ce sont des méthodes for-
melles [?]).

Pour la résolution des équations polynomiales, nous sommes amenés à l’utilisation des méthodes
numériques. Il en est de même pour les équations contenant des fonctions transcendantes (i.e.
non-algébriques) telles que exp(x), ch(x), etc... En étudiant les méthodes de résolution d’une
équation f , on se posera deux questions :
1. La méthode converge-t-elle vers x∗ où f (x∗ ) = 0. En d’autre termes la suite (xn )n∈N
générée par la méthode converge-t-elle vers la solution x∗ .
2. Dans l’affirmative, avec quelle rapidité (ceci permet de comparer les méthodes entre-
elles).

4.1 Quelques algorithmes classiques

On désigne par (xn ) une suite de nombres réels x0 , x1 , . . . , xn , . . . . Si une telle suite converge
vers un point x∗ , on écrira (xn ) −→ x∗ . En analyse numérique, on construit les suites à l’aide
d’algorithmes.

4.1.1 Méthode de la bissection

On considère une fonction f continue et on cherche x∗ solution de f (x) = 0. On suppose
qu’on a localisé par tâtonnement un intervalle [a, b] dans lequel la fonction f change de signe :
f (a).f (b) < 0.

On est donc certain qu’il y a entre a et b au moins un zéro de f . Pour approcher de façon
précise l’un de ces zéros, on utilise l’algorithme suivant :

Algorithme 1 (Dichotomie).

36
Cours d’Analyse Numérique Ines Abdeljaoued Tej

Entrée : les extrémitées a et b,

la fonction f, la précision  souhaitée et
N le nombre d’itérations maximum.
Sortie : Une approximation c de la racine x∗ de f
ou bien un mesSAGEmath d’erreur.
1. Si f (a).f (b) > 0 Alors Imprimer(pas de solution).
Sinon
n := 1 ;
2. Tant que n ≤ N Faire
3. c := a+b
2
;
4. Si f (c) = 0 ou b−a2
≤  Alors
5. Imprimer(c) ; n := N + 2 ;
Fin Si ;
6. n := n + 1 ;
7. Si f (a).f (c) > 0 Alors a := c ; Sinon b := c ;
Fin Tant que ;
Si n = N + 2 Alors Imprimer(c)
Sinon Imprimer(Après N itérations l’approximation de x∗
obtenue est c et l’erreur maximale est b−a
2
);
Fin Si ;
Fin Si ;
Fin.

A chaque itération, l’algorithme construit un nouvel intervalle, autour de x∗ , qui est de longueur
égale à la moitié de la longueur de l’intervalle précédent. Au bout de n itérations, la longueur
de l’intervalle est alors de b−a
2n
.

Le nombre d’itérations nécessaire pour obtenir une précision de calcul égal à  est définie par :
log (b − a) − log 
n≥ .
log 2
Remarque 4.1.1. Les avantages de cet algorithme sont la convergence assurée et une marge
d’erreur sûre. Par contre il est relativement lent. Il est souvent utilisé pour déterminer une
approximation initiale à utiliser avec un algorithme plus rapide.

4.1.2 Méthode de Newton

Cet algorithme se base sur la suite suivante (donner une explication géométrique) :
(
x0 donné
xn+1 = xn − ff0(x n)
(x )
.
n

La signification géométrique consiste à mener par le point (xn , f (xn )) la tangente à la courbe
f (x). La pente de cette tangente est
0
T (x) = f (xn ) + (x − xn )f (xn ) .

ESSAI - 2017/2018 37 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Si on cherche le point d’intersection de la tangente avec l’axe des x, i.e. si on résoud T (x) = 0,
on retrouve xn+1 tel que défini par la suite.

Théorème 4.1.1. [Convergence globale de la méthode de Newton] Soit f une fonction de Classe
C 2 sur [a, b] et g(x) = x − ff0(x)
(x)
. Si f vérifie :
1. f (a).f (b) < 0
0
2. ∀x ∈ [a, b], f (x) 6= 0 (c’est la stricte monotonie),
00
3. ∀x ∈ [a, b], f (x) 6= 0 (concavité dans le même sens).
00
Alors en choisissant x0 ∈ [a, b] tel que f (x0 ).f (x0 ) > 0, la suite (xn ) définie par x0 et xn+1 =
g(xn ) converge vers l’unique solution de f (x) = 0 dans [a, b].

Preuve : Les deux premières conditions du théorème assurent l’existence et l’unicité d’une
racine x∗ unique de f dans [a, b]. D’autre part, comme xn+1 = xn − ff0(x n)
(xn )
on obtient :

f (xn )
xn+1 − x∗ = xn − x∗ −
f 0 (xn )
f (x∗ ) − f (xn )
= xn − x∗ +
f 0 (xn )
∗ −x 2
(x∗ − xn )f 0 (xn ) + (x
∗ 2
n)
f 00 (ζn )
= xn − x +
f 0 (xn )

Ainsi :
∗
∗ ∗ f 0 (xn ) + x −x
2
n 00
f (ζn )
xn+1 − x = (xn − x )(1 − 0
)
f (xn )
(xn − x∗ )2 f 00 (ζn )
=
2 f 0 (xn )

Si f 00 et f 0 sont de même signe, alors xn+1 − x∗ > 0 et la suite est alors minorée par x∗ à partir
du rang 1. Sinon, si f 00 et f 0 sont de signes contraires, alors xn+1 − x∗ < 0 et la suite est alors
majorée. De plus, et selon les cas, la suite va être soit décroissante (f 00 > 0 et f 0 > 0 ou bien
f 00 < 0 et f 0 < 0) donc convergente soit croissante (f 00 > 0 et f 0 < 0 ou bien f 00 < 0 et f 0 > 0)
et donc convergente.

Algorithme 2 (Newton).

ESSAI - 2017/2018 38 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Entrée : Une approximation initiale x0 , la précision  souhaitée

et N le nombre maximum d’itérations.
Sortie : Une approximation c de la racine x∗ de f ou bien un mesSAGEmath d’erreur.
1. n := 1 ;
2. Tant que n ≤ N Faire
c := x0 − ff0(x0)
(x0 )
;
3. Si | c − x0 |≤  Alors c ; n := N + 2 Fin Si ;
n := n + 1 ;
x0 := c ;
Fin Tant que ;
Si n = N + 2 Alors Imprimer(c) Sinon
Imprimer(La méthode a échoué après N itérations) ;
Fin Si ;
Fin.

Remarque 4.1.2. Il est essentiel de fixer une limite au nombre d’itérations car la convergence
n’est pas assurée. On verra dans la suite qu’on peut construire des exemples pour lesquels cet
algorithme ne converge pas. Cependant lorsque la convergence aura lieu, elle sera trés rapide
comme le montre le théorème 4.1.2. D’autre part, le test d’arrêt n’assure pas que la racine x∗
soit à une distance  de x.

en+1
Définition 4.1.1. La méthode définie par xn+1 = g(xn ) est dite d’ordre p si la limite de | epn
|
lorsque n −→ +∞ est une constante réelle strictement positive.

Théorème 4.1.2. Soit f une fonction de Classe C 2 . Si x∗ est une racine simple de f , alors la
méthode de Newton est au moins d’ordre 2.
f (x).f 00 (x)
Preuve : On remarque que g 0 (x) = (f 0 (x))2
donc g 0 (x∗ ) = 0. D’autre part :

xn+1 − x∗ = g(xn ) − g(x∗ )

00 ∗
∗ 2 g (x )
= (xn − x ) + (xn − x∗ )2 (en )
2
00 ∗
g (x )
= e2n + e2n (en )
2
Comme (en ) est négligeable au voisinage de +∞ on obtient le résultat. A chaque itération, on
multiplie par 2 le nombre de chiffres exacts de l’approximation.
2 p
D’une façon générale, si f est de Classe C p alors : en+1 = en g 0 (x∗ ) + e2!n g 00 (x∗ ) + · · · + ep!n g (p) (x∗ ) +
epn (en ). La méthode est alors d’ordre p si g 0 (x∗ ) = · · · = g (p−1) (x∗ ) = 0 et g (p) (x∗ ) 6= 0.

Il est parfois ennuyeux dans l’utilisation de la méthode de Newton-Raphson d’avoir à calculer

0
f (x). Il se peut par exemple qu’on ne dispose pas du programme d’ordinateur pour en effectuer
le calcul alors qu’on peut facilement calculer f (x). L’algorithme suivant peut être considéré
comme une approximation de la méthode de Newton.

ESSAI - 2017/2018 39 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

4.1.3 Méthode de la sécante

Au lieu d’utiliser la tangente au point xn , nous allons utiliser la sécante passant par les points
d’abscisses xn et xn−1 pour en déduire xn+1 (donner une explication géométrique). L’équation
de la sécante s’écrit sous la forme :
f (xn ) − f (xn−1 )
S(x) = f (xn ) + (x − xn )τn avec τn = .
xn − xn−1

On définit xn+1 comme étant l’intersection de la sécante S(x) avec l’axe des x : S(xn+1 ) = 0 et
on obtient
xn − xn−1
xn+1 = xn − f (xn ) .
f (xn ) − f (xn−1 )
Remarque 4.1.3. La méthode de la sécante est une méthode itérative à deux pas en ce sens que
xn+1 dépend de xn et de xn−1 . Dans la méthode de Newton-Raphson xn+1 ne dépend que de
xn . C’est une méthode à un pas.
Une façon commode de procéder est de choisir x0 arbitrairement puis de choisir x1 = x0 + h
où h est un accroissement petit. Si xn et xn+1 sont proches (on peut espérer que ce sera le cas
0
si la méthode converge), on peut considérer que f (xxnn)−f (xn−1 )
−xn−1
est une approximation de f (xn ).
Dans ce sens l’algorithme de la sécante est une approximation de l’algorithme de Newton.

Algorithme 3 (Algorithme de la sécante).

Entrée : Deux approximations initiale x0 et x1 ,

la précision  souhaitée et N le nombre maximum d’itérations
Sortie : Une approximation c de la racine x∗ de f ou bien un mesSAGEmath d’erreur.

1. n := 2 ;
q0 := f (x0 ) ;
q1 := f (x1 ) ;
2. Tant que n ≤ N + 1 Faire
c := x1 − q1 xq11 −x
−q0
0
;
3. Si | c − x1 |≤  Alors Imprimer(c) ; Fin.
4. n := n + 1 ;
x0 := x1 ; x1 := c ;
q0 := q1 ; q1 := f (c) ;
Fin Tant que ;
Imprimer(La méthode a échoué après N itérations) ;
Fin.

4.1.4 Méthode du point fixe

On peut toujours transformer un problème du type f (x) = 0 en un problème de la forme
g(x) = x et ce d’une infinité de façons. Par exemple, la méthode de Lagrange est donnée
par xn+1 = g(xn ) où g(x) = aff (x)−xf (a)
(x)−f (a)
(on pose alors x0 = b). Il faut toutefois noter que de
telles transformations introduisent parfois des solutions parasites. Par exemple, on peut écrire
f (x) = x1 − 3 = 0 ou encore x = 2x − 3x2 = g(x) et 0 est solution de la seconde équation mais

ESSAI - 2017/2018 40 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

pas de la première.

Le théorème suivant caractérise les fonctions g qui donnent des suites convergentes :

Théorème 4.1.3. Si dans l’intervalle [a, b], g vérifie les conditions suivantes :
1. x ∈ [a, b] alors g(x) ∈ [a, b],
0
2. l’application g est strictement contractante, c’est-à-dire maxx∈[a,b] | g (x) |= L < 1 .
Alors pour tout x0 ∈ [a, b], la suite définie par xn+1 = g(xn ) converge vers l’unique point fixe
x∗ de g sur [a, b].

Il est souvent délicat de déterminer un intervalle [a, b] dans lequel les hypothèses du théorème
0
4.1.3 sont vérifiées. Si on peut estimer | g (x∗ ) |, on a le résultat suivant :

0 0
Théorème 4.1.4. Soit x∗ une solution de l’équation x∗ = g(x∗ ). Si g est continue et | g (x∗ ) |<
1 alors il existe un intervalle [a, b] contenant x∗ pour lequel la suite définie par x0 ∈ [a, b] et
xn+1 = g(xn ) converge vers x∗ .

0
Proposition 4.1.5. Soit x∗ une solution de l’équation x∗ = g(x∗ ). Si g est continue au voisi-
0
nage de x∗ et si | g (x∗ ) |> 1 alors pour tout x0 ∈ [a, b] différent de x∗ , la suite xn+1 = g(xn )
ne converge pas vers x∗ .

4.2 Méthodes d’accélération de la convergence

On veut savoir si la suite définie par xn+1 = g(xn ) converge rapidement et comment se comporte
l’erreur en = xn − x∗ d’une itération à la suivante ?

4.2.1 Procédé ∆2 d’Aitken

Dans les méthodes d’ordre 1 (i.e. g 0 (x∗ 6= 0)), on considère la suite
(xn+1 − xn )2
y n = xn − .
xn+2 − 2xn+1 + xn
C’est une meilleure approximation de x∗ que xn . Lorsque xn est proche de la solution x∗ ,
(xn+1 − xn )2 et xn+2 − 2xn+1 + xn sont voisins de 0. Le calcul de yn est alors instable. On utilise
alors une autre écriture équivalente à la précédente mais plus stable :
1
yn = xn+1 + 1 .
xn+2 −xn+1
− xn+11−xn

4.2.2 Méthode de Steffensen

On poursuit les itérations en remplaçant xn par yn et on calcule ainsi g(yn ) = zn et g(g(yn )) =
xg(g(x))−g(x)2
g(zn ). On construit formellement la méthode de Steffensen zn+1 = G(zn ) où G(x) = g(g(x))−2g(x)+x .
La méthode de Steffensen est d’ordre au moins 2 si (xn ) est d’ordre 1 et si (xn ) est d’ordre p
alors (zn ) est d’ordre 2p − 1.

ESSAI - 2017/2018 41 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

4.3 Exercices d’Applications

1. Soit f (x) = x2 − exp (−1 − x2 ) une fonction définie sur [0, 5].
(a) Tracer f , f 0 et f 00 sur [0, 5].
(b) Montrer qu’il existe une racine unique de f sur [0, 5].
(c) Appliquer la méthode de Newton avec x0 = 5.
(d) Est-ce que la méthode de Newton converge pour x0 = 0 ? justifier votre réponse.
(e) Appliquer la méthode de la sécante à f avec x0 = 0. Que remarquez-vous ?
(f) Trouver une fonction g vérifiant g(x) = x si et seulement si f (x) = 0 qui assure les
conditions du théorème du point fixe. Appliquer la méthode des itérations successives.
(g) Citer deux procédés pour accélérer la convergence des algorithmes.
2. Refaire l’exercice précédent avec f (x) = x3 − 4x + 1 sur [0, 0.5] et g(x) = x − f (x).
3. (a) Donner le théorème de convergence globale de la méthode de Newton pour une fonc-
tion de classe C 2 .
(b) Ecrire l’algorithme de Newton qui prend en entrée les points x0 , a < b et une fonction
f et rend la racine de f (x) = 0 sur [a, b] ou bien un mesSAGEmath d’erreur.
(c) Calculer la racine de f (x) = −2x3 + x2 − 2x + 2 sur l’intervalle [0.5, 1] avec x0 = 1.
La précision des calculs est à 10−6 près.
4. Soit f une fonction de R dans R telle que
f (x) = x3 + 3x2 − 1 pour x ∈ R.
(a) Montrer que l’équation f (x) = 0 admet une unique solution sur [0, 1].
(b) Nous désirons déterminer la solution de l’équation f (x) = 0 avec la méthode de
Newton sur [0, 1]. En quoi consiste cette méthode ?
(c) Ecrire l’algorithme général de la méthode de Newton qui prend en entrée les points
x0 , a < b et une fonction f et rend la racine de f (x) = 0 sur [a, b] ou bien un
mesSAGEmath d’erreur.
(d) Est-ce que la fonction f vérifie les conditions du théorème de convergence globale de
la méthode de Newton pour x0 = 1. Donner le théorème.
(e) Calculer la racine de f sur l’intervalle [0, 1] avec une erreur de 10−6 .
5. On cherche le zéro réel de f (x) = x3 + 4x2 − 10.
(a) En choisissant x0 = 1.5, appliquer la méthode du point fixe aux transformations
suivantes :
i. g1 (x) = x − x3 − 4x2 + 10
√
ii. g2 (x) = 12 10 − x3
q
10
iii. g3 (x) = 4+x
On remarquera que l’algorithme ne converge pas pour g1 , alors qu’il converge mais à
des vitesses différentes pour les autres.
(b) Approximer la racine de f avec la méthode de Newton et donner un nombre d’itéra-
tions n qui donne une précision de x∗ à 10−9 près.
6. Implémenter l’ensemble des algorithmes de ce chapitre sur SAGEmath.

ESSAI - 2017/2018 42 inestej@[Link]

Chapitre 5

Résolution de systèmes d’équations

linéaires

Un des points essentiels dans l’efficacité des méthodes envisagées concerne la taille des systèmes
à résoudre. Entre 1980 et 2000, la taille de la mémoire des ordinateurs a augmenté. La taille
des systèmes qu’on peut résoudre sur ordinateur a donc également augmenté, selon l’ordre de
grandeur suivant :

1980 Matrice pleine n = 102

Matrice creuse n = 106
2000 Matrice pleine n = 106
Matrice creuse n = 108
Le développement des méthodes de résolution de systèmes linéaires est lié a l’évolution des
machines informatiques. Un grand nombre de recherches sont d’ailleurs en cours pour profiter
au mieux de l’architecture des machines (méthodes de décomposition en sous domaines pour
profiter des architectures parallèlles, par exemple).

De nombreux efforts ont été consacrés au développement de méthodes de résolution de sys-

tèmes d’équations linéaires. Les méthodes standards incluent l’élimination de Gauss-Jordan, et
la décomposition LU. Les méthodes itératives telles que la méthode du gradient conjugué sont
généralement préférées sur les larges systèmes d’équations.

L’analyse numérique matricielle s’intéresse principalement à deux types de problèmes :

— Résoudre un système linéaire Ax = b ;
— Calculer les valeurs propres et vecteurs propres d’une matrice donnée.
Nous étudierons pour cela des algorithmes destinés à des matrices de grande taille. Comme il
s’agit de calculs volumineux sur machine, à chaque opération des erreurs systématiques dites
d’arrondi son générées. Il convient donc de se préoccuper de l’influence de ces petites erreurs
sur la solution du problème et donc d’étudier la stabilité. Nous nous intéresserons d’autre part,
au coût de ces algorithmes en termes d’opérations élémentaires.

Il s’agira donc de construire et d’analyser des méthodes numériques pour la solution de pro-
blèmes d’algèbre linéaire et d’analyse. On étudiera les notions de stabilité, de consistance, de
convergence des méthodes numériques.

43
Cours d’Analyse Numérique Ines Abdeljaoued Tej

5.1 Rappels et compléments sur les matrices

Les matrices diagonales et triangulaires jouent un rôle particulier en analyse numérique matri-
cielle pour la raison évidente que les systèmes linéaires avec de telles matrices sont de résolution
très simple. Aussi nous nous intéresserons ici aux changement de base permettant d’accéder à
des structures de ce type.

Nous considérons un espace vectoriel de dimension n sur R ou sur C. Toutes les matrices consi-
dérées dans la suite sont des matrices carrées.

Rappelons qu’une matrice est dite régulière si elle est inversible, sinon elle est singulière. Une
matrice A est dite normale si AA∗ = A∗ A où A∗ est la matrice adjointe définie par < Au, v >=<
u, A∗ v > pour tout u, v ∈ Cn . Le signe <, > dénote le produit scalaire dans Cn , soit < u, v >=
t
ūv. Une matrice unitaire est une matrice vérifiant A∗ A = AA= In . Elle est hermitienne si
A∗ = A.
Une matrice est dite symétrique si elle est réelle et si t A = A (dans ce cas, t A = A∗ ). La matrice
A est orthogonale si t AA = A t A = In .
Théorème 5.1.1. Etant donné une matrice carrée A, il existe une matrice unitaire U telle que
U −1 AU soit triangulaire. Etant donné une matrice normale A, il existe une matrice unitaire U
telle que U −1 AU soit diagonale.

5.2 Normes matricielles

Soit V un espace vectoriel de dimension n sur R ou C, on peut définir plusieurs normes
équivalentes sur V dont :
Pi=n
k x k1 = i=1 |xi |, x = t (x1 , x2 , . . . , xn ) ∈ V
k x k∞ = max
qPi=1..n |xi | √
i=n 2 t x̄x,
k x k2 = i=1 |xi | = norme de schur
Pi=n 1
p p
k x kp = ( i=1 |xi | ) .

Si V est muni d’une norme k . k, on définit sur l’ensemble des matrices carrées d’ordre n la
norme suivante :
k Ax k
k A k= supx6=0 = maxkxk=1 k Ax k
kxk
Cette quantité définit bien une norme puisqu’elle vérifie
k A k= 0 ⇒ A=0
k λA k = |λ| k A k
kA+B k ≤ k A kk B k
Nous considérons un système de n équations linéaires à n inconnues
 
a1,1 x1 +a1,2 x2 +. . . +a1,n xn = b1
 a1,1 x1 +a1,2 x2 +. . . +a1,n xn = b1 
 
 .. .. .. .. .. 
 . . . . . 
a1,1 x1 +a1,2 x2 +. . . +a1,n xn = b1

ESSAI - 2017/2018 44 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Les ai,j et les bi étant fixés.

Ce système est équivalent à Ax = b avec
   
x1 b1

a1,1 a1,2 ... a1,j ... a1,n
 x2
  b2 
 a2,1
 ..
  .. 
a2,2 ... a2,j ... a2,n  . .
   
A= =
   
 ai,1 ai,2 ... ai,j ... ai,n   xj   bi 

 .   .. 
an,1 an,2 ... an,j ... an,n  ..   . 
xn bn

5.2.1 Conditionnement d’une matrice

Un problème est bien conditionné si une petite variation des données entraîne une petite varia-
tion sur la solution, mal conditionné dans le cas contraire.

Définition 5.2.1. On appelle conditionnement d’une matrice le nombre le nombre

c(A) =|| A || . || A−1 ||

On se propose d’examiner la sensibilité d’une méthode de résolution :

1. aux erreurs d’arrondis
2. aux ”petites perturbations” des données du problème.
On mesurera cette sensibilité par un nombre appelé conditionnement, qui dépend d’une norme
sur Kn . Lorsqu’il est très élevé, cela peut mettre cause la pertinence du résultat numérique
obtenu.
On part d’un système Ax = b, de taille n×n, et on se donne par ailleurs ∆A et ∆b, présumés
petits par rapport à A et b et tel que A + ∆A = A soit aussi inversible. On se propose de
comparer les solutions x et x des systèmes : Ax = b et le système perturbé Ax = b. Soient
le système Ax = b et le système perturbé Āx̄ = b̄. Les erreurs relatives sur A et b sont :
A = ||Ā−A||
||A||
, b = ||b̄−b||
||b||
.
||x̄−x|| c(A)
Théorème 5.2.1. Si A c(A) < 1 : ||x||
≤ (
1−A c(A) A
+ b ).
Le théorème 5.2.1 donne une borne de l’erreur relative sur x : la solution du système Ax = b
risque d’être très sensible aux variations sur les données si le conditionnement est élevé.
Un problème est bien conditionné si une petite variation des données entraîne une petite varia-
tion sur la solution, mal conditionné dans le cas contraire.
Exemple
 5.2.1. On souhaite résoudre
 le  linéaire Ax = b, où A et b sont donnés
système   par :
10 7 8 7 32 1
 7 5 6 5   et b =  23 . La solution du système est : x =  1 .
   
A=  8 6 10 9   33   1 
7 5 9 10 31 1
   
32.1 9.2
 22.9   −12.6 
On résout le système perturbé Ax = b avec b =   33.1 . On trouve x =  4.5  .
  

30.9 −11

ESSAI - 2017/2018 45 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Autrement dit, de très petites variations sur b ont conduit à de grandes variations sur x. Le
conditionnement de A est c(A)=4488 >>> 1. Ce phénomène de mauvais conditionnement
explique pour partie la difficulté de prévoir certains phénomènes. Les appareils de mesure ne
sont jamais parfaits, et il est impossible de connaître exactement b. Cela peut entrainer une
très grande imprécision sur la valeur de x.
Les conditionnements de matrice utilisés dans la pratique correspond aux trois normes usuelles
[Link] , p=1, 2, 8. On note c (A) = kAkp ∗ kA−1 kp .
Si le conditionnement de A est de l’ordre 1 alors on dit que A est bien conditionnée sinon, elle
est mal conditionnée.

5.2.2 La méthode de Gauss

Exercice 1. On considère le système Ax = b suivant :
   
3 −1 0 5
 −2 1 1 x =  0 
2 −1 4 15

3 -1 0 5 3 -1 0 5
-2 1 1 0 l2 ⇒ l2 + 2/3l1 ↔ 0 1/3 1 10/3
2 -1 4 15 l3 ⇒ l3 − 2/3l1 0 -1/3 4 35/3 l3 ⇒ l3 + l2
Nous effectuons une élimination ascendante (back-substitution) ou une remontée afin de ré-
soudre le nouveau système triangulaire obtenu par la méthode de Gauss :
3 -1 0 5 x3 = 3
↔ 0 1/3 1 10/3 x2 = 1
0 0 5 5 x1 = 2
Nous considérons un système de n équations linéaires à n inconnues
a1,1 x1 + a1,2 + · · · + a1,j xj + · · · + a1,n xn
 
= b1
 a2,1 x1 + a2,2 + · · · + a2,j xj + · · · + a2,n xn = b2 
 . .. 
 .
 . . 

A=
 ai,1 x1 + ai,2 x2 + · · · + ai,j xj + · · · + ai,n xn = bi 

 . .. 
 .. . 
an,1 x1 + an,2 x2 + · · · + an,j xj + · · · + an,n xn = bn
les ai,j et les bi,j étant fixés. Ce système est équivalent à Ax = b et par combinaisons linéaires
successives des lignes du système, on triangularise supérieurement. La résolution du système
(1)
A(1) x = b(1) avec A(1) = A une matrice inversible telle que a1,1 6= 0 :
 
(1) (1) (1) (1)
a1,1 a1,2 . . . a1,j . . . a1,n  (1) 
 (1)
 a2,1 a(1) (1) (1)  b1
2,2 . . . a2,j . . . a2,n 
 b(1)
 .  
 . .. .. .. .. ..  2 
 . . . . . . 

(1) (1)
 (1) 
A =  (1) et b =  b3 

.
 ai,1 a(1) . . . a
(1)
. . . a
(1) 
.
 .. 
i,2 i,j i,n
  
 . .. .. .. .. .. 
 .
 . . . . . . 

(1)
(1) (1) (1) (1) bn
an,1 an,2 . . . an,j . . . an,n

ESSAI - 2017/2018 46 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Annuler l’élément en position (i, 1) pour i = 1, . . . , n revient à soustraire de la ième ligne, la

(1) (1)
première ligne multipliée par mi,1 = aj,1 /a1,1 .
Cet ensemble d’éliminations revient à prémultiplier la matrice A(1) par une matrice M1 vérifiant
 
1 0 0 ··· 0
 -m1 1 0 · · · 0 
 
M1 =  -m2 0 1 · · · 0 
 
 .. .. .. .. .. 
 . . . . . 
-mn 0 0 · · · 1

a1,1 est appelé élément pivot de la transformation M1 et on obtient un nouveau système équi-
valent au système de départ : M1 A(1) x = M1 b(1) où
 
(1) (1) (1) (1)
a1,1 a1,2 . . . a1,j . . . a1,n  (1) 
 (2) (2) (2)  b1
 0 a2,2 . . . a2,j . . . a2,n   (2) 
 .
 . .. .. .. .. ..   b2 
(2)  . . . . . .  (2) (1)
 (2) 
A = et b = M1 b =   b3 
(2) (2)

(2) 
 0 a . . . a . . . a .
 .. 
i,2 i,j i,n
  
 . .. .. .. .. .. 
 .
 . . . . . . 

(2)
(2) (2) (2) bn
0 an,2 . . . an,j . . . an,n

(2) (2)
les seuls éléments à calculer sont les ai,j , pour i = 2..n et bi , i = 2..n. Ce calcul s’effectue
grâce aux formules suivantes :

(1) (1)
 mi,1 = ai,1 /a1,1

(2) (1)
bi = mi,1 bi
 a(2) = a(1) − m a(1)

i,j i,j i,1 i,j

(2)
De la même manière, on réitère la procédure de l’élimination de Gauss en supposant que a2,2 6= 0
selon les formules 
(2) (2)
 mi,2 = ai,2 /a2,2

(3) (2)
bi = mi,2 bi
 a(3) = a(2) − m a(2)

i,j i,j i,2 i,j

par application de n − 1 transformations de même type de M1 et M2 on obtient finalement,

pour autant que les pivots successifs soient non nuls, le système : Mn−1 Mn−2 . . . M2 M1 A(1) x =
Mn−1 Mn−2 . . . M2 M1 b(1) équivalent au système A(1) x = b(1) , mais dont la matrice est triangu-
laire supérieure.

Complexité de calcul
Le coût de la méthode de Gauss est donné par : n−1
P Pn−1
k=1 (2(n−k)(n−k +1)+n−k) = i=1 (2i(i+
1) + i) opérations soit 23 n3 + 12 n2 − 76 n auxquelles il faut encore ajouter n 2 opérations pour
l’élimination ascendante (algorithme de remontée). La méthode de Gauss exige donc finalement
le calcul de 32 n3 + 32 n2 − 76 n opérations en tout. Tout ceci n’est évidemment valable que si les
pivots successifs des transformations sont non nuls.

ESSAI - 2017/2018 47 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Remarque 5.2.1. Si, au cours des calculs, on rencontre un pivot nul, il faut permuter la ligne de
ce pivot avec une ligne d’indice supérieur, présentant un pivot non nul dans la même colonne.
Ceci sera toujours possible si le système donné est compatible (detA 6= 0). En pratique, on ne
permute pas les lignes entre-elles, on change simplement le numéro des lignes à permuter en
stockant dans un vecteur de travail les indices des pivots successifs. Pour la "back-substitution"
on tient compte de l’ordre utilisé lors de la triangularisation.

5.2.3 La factorisation LU d’une matrice

On dit qu’une matrice A(n, n) possède une factorisation LU si A = LU avec L une matrice
(n, n) triangulaire inférieure à diagonale unité et U une matrice triangulaire supérieure.

Proposition 5.2.2 (Unicité). On suppose que A est régulière et admet une factorisation A=
LU, alors U est régulière et la factorisation est unique.

Preuve : det(A)=det(LU)=det(L)det(U)=det(U).
Si A est régulière (det(A) 6= 0) alors U l’est aussi. Supposons que A=L1 U1 =L2 U2 alors L−1
2 L1 =
−1 −1
U2 U1 = D. Comme L2 est une matrice triangulaire inférieure à diagonale unité donc L2 trian-
gulaire inférieur à diagonale unité. De même, comme L1 est une matrice triangulaire inférieure
à diagonale unité donc D = L−1 2 L1 est une matrice triangulaire inférieure à diagonale unité.
D’autre part, U1 est une matrice triangulaire supérieure donc U1−1 est une matrice triangulaire
supérieure et U2 est une matrice triangulaire supérieure donc D = U2 U1−1 est une matrice
triangulaire supérieure. Ainsi D est la matrice unité, L1 = L2 et U1 = U2 .

Proposition 5.2.3 (Existence). On suppose que l’élimination de Gauss est faisable sans per-
mutation de lignes sur un système de matrice A supposée régulière. Alors A possède une facto-
risation LU avec U régulière.

La méthode de Gauss fournit un système à matrice triangulaire supérieure :

Mn−1 Mn−2 . . .M2 M1 Ax = Mn−1 Mn−2 . . .M2 M1 b = b 1 . Posons Mn−1 Mn−2 . . .M2 M1 A = U.
Comme les Mi sont triangulaires inférieures à diagonale unité le produit
Mn−1 Mn−2 . . .M2 M1 est aussi triangulaire inférieure à diagonale unité et l’on a :
A = (Mn−1 Mn−2 . . .M2 M1 )−−1 U = LU où L = (Mn−1 Mn−2 . . .M2 M1 )−−1 est une matrice
triangulaire inférieure à diagonale unitaire.

Théorème 5.2.4. La matrice A(n,n) posséde une factorisation LU avec U régulière si et seule-
ment si ses sous-matrices principales Ak k=1. . .n sont régulières.

Construction de A : Si on connaît une factorisation (on dit aussi décomposition) LU de la

matrice A du système Ax = b , sa solution s’obtient de la façon suivante : 1◦ ) on calcule la
solution du système Ly = b 2◦ ) on résout ensuite Ux = y . En effet : Ax = LUx = Ly =
b .
La connaissance de la décomposition LU permet de résoudre le système de départ au prix de 2n2
opérations arithmétiques (supplémentaires) puisqu’il faut résoudre deux systèmes triangulaires
à n inconnues.
Remarque 5.2.2. Le coût de la factorisation LU est du même ordre que celui de la méthode de
Gauss soit 2n2 /3.

ESSAI - 2017/2018 48 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

La méthode ne fonctionne qu’avec les matrices carrés régulières (det A 6= 0). La méthode fournit
la matrice inverse de A grâce à la relation A−1 = U−1 L−1 . La matrice U est identique à la ma-
trice triangulaire fournie par la méthode de Gauss. La matrice L = (Mn−1 Mn−2 . . .M2 M1 )−−1
s’écrit: 
1 0 0 ··· 0
 m21
 1 0 ··· 0 
L=
 m31 m32 1 · · · 0  où les Mik sont les facteurs multiplicatifs rencontrés lors

 .. .. .. . . .. 
 . . . . . 
mn1 mn2 mn3 · · · 1
de la triangularisation de Gauss.
 
4 −9 2
Exemple 5.2.2. Effectuer la factorisation LU de la matrice : A =  2 −4 4 .
−1  2 2 
4 −9 2
D’après les formules : m21 = 0.5 et m31 = -0.25 ce qui nous donne A =  0 0.5 3 .
 0 −0.25  2.5
1 0 0
De la même manière m32 = -0.5 et nous obtenons que L =  0.5 1 0  et U =
  −0.25 −0.5 1
4 −9 2
 0 0.5 3 .
0 0 4

5.2.4 La factorisation de Cholesky

Cette méthode s’applique uniquement aux matrices symétriques définies positifs et permet
comme déjà mentionné la factorisation A= RRT où R est une matrice triangulaire inférieure.

Définition 5.2.2. Soit A une matrice

 symétrique, on dit que A est définie positive si :
n n
 T
∀ x ∈ IR x A x > 0 ∀ x ∈ IR xT A x > 0
T T
x A x = xA x = 0 ⇔ x = 0 x A x = xA xT = 0 ⇔ x = 0
T

Théorème 5.2.5. Une matrice A symétrique est définie positive si tous ses déterminants sont
positifs

Théorème 5.2.6. Si les valeurs propres d’une matrice A symétrique sont toutes strictement
positives alors A est définie positive. De même si A est définie positive alors toutes ses valeurs
propres sont strictement positives.

Rappel : Les valeurs propres d’une matrice triangulaire sont les éléments de sa diagonale.

Théorème 5.2.7. Si A est une matrice symétrique définie positive, alors il existe au moins
une matrice réelle triangulaire inférieure R tel que A=RRT
de plus si Rii > 0 ∀ i=1,2,. . .,n alors la factorisation est unique.

Construction de R par la méthode de Cholesky :

√
r11 = a11
ai1
ri1 = r11
pour i = 2 . . . n

ESSAI - 2017/2018 49 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

q
rii = aii − i−1 2
P
k=1 rik pour i = 2 . . . n
h i
rji = r1ii aji − i−1
P
r
k=1 jk ∗ rik pour i = 2 . . . n , j = i + 1 . . . n
 
1 1 1 1
 1 5 5 5 
Exemple 5.2.3. Soit A =   1 5 14 14 . La matrice symétrique A décomposée par la mé-


1 5 14 15
 
1 0 0 0
 1 2 0 0 
thode de Cholesky tel que A=RRT avec R =   1 2 3 0 


1 2 3 1
Remarque 5.2.3. Le coût du calcul est de l’ordre de n2 /2 opérations arithmétiques et n ex-
tractions de racine carré, c’est plus économique que la factorisation LU pour une matrice
symètrique.
Cette méthode fournit un test efficace pour déterminer si une matrice symétrique est définie
positive ou non :
Dans de très nombreux problèmes concrets, les systèmes à résoudre ont des matrices symétriques
définies positives en raison même des phénomènes qu’ils modélisent. La méthode de factorisation
de Cholesky revêt donc une importance particulière.
Toutefois, comme toutes les méthodes directes, son coût de calcul et les erreurs d’arrondis
deviennent un handicap quasiment insurmontable pour des systèmes comptant de l’ordre de
200 équations. En pratique on est surtout amené à résoudre des systèmes de plusieurs dizaines
de milliers d’équations dont les matrices, fort heureusement, sont généralement creuses. La
rèsolution de tels systèmes nécessite l’emploi de méthodes itératives.

5.2.5 Valeurs propres et vecteurs propres

Aujourd’hui, le calcul des valeurs et vecteurs propres est indispensable dans toutes les branches
de la science, en particulier pour la solution des systèmes des équations linéaires, en théorie
de stabilité, pour les questions de convergence de processus itératifs, et en physique et chimie
(mécanique, circuits, cinétique chimique, équation de Schrödinger).
Il existe plusieurs méthodes pour la détermination des valeurs et vecteurs propres d’une matrice :
1. Méthode de la puissance
2. Méthode de la puissance inverse de Wielandt
3. Transformation sous forme de Hessenberg (ou tridiagonale)
4. L’itération orthogonale
5. L’algorithme QR

5.2.6 La méthode de la puissance

Définition 5.2.3. Une matrice A admet une valeur propre Dominante lorsqu’elle est la plus
grande valeur propre en module :
|λ1 | ≥ |λ2 | ≥ · · · ≥ |λn |
Le vecteur propre associé est dit dominant.

ESSAI - 2017/2018 50 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Remarque 5.2.4. La méthode de la puissance permet de calculer la valeur et le vecteur propre

dominants de A.

Théorème 5.2.8. Supposons que λ et ν sont respectivement la valeur propre et le vecteur

propre d’une matrice A. Pour toute constante α, on a (λ - α) est la valeur propre de (A - αIn )
et de vecteur propre A

Preuve : Immédiate.

Théorème 5.2.9. si λ et ν sont respectivement la valeur propre et le vecteur propre d’une

matrice A inversible, alors si λ ? 0, on a 1/λ et ν sont respectivement la valeur propre et le
vecteur propre de A−1 .

Preuve : On a Av=λv, λ A−1 v = λ A−1 λv = A−1 Av = v et donc A−1 v = 1/λ v.

Définition 5.2.4. . Un vecteur propre est dit normalisé si la plus grande de ses coordonnées
est égale à 1.
 
v1
 v2 
Remarque 5.2.5. Il est facile de transformer un vecteur propre v =  ..  en formant un
 
 . 
vn
 
v1
 v2 
nouveau vecteur v 0 = 1c  ..  avec c = vj et
 
 . 
vn

|vj | = max |vi | .

i=1..n

Théorème 5.2.10 (Théorème de la méthode de puissance). Supposons que la matrice A ad-

mette une valeur propre dominante et un vecteur propre normalisé v. Ces deux quantités peuvent
être déterminées à partir laméthode
 de la puissance.
1
 1 
On commence par X0 =  ..  ; il est important que X0 soit différent des autres vecteurs
 
 . 
1
propres.
X0 = ni=1 αi vi =αv+ ni=2 αi vi où v est le vecteur normalisé
P P
On génère alors une suite (Xk ) de manière récursive en utilisant la définition suivante :

Yk = AX k
Xk+1 =Ck+1 Yk

où Ck+1 est la plus grande des coordonnées en valeur absolue de Yk

alors la suite (Xk ) converge vers v (vecteur normalisé) et la suite Ck converge vers (valeur
propre dominante). 
limk→+∞ Xk =v
limk→+∞ Ck =

ESSAI - 2017/2018 51 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Preuve : Sachant que A admet n valeurs propres associées à n vecteurs propres distincts
vj 1 ≤ j ≤ n, linéairement indépendants normalisées et formant une base de dimension n.
On peut alors écrire X0 dans cette base :

X0 =α1 v1 +α2 v2 + · · · +αn vn , α1 6= 0

 
x1
 x2 
on suppose que pour X0 =  ..  , maxi=1..n |xi | = 1
 
 . 
xn

Y0 = AX0 = A [α1 v1 +α2 v2 + · · · +αn vn ]

= α1 Av 1 +α2 Av 2 + · · · +αn Avn

= α1 λ1 v 1 +α2 λ2 v 2 + · · · +αn λn vn
 
λ1 λ2 λn
= λ1 α1 v +α2 v + · · · +αn vn
λ1 1 λ1 2 λ1
     
λ2 λn
= λ1 α1 v1 +α2 v + · · · +αn vn
λ1 2 λ1
en normalisant, on obtient :
     
λ1 λ2 λn
X1 = α1 v1 +α2 v + · · · +αn vn
C1 λ1 2 λ1

après k itérations
Yk−1 = AXk−1
" "  k−1  k−1 ##
k−1
(λ1 ) λ2 λn
=A α1 v1 +α2 v + · · · +αn vn
C1 C2 . . .Ck−1 λ1 2 λ1
" "  k−1  k−1 ##
(λ1 )k−1 λ2 λn
= α1 Av1 +α2 Av + · · · +αn Av n
C1 C2 . . .Ck−1 λ1 2 λ1
" "  k−1  k−1 ##
(λ1 )k−1 λ2 λn
= α1 λ1 v1 +α2 v + · · · +αn nvn
C1 C2 . . .Ck−1 λ1 2 2 λ1
" "  k  k ##
(λ1 )k λ2 λn
= α1 v1 +α2 v + · · · +αn vn
C1 C2 . . .Ck−1 λ1 2 λ1
 k
λj λ
or λ1 <1 ∀ j ∈ [2, n] alors limk→+∞ αj λ1j =0 ∀ j ∈ [2, n]
ainsi
α1 (1 )k
lim Xk = v1
k→+∞ C1 C2 . . .Ck
comme Xk et v1 sont normalisés, on a :
)k
limk→+∞ C1αC1 (21...C k
= 1 (éq 3.3.1)

ESSAI - 2017/2018 52 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

donc
lim Xk =v1
k→+∞

on écrit l’équation 3.3.1 pour l’étape (k-1), on aura :

)k−1
limk→+∞ C1αC1 (21...C k−1
= 1 (eq 3.3.2)
en effectuant la division de l’eq 3.3.1 par l’eq 3.3.2, on aura :
1
lim =1
k→+∞ Ck
lim Ck =1
k→+∞

ESSAI - 2017/2018 53 inestej@[Link]

Chapitre 6

Equations différentielles ordinaires

6.1 Rappel
On dit qu’une fonction g définie sur un intervalle [a, b] est lipschitzienne de constante L si

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

En pratique, la fonction g est de classe C 1 sur [a, b] et cette inégalité est une conséquence de la
formule des accroissements finis :

g(x) − g(x0 ) = g 0 (c)(x − x0 ) c ∈ (x, x0 )

et on pose L = maxx∈[a,b] | g 0 (x) |.

Lemme 6.1.1. 1. On montre par récurrence que si un+1 ≤ un α + β pour tout n alors
n −1
un ≤ u0 α + β αα−1
n
.
2. On montre grâce à une étude de fonction que si X > 0 alors (1 + X)n ≤ enX .

6.2 Problème de Cauchy

Soient f une fonction de [a, b] × R dans R continue et τ ∈ R. Le problème de Cauchy ou
problème de la condition initiale est le suivant : existe-t-il une fonction y(x) définie sur [a, b] à
valeurs dans R dérivable sur [a, b] et vérifiant :

y(a) = τ
(6.2.1)
y(x) = f (x, y(x))

Théorème 6.2.1 (admis). On supprose que f est lipschitzienne par rapport à sa deuxième
variable uniformément par rapport à sa première variable. Il existe donc L > 0 tel que

| f (x, z1 ) − f (x, z2 ) |≤ L | z1 − z2 | , ∀x ∈ [a, b] , z1 , z2 ∈ R. (6.2.2)

Alors le probmème de Cauchy admet une solution unique.

La condition (6.2.2) sera satisfaite si f admet une dérivée partielle par rapport à sa deuxième
variable qui soit bornée dans [a, b] × R. On obtient grâce au Théorème des Accroissements

54
Cours d’Analyse Numérique Ines Abdeljaoued Tej

Finis :
∂f
f (x, z1 ) − f (x, z2 ) = (x, zc )(z1 − z2 )
∂z
∂f
⇒ | f (x, z1 ) − f (x, z2 ) | ≤ maxx∈[a,b] | (x, zc ) | | z1 − z2 | .
∂z

On présente dans les sections suivantes des méthodes d’approximation de la solution y(x) de
(6.2.1). Pour cela, on subdivise l’intervalle [a, b] uniformément et on approche la fonction y en
ces points.
Soit n un entier naturel non nul. On pose xi = a + ih avec h = b−a n
et i = 0..n. On va calculer
une approximation yi de y(xi ) pour i = 0..n et estimer l’erreur d’approximation ei = yi − y(xi ).

6.3 Méthode d’Euler

Partant de x0 = a, on connaît y0 = y(x0 ) = τ . L’idée de la méthode d’Euler est de considérer que
la courbe y sur [x0 , x1 ] n’est pas très éloignée de sa tangente en x0 : une estimation rationnelle de
y(x1 ) est y1 = y0 + hf (x0 , y0 ). En effet, pour x1 , l’équation de la tangente est y 0 (x0 )(x − x0 ) + y0 .
On estime alors que y 0 (x1 ) n’est pas très différent de f (x1 , y1 ). Sur [x1 , x2 ], on remplace la courbe
par sa tangente approchée en x1 ce qui donne l’estimation y2 de y(x2 ) : y2 = y1 + hf (x1 , y1 ) et
on poursuit de la même manière.
Algorithme d’Euler :
Entrée : n, τ , a et b, et la fonction f .
Sortie : Une approximation de yx.
x0 := a ;
y0 := τ ;
h := b−an
;
Pour i allant de 0 à n Faire
xi+1 := xi + h ;
yi+1 := yi + hf (xi , yi ) ;
Fin Pour ;
Imprimer(La méthode d’Euler donne une approximation de y en (xi , yi ), i = 0..n).
Fin.
Théorème 6.3.1. Si f est continue sur [a, b] × R, lipshitzienne par rapport à sa deuxième
variable uniformément par rapport à sa première variable et si la soltion y de l’EDO est de
Classe C 2 sur [a, b], alors :

(eL(b−a) − 1) maxx∈[a,b] | y 00 (x) |

| en |≤ h .
L 2
Preuve : Voir la démonstration du théorème 6.5.1.
Définition 6.3.1. Une méthode de calcul numérique qui fournit une erreur | en |≤ khp est dite
d’ordre p. La méthode d’Euler est d’ordre 1.
Une méthode d’ordre 1 ne converge pas assez vite pour donner des résultats pratiques intéres-
sants. On va étudier dans les deux sections suivantes deux méthodes d’ordre plus élevé.

ESSAI - 2017/2018 55 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

6.4 Méthode de Taylor d’ordre 2

L’idée consiste à remplacer la courbe y sur un intervalle [xi , xi+1 ] par une parabole (et non pas
une droite comme dans la méthode d’Euler).
En utilisant la formule de Taylor sur y au voisinage de xn , on obtient : y(xn+1 ) = y(xn ) +
2 3
hy 0 (xn ) + h2 y 00 (xn ) + h6 y (3) (cn ) avec cn ∈ [xn , xn+1 ]. Comme y 0 (x) = f (x, y(x)), on a : y 00 (x) =
2
( ∂f
∂x
+ f ∂f
∂z
)(x, y(x)). On pose yn+1 = yn + hf (xn , yn ) + h2 ( ∂f ∂x
+ f ∂f
∂z
)(xn , yn ).
Algorithme de Taylor d’ordre 2 :

Entrée : n, τ , a et b, et la fonction f .
Sortie : Une approximation de yx.
x0 := a ;
y0 := τ ;
h := b−an
;
Pour i allant de 0 à n Faire
xi+1 := xi + h ;
2
yi+1 := yi + hf (xi , yi ) + h2 ( ∂f
∂x
+ f ∂f
∂z
)(xi , yi ) ;
Fin Pour ;
Imprimer(La méthode de Taylor d’ordre 2 donne une
approximation de y en (xi , yi ), i = 0..n).
Fin.

Théorème 6.4.1. Si f est de Classe C 2 sur [a, b] × R, lipshitzienne par rapport à sa deuxième
variable uniformément par rapport à sa première variable (avec la constante L) et si ∂f
∂x
+ f ∂f
∂z
est lipshitzienne par rapport à sa deuxième variable uniformément par rapport à sa première
variable (avec la constante L1 ), alors pour h ∈ [0, h0 ] on a :
h0
(e(b−a)(L+ 2 L1 )
− 1) maxx∈[a,b] | y (3) (x) | 2
| en |≤ h .
L 6
Preuve : Voir la démonstration du théorème 6.5.1.

6.5 Méthode de Runge-Kutta

Définition 6.5.1. Une méthode générale à un pas est définie par yn+1 = yn + h.φ(xn yn , h) où φ
est continue sur [a, b]×R×[0, h0 ] lipshitzienne par rapport à sa seuxième variable uniformément
par rapport à la première et à la troisième de constante Lφ .

Théorème 6.5.1. Soit une méthode à un pas définie par yn+1 = yn + h.φ(xn yn , h) où φ est
continue sur [a, b] × R × [0, h0 ] lipshitzienne par rapport à sa seuxième variable uniformément
par rapport à la première et à la troisième de constante Lφ et telle que

y(x + h) − y(x)
− φ(x, y(x), h) ≤ khp
h
où y est solution de l’EDO y 0 (x) = f (x, y(x)) et f de Classe C p sur [a, b] × R. Alors

| en | ≤ K hp .

ESSAI - 2017/2018 56 inestej@[Link]

Chapitre 7

Application aux réseaux de neurones

Un réseau de neurones artificiel est une structure composée d’entités (les neurones) capables de
calcul interagissant entre eux. Il permet de traiter, par le biais de l’informatique, des problèmes
de différentes natures que les outils classiques ont du mal à résoudre. Les explications et autres
définitions sont tirés du polycope de Marie Cottrel de l’Université de la Sorbonne 1 .

7.1 Historique
Les recherches sur le cerveau ont fait d’énormes progrès, tant du point de vue de la compréhen-
sion du fonctionnement qu’à l’échelle cellulaire et membranaire. On sait désormais mesurer et
enregistrer l’activité électrique d’un neurone ou d’un ensemble de neurones. Le cerveau humain
dispose de l’ordre de 1011 avec plus de 1015 connexions. Il y a des propriétés collectives des
neurones qui agissent en relations excitatrices ou inhibitrices. On distingue dans le cerveau des
sous-ensembles, des connexions internes et externes, des entrées et des sorties, pour réaliser
certaines tâches (parole, vision, émotion,...).

Les années 1940 sont connues pour l’apprentissage du codage et de la recherche opérationnelle :
A. Turing avait mis en place la machine de Turing, premier ordinateur conçu pour déchiffrer
le code Enigma. Depuis, le besoin de proposer des machines ou des modèles reproduisant des
fonctions des cerveaux humains ont été développés. Il s’agit d’outils de calcul, de mémoire
pour reproduire des comportements intelligents. Des chercheurs dans différents domaines ont
travaillé sur l’intelligence artificielle.
Le tableau suivant donne brièvement un historique des réseaux de neurones artificiels :
1940 - 1960 Mac Cullogh & Pitts 43 modèle de neurone
concepts Hebb loi d’adaptation
Rosenblatt 58 perceptron
1960 - 1980 Widrow-Hoff 58 adaline
transition et déclin Minsky & Papert 69 limites aux perceptron
1980 - ... Hopfield 82 réseaux dynamiques
renouveau Kohonen 72, 82, 84 auto-organisation
Rumelhart, Le Cun 86 rétro-propagation
Werbos 74
1. Marie Cottrel, Les réseaux de neurones historique, méthodes et applications SAMOS-MATISSE, Université
Paris 1, Sorbonne

57
Cours d’Analyse Numérique Ines Abdeljaoued Tej

7.2 Définitions
Un neurone biologique est une cellule nerveuse dont la fonction est de transmettre, dans cer-
taines conditions, un signal électrique (0.01s). C’est comme un relai entre une couche de neu-
rones et celle qui la suit. Le corps

Le corps d’un neurone est relié d’une part par un ensemble de dendrites (entrées du neurone)
et d’autre par à un axone, qui est la partie étirée de la cellule. L’axone représente pour nous la
sortie du neurone. Les charges électriques, reçues par les connexions synaptiques (en amont),
s’accumulent dans le noyau du neurone jusqu’à atteindre un certain seuil. Au delà de ce seuil,
la transmission du signal électrique se déclenche via son axone vers d’autres neurones (en aval).
On remarque que les liaisons axone-dendrite entre deux neurones (connexions synaptiques) ne
sont pas toutes de la même efficacité. Ainsi l’entrée associée à une certaine dendrite du neurone
pourra avoir plus d’influence qu’une autre sur la valeur de sortie. On peut représenter la qualité
de la liaison par un poids, sorte de coefficient s’appliquant au signal d’entrée. Le poids sera
d’autant plus grand que la liaison est bonne. Un poids négatif aura tendance à inhiber une
entrée, tandis qu’un poids positif viendra l’accentuer 2 .

7.2.1 Neurone formel

Le neurone formel reproduit le neurone biologique. Il doit être capable de recevoir, en entrée,
différentes informations provenant des neurones environnants. Il doit analyser et ajuster ces
informations, de manière à envoyer en sortie une réponse. On assimile donc à un neurone un
quadruplet (poids w1 , wn , biais w0 , fonction d’activation f , seuil θ) :

1 si w0 + w1 x1 + · · · + wn xn > θ
f (x1 , . . . , xn ) =
−1 sinon.

On utilise aussi une notation plus simple :


~ · ~x > θ
1 si w
f (~x) =
−1 sinon.

7.2.2 Fonctions d’activation

Il y a plusieurs exemples de fonctions d’activation utilisées pour les réseaux de neurones. Elle
ressemblent le plus souvent à des fonctions Tout ou rien :

2. Texte et dessin d’Adrien Vergé, Réseaux de neurones artificels, Lycée Michelet, TIPE 2009.

ESSAI - 2017/2018 58 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Nom Valeur Représentation

Seuil f (x) = 0 si x ≤ 0 et f (x) = 1 sinon

Seuil symétrique f (x) = −1 si x < 0

f (0) = 0
f (x) = 1 si x > 0

Linéaire saturée f (x) = 0 si x < 0

f (x) = x pour 0 ≤ x ≤ 1
f (x) = 1 si x > 0

Linéaire f (x) = x

1
Sigmoïde f (x) =
1 + e−x

ex − e−x
Tangente hyperbolique f (x) =
ex + e−x

7.2.3 Perceptron multi-couches

On souhaite classer des images de 4x6 pixels dans deux catégories, selon des critères dépendant
de la couleur et de la disposition de ces pixels. La première couche de notre réseau devra com-
porter 24 entrées, une entrée pour la valeur de chaque pixel. La couche de sortie fournira deux
valeurs selon la catégorie. Le classement des images se fera selon les «affinités» : on placera
l’image dans la catégorie avec laquelle l’affinité est grande.

Chaque neurone comportant plusieurs entrée et une seule sortie, on peut d’une façon générale,
superposer dans une même couche, la couche 1 qui a s1 neurones formels à e1 entrées. On
obtient alors un système admettant e1 entrées et s1 sorties. Puis, lorsqu’on concatène N couches
connectées les unes à la suite des autres (le nombre d’entrées d’une couche étant égal au nombre
de sorties de la couche précédente : sn−1 = en ), on obtient un réseau de N couches de neurones.

7.3 Fonctionnement du réseau de neurones

Les applications des réseaux neuronaux vont de la reconnaissance vocale, en passant par la
classification d’image ou de la prédiction financière. Les réseaux de neurones sont bien indiqués
lorsqu’il s’agit de modéliser des structures complexes et des données irrégulières. En particulier
quand les relations entre les variables sont non linéaires. Ils requièrent un pré-traitement sur

ESSAI - 2017/2018 59 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

les données aberrantes auxquelles ils sont sensibles. Pour les PMC, les données doivent être
normalisées. La particularité du réseau de neurones comme outil mathématique repose sur sa
capacité d’apprentissage. Grâce à un processus d’entraînement, on peut améliorer les perfor-
mances du réseau, c’est-à-dire qu’en réponse à une certaine entrée, il fournit la bonne sortie. Il
s’agira de modifier les poids des liaisons en réaction aux entrées qu’on lui soumet, de manière à
ce que la sortie corresponde au résultat désiré. Le réseau s’améliore en réduisant l’erreur qu’il
commet à chaque simulation. On utilise pour cela un algorithme d’apprentissage. Il y a deux
types d’algorithme d’apprentissage :
— Les apprentissages supervisés : on fournit un grand nombre de couples (x, c) où x est
l’entrée et c est la sortie désirée. La correction s’effectue selon l’erreur obtenue pour
chaque couple. Dans le cas d’un perceptron, la correction de l’erreur se fait comme suit :

wi ← wi + ∆wi

avec
∆wi = η(t − c)xi
où
t = f (~x) est la valeur cible
c est la sortie désirée du perceptron
η une constante de petite taille (e.g., .1) appelée taux d’apprentissage
On prouve la convergence lorsque les données d’apprentissage sont linéairement sépa-
rables et η suffisamment petit. On parle en général d’apprentissage par l’exemple.
— Les apprentissages non supervisés : ne connaissant pas la sortie désirée pour une entrée
donnée, on n’a pas de signal d’erreur le réseau apprendra par lui-même à classer des
entrées similaires en trouvant des points communs aux stimuli appliquées.
La correction de l’erreur en réaction à des stimuli en entrée constitue le processus d’apprentis-
sage
To understand, consider simpler linear unit, where

o = w0 + w1 x1 + · · · + wn xn
Let’s learn wi ’s that minimize the squared error
1X
~ ≡
E[w] (td − od )2
2 d∈D
Where D is set of training examples.
 
∂E ∂E ∂E
∇E[w]
~ ≡ , ,···
∂w0 ∂w1 ∂wn

Training rule :
~ = −η∇E[w]
∆w ~
i.e.,
∂E
∆wi = −η
∂wi

ESSAI - 2017/2018 60 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

∂E ∂ 1X
= (td − od )2
∂wi ∂wi 2 d
1X ∂
= (td − od )2
2 d ∂wi
1X ∂
= 2(td − od ) (td − od )
2 d ∂wi
X ∂
= (td − od ) (td − w~ · x~d )
d
∂w i

∂E X
= (td − od )(−xi,d )
∂wi d

Gradient-Descent(training_examples, η)
Each training example is a pair of the form h~x, ti, where ~x is the vector of input
values, and t is the target output value. η is the learning rate (e.g., .05).
— Initialize each wi to some small random value
— Until the termination condition is met, Do
— Initialize each ∆wi to zero.
— For each h~x, ti in training_examples, Do
— Input the instance ~x to the unit and compute the output o
— For each linear unit weight wi , Do

∆wi ← ∆wi + η(t − o)xi

— For each linear unit weight wi , Do

wi ← wi + ∆wi

σ(x) is the sigmoid function

1
1 + e−x
dσ(x)
Nice property : dx
= σ(x)(1 − σ(x))
L’erreur :
∂E ∂ 1X
= (td − od )2
∂wi ∂wi 2 d∈D
1X ∂
= (td − od )2
2 d ∂wi
1X ∂
= 2(td − od ) (td − od )
2 d ∂wi
 
X ∂od
= (td − od ) −
d
∂wi
X ∂od ∂netd
= − (td − od )
d
∂netd ∂wi

ESSAI - 2017/2018 61 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

But we know :
∂od ∂σ(netd )
= = od (1 − od )
∂netd ∂netd
∂netd ~ · ~xd )
∂(w
= = xi,d
∂wi ∂wi
So :
∂E X
= − (td − od )od (1 − od )xi,d
∂wi d∈D

Backpropagation Algorithm :
Initialize all weights to small random numbers.
Until satisfied, Do
— For each training example, Do
1. Input the training example to the network and compute the network outputs
2. For each output unit k

δk ← ok (1 − ok )(tk − ok )
3. For each hidden unit h
X
δh ← oh (1 − oh ) wh,k δk
k∈outputs

4. Update each network weight wi,j

wi,j ← wi,j + ∆wi,j

where

∆wi,j = ηδj xi,j

La faiblesse des réseaux de neurones réside dans le caractère non explicite des résultats. La
convergence n’est pas forcément vers une solution optimale globale. Le choix peu guidé de la
structure du PMC(nombre de couches, nombre de neurones par couches). Sensible à un très
grand nombre de variables discriminantes.

Input Output
10000000 → 10000000
01000000 → 01000000
00100000 → 00100000
00010000 → 00010000
00001000 → 00001000
00000100 → 00000100
00000010 → 00000010
00000001 → 00000001

Learned hidden layer representation :

ESSAI - 2017/2018 62 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Input Hidden Output

Values
10000000 → .89 .04 .08 → 10000000
01000000 → .01 .11 .88 → 01000000
00100000 → .01 .97 .27 → 00100000
00010000 → .99 .97 .71 → 00010000
00001000 → .03 .05 .02 → 00001000
00000100 → .22 .99 .99 → 00000100
00000010 → .80 .01 .98 → 00000010
00000001 → .60 .94 .01 → 00000001

7.4 Statistique et PMC

ESSAI - 2017/2018 63 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

7.5 Logiciels utilisés

7.5.1 nnet de R pour le perceptron multi-couches
Le logiciel R est un logiciel de statistiques libre. Il est constitué de packages permettant de
réaliser des études statistiques. Le package nnet est utilisé pour le Percepteur Multi Couches.
Nous allons le comparer avec le modèle polynomial d’interpolation. Pour cela, nous considérons
l’exemple introductif donné par Annick Valibouze à l’ISUP lors du cours sur les Réseaux de
Neurones Artificiels sous le Logiciel R 3 .

Nous considérons une fonction sin(x) bruitée et nous comparons avec le modèle polynomial
(interpolation) ainsi que le PMC.

library(nnet)
library(MASS)
# Les donn\’ees
x <-sort(10*runif(50))
y <- sin(x)+0.2*rnorm(x)
data <- [Link](x,y)
# Construction du PMC linout=TRUE pour \^etre entre 0 et 1
nn <- nnet(x,y,size=4,maxit=100,linout=TRUE)
plot(x,y,col="blue",pch=16)
# Le PMC en mode de fonctionnement
lines(x1,predict(nn,[Link](x=x1)),col="red")
# La courbe sin(x)
x1<-seq(0,10,by=0.01)
lines(x1,sin(x1),col="green")
3. Troisième année ISUP, Filière GRIE, Master 2 De Statistique, Major Data Science

ESSAI - 2017/2018 64 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

# Comparaison avec le mod\‘ele polynomial

modelPoly <- lm(y$\tilde$ x + I(x$^$2) + I(x$^$3) + I(x$^$4))
lines(x,predict(modelPoly), col="orange")
# L\’egende
legend(3,1.3,c("donnees simulees avec bruit","fonction sin(x)","PMC",
"modele polynomial"),lty=c(0,1,1,1),pch=c(16,-1,-1-1),
col=c("blue","green","red","orange"))
x #Les 50 valeurs de x
summary(nn) # le réseau de neurones

Si vous disposez d’une licence SAS, alors il y a moyen de manipuler des réeaux de neurones
artificiels via le package Proc Neural.

7.5.2 SageMath, Python

7.6 Conclusion
Un réseau de neurones est avant tout un outil de traitement de l’information. En le formant via
l’apprentissage, il peut être capable d’effectuer des tâches de econnaissance, de classification,
d’approximation, de prévision, etc. Nous utilisons par exemple les réseaux de neurones pour
réduire l’écho dans les communications téléphoniques outre-mer (réseau Adaline). Ils peuvent
aussi servir à éliminer du bruit, dans les prévisions météorologiques, dans la compression de
données, dans la modélisation du marché monétaire, dans le développement de l’Intelligence
Artificielle, etc. Il y a aussi moyen de leur apprendre à reconnaître des faces, le son de la voix
et autres lésions cancéreuses. Voir par exemple le projet NetTalk qui est un réseau de 309
neurones relié en entrée à un scanner et branché en sortie à un haut parleur. NetTalk a appris
à lire un texte à haute voix après avoir appris à reconnaître 50000 mots d’un dictionnaire toute
une nuit. Une fois entraîné, le réseau est capable de lire de nouveaux mots qu’il n’ a jamais vu
auparavent. De nos jours, ce sont de vrais neurones animales qui sont utilisés : en Angleterre ce
sont des neurones de rats, mises en culture et connectés par des électrodes au système moteur
d’une machine, qui apprend à éviter les murs à force de se cogner.

ESSAI - 2017/2018 65 inestej@[Link]

Annexe A

Travaux Pratiques sur SAGEmath

Les codes suivants représentent des exercices corrigés ou des Projets de Fin d’années que nous
avons encadrés en Analyse Numérique. Il s’agit de problèmes liés à l’interpolation polynomiale,
à la recherche de solutions d’équations non linéaires et aux problèmes d’Analyse Numérique
matricielle.

A.1 Interpolation
TP2 : Estimation de l’erreur

Théorème : Soient f une fonction de classe C n+1 sur un intervalle [a, b] et P le polynôme de
Lagrange de f en les points x0 < x1 < · · · < xn ∈ [a, b]. Alors

(x − x0 )(x − x1 ) . . . (x − xn ) (n+1)
f (x) − P (x) = f (ζ) (A.1.1)
(n + 1)!

où a ≤ min(x, x0 ) < ζ < max(x, xn ) ≤ b.

Démonstration : La relation est triviale pour x = xi ∀i = 0..n. Soit x 6= xi , pour tout i = 0..n.
f (x)−P (x)
On pose k(x) = (x−x0 )(x−x 1 )...(x−xn )
et w(t) un polynôme qui s’annule en x0 , x1 , . . ., xn , et x
définit par : w(t) = f (t) − P (t) − (t − x0 )(t − x1 ) . . . (t − xn )k(x). Il s’annule en (n + 2) points.
En appliquant le thórème de Rolle, on montre que wn+1 (t) s’annule en au moins un point ζ :

f (n+1) (ζ)
w(n+1) (ζ) = f (n+1) (ζ) − (n + 1)! k(x) = 0 et k(x) = .
(n + 1)!

Cette quantité ne tend pas nécessairement vers 0 car les dérivées f (n+1) de f peuvent grandir
très vite lorsque n croît. Dans l’exemple suivant, nous présentons deux cas sur la convergence
de l’interpolation de Lagrange.
Exemple : Soit f (x) = sin(x), | f k (x) | ≤ 1 : l’interpolé P converge vers f quelque soit
1
le nombre de points d’interpolation et leur emplacement. Par contre, pour f (x) = 1+14x 2 sur

[−1, 1] et bien que f soit indéfiniment continûment dérivable sur [−1, 1], les grandeurs

max−1≤x≤1 | f k (x) | , k = 1, 2, 3, . . .

66
Cours d’Analyse Numérique Ines Abdeljaoued Tej

explosent très rapidement et ceci ne nous assure plus la convergence de l’interpolation.

Le choix des points xi apporte une amélioration sensible de l’interpolation : Lorsque a ≤ x ≤ b,

on choisit les abscisses d’interpolation
a+b b−a 2i + 1
xi = + cos( π)
2 2 2(n + 1)
pour i = 0..n, on obtient :
(b − a)n+1
| f (x) − P (x) | ≤ 2n+1
maxx∈[a,b] | f (n+1) (x) | .
(n + 1)!2
Les points xi = a+b
2
+ b−a
2
2i+1
cos( 2(n+1) π) pour i = 0..n sont calculés à partir des zéros des poly-
nômes de Tchebychev (Tn (x) = vérifiant Tn+1 (x) = 2xTn (x) − Tn−1 (x)).

Exercice à faire sur SAGEmathmath :

1. Soit f (x) = sin(x) et une subdivision de l’intervalle [−1, 1] en n = 5 sous-intervalles
equidistants. Posons x0 = −1, x1 = −0.5, x2 = 0, x3 = 0.5 et x4 = 1.
Implémenter l’algorithme de Lagrange se basant sur les polynômes d’interpolation de
base. On prendra fi = f (xi ) pour i = 0..n.
Dessiner dans un même graphique la fonction f , le polynôme P ainsi que les points
(xi , fi ). Reprendre les calculs pour n = 13.
1 2i+1
2. Refaire le même travail avec la fonction f (x) = 1+14x 2 , puis avec xi = cos( 2n+2 π) pour

i = 0..n et n = 5 puis n = 13.

Correction sur SAGEmathmath :

1. Le code sous SAGEmathmath est donné ci-dessous :

SAGEmath: from pylab import arange

SAGEmath: xdata=arange(-1.0,1.1,0.5)
SAGEmath: f=sin
SAGEmath: P=lagrange(f, xdata)
SAGEmath: R=plot(f, -1,1, rgbcolor=(0,2,1))
SAGEmath: Q=plot(P, -1,1, rgbcolor=(0,1,2))
SAGEmath: n=len(xdata)-1
SAGEmath: s=Graphics()
SAGEmath: for i in range(n+1):
... xi = xdata[i]
... fi = f(xdata[i])
... s = s + point((xi,fi), rgbcolor = (1,0,0))
SAGEmath: R+Q+s
Les graphiques suivants donnent la fonction sin et son polynoôme de Lagrange sur n = 5 puis
sur n = 13 points :
1
2. Si on calcule le polynôme de Lagrange sur la fonction f = 1+14x 2 , nous obtenons :

ESSAI - 2017/2018 67 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Figure A.1 – Interpolation pour n =

5 Figure A.2 – Interpolation pour n = 13

SAGEmath: from pylab import arange

SAGEmath: xdata=arange(-1.0,1.1,0.5)
SAGEmath: f=1/(1+14*x^2)
SAGEmath: P=lagrange(f, xdata)
SAGEmath: R=plot(f, -1,1, rgbcolor=(0,2,1))
SAGEmath: Q=plot(P, -1,1, rgbcolor=(0,1,2))
SAGEmath: n=len(xdata)-1
SAGEmath: s=Graphics()
SAGEmath: for i in range(n+1):
... xi = xdata[i]
... fi = f(xdata[i])
... s = s + point((xi,fi), rgbcolor = (1,0,0))
...
SAGEmath: R+Q+s

Figure A.3 – Interpolation pour n =

5 Figure A.4 – Interpolation pour n = 13

Enfin, si on considère à la place des abscisses xi équidistants, les zéros de Tchebychev, nous
obtenons la figure 5 pour n = 5 points et ensuite la figure 6 pour n = 13 :
SAGEmath: a=-1; b=1; n=5;
SAGEmath: ydata=[simplify((a+b)/2+(b-a)/2*cos((2*k+1)*3.14/(2*n+2)))
for k in range(n+1)]

ESSAI - 2017/2018 68 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

SAGEmath: print(ydata)
[0.965960168538399, 0.707388269167200, 0.259459981914882, -0.257921542147459,
-0.706261644820005, -0.965546938710468]
SAGEmath: f=1/(1+14*x^2)
SAGEmath: P=lagrange(f, ydata)
SAGEmath: A=plot(f, -1,1, rgbcolor=(0,2,1))
SAGEmath: B=plot(P, -1,1, rgbcolor=(1,0,0))
SAGEmath: A+B

Figure A.5 – Interpolation sur les zé- Figure A.6 – Interpolation sur les zéros de
ros de Tchebychev pour n = 5 Tchebychev pour n = 13

SAGEmath: a=-1; b=1; n=13;

SAGEmath: ydata=[simplify((a+b)/2+(b-a)/2*cos((2*k+1)*3.14/(2*n+2)))
for k in range(n+1)]
SAGEmath: print(ydata)
[0.993718576879459, 0.943939675865892, 0.846875476196380, 0.707388269167200,
0.532465465533183, 0.330869571228815, 0.112699242252689, -0.111116592961247,
-0.329366202474389, -0.531116686408463, -0.706261644820005, -0.846027443250886,
-0.943412715311215, -0.993539086045688]
SAGEmath: f=1/(1+14*x^2)
SAGEmath: P=lagrange(f, ydata)
SAGEmath: A=plot(f, -1,1, rgbcolor=(0,2,1))
SAGEmath: B=plot(P, -1,1, rgbcolor=(1,0,0))
SAGEmath: A+B

A.2 Résolution d’équations

Considérons le polynôme f (x) = x4 + 6x2 − 60x + 36 sur [0, 4]. Il s’agira d’appliquer la méthode
du point fixe pour les fonctions suivantes :
4 2 +36
— g1 (x) = x +6x60

−36
— g2 (x) = x3 +6x−60
1
— g3 (x) = (−6x2 + 60x − 36) 4

ESSAI - 2017/2018 69 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

1. Ecrire une procédure qui prend en entrée g, x0 et N et qui calcule le Nème terme de la
liste xn+1 = g(xn ).
2. Vérifier ensuite la convergence de la méthode du point fixe pour chacun des gi sur [0, 4].
Pour cela, il faut chercher à satisfaire toutes les conditions du théorème de convergence
de la méthode du point fixe :
3. Ecrire un script qui permet d’accélérer la convergence de g3 .

A.3 Méthode d’Euler

Partant de x0 = a, on connaît y0 = y(x0 ) = τ . L’idée de la méthode d’Euler est de considérer que
la courbe y sur [x0 , x1 ] n’est pas très éloignée de sa tangente en x0 : une estimation rationnelle de
y(x1 ) est y1 = y0 + hf (x0 , y0 ). En effet, pour x1 , l’équation de la tangente est y 0 (x0 )(x − x0 ) + y0 .
On estime alors que y 0 (x1 ) n’est pas très différent de f (x1 , y1 ). Sur [x1 , x2 ], on remplace la courbe
par sa tangente approchée en x1 ce qui donne l’estimation y2 de y(x2 ) : y2 = y1 + hf (x1 , y1 )
et on poursuit de la même manière. La méthode d’Euler est d’ordre 1. Une telle méthode ne
converge pas assez vite pour donner des résultats pratiques intéressants. On va étudier dans les
deux sections suivantes deux méthodes d’ordre plus élevé. La méthode d’Euler est implémentée
dans SAGEmath comme suit :

def euler(f,n,tau,a,b):
h = (b-a)/n
x = [x0]
y = tau
for i in range(n):
x += [x[i]+h]
y += [y[i]+h*f(x[i],y[i])]
return x,y

A.4 Analyse Numérique Matricielle : SVD

«L’imagerie numérique a longtemps été un exemple d’application de la SVD. Maintenant, il
existe des moyens plus élaborés et plus efficaces. Reste que la SVD est un outils intéressant
pour le transfert d’image.
Une photo en noir et blanc peut être assimilée à une matrice A dont les éléments correspondent
à des niveaux de gris pour les différents pixels qui constituent l’image. Chaque ai,j donne le
niveau de gris du pixel (i, j) de l’image : ai,j ∈ [0, 1] avec par exemple ai,j = 0 pour un pixel
(i, j) blanc et ai,j = 1 pour un pixel (i, j) noir».

Ci-dessous une image en noir et blanc de taille 2500 X 4000 pixels dont la matrice A représen-
tative est générée par un logiciel de lecture de l’image (via un scanner ou un appareil photo
numérique). le résultat renvoyé correspond à une matrice A de rang r. la matrice admet donc
r valeurs singulières non nulles. Soit σ1 la valeur maximale de ces valeurs singulières et soient
σ1 , σ2 , ...σr ces valeurs singulières triées par ordre décroissant. On sait alors qu’il existe une
matrice U de taille 2500, une matrice V de taille 4000 telsque A = U diag(σ1 , σ2 , ...σr ) V t . En
désignant par ui les vecteurs colonnes de U , et par vi les vecteurs colonnes de V , on obtient

ESSAI - 2017/2018 70 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Pr
que A = i=1 σi ui vit .

La connaissance de A elle-même correspond à la connaissance des 2500 X 4000 coefficients de

cette matrice. Si on regarde l’expression de droite, la connaissance de A est donnée par les r
valeurs singulières de A, plus r vecteurs de taille 2500 et r vecteurs de taille 4000.

Mais si on regarde de plus près les valeurs singulières non nulles σ1 , σ2 , ...σr , on constate que
nombreuses d’entre-elles sont très petite (voir négligeable) devant les autres. On peut alors
décider que σσ1i > , où  est un filtre choisi par l’utilisateur, par exemple  = 10−3 . Dans ce
cas, σi est considérée comme nulle : on réduit ainsi le nombre de valeurs singulières à conserver
dans l’expression de A selon la décomposition SVD.

Notons σ1 , σ2 , ...σk les valeurs singulières conservées. Dans l’exemple proposé ci-dessous, avec
k = 3, l’image est floue mais on constate facilement que c’est un paysage. Pour k = 13 l’image
est déjà ressemblante à l’originale. pour k = 80 on peut dire que l’image obtenue est identique à
l’originale. La connaissance de A est alors donnée par la connaissance de 80 valeurs singulières,
plus 80 vecteurs de taille 2500, plus de 80 valeurs de taille 4000. Soit 520 080 au lieu de 10 000
000.

A.4.1 A faire sur Python :

1. Donnez la décomposition en valeurs sigulières de la matrice A = 1, 1, 1, −1
2. Trouvez la commande qui permet d’effectuer la décomposition en valeurs singulières de
A sur SAGEmath
3. Téléchargez une image sous format .png

ESSAI - 2017/2018 71 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

4. Calculez la SVD de cette image en utilisant les logiciels numpy et pylab.

A.4.2 Correction
%matplotlib inline
import [Link] as plt
import numpy as np
import time

from PIL import Image

img = [Link](’input/img_6847.jpg’) #ici l’image chargée est input/img_6847.jpg

imggray = [Link](’LA’)
[Link](figsize=(9, 6))
[Link](imggray);

#On la transforme en une matrice numpy

imgmat = [Link](list([Link](band=0)), float)
[Link] = ([Link][1], [Link][0])
imgmat = [Link](imgmat)
[Link](figsize=(9,6))
[Link](imgmat, cmap=’gray’);

#Le calcul de la SVD se fait comme suit :

U, sigma, V = [Link](imgmat)

#Ici la reconstruction de l’image

reconstimg = [Link](U[:, :1]) * [Link](sigma[:1]) * [Link](V[:1, :])
[Link](reconstimg, cmap=’gray’);

for i in xrange(2, 4):

reconstimg = [Link](U[:, :i]) * [Link](sigma[:i]) * [Link](V[:i, :])
[Link](reconstimg, cmap=’gray’)
title = "n = %s" % i
[Link](title)
[Link]()

for i in xrange(5, 51, 5):

reconstimg = [Link](U[:, :i]) * [Link](sigma[:i]) * [Link](V[:i, :])
[Link](reconstimg, cmap=’gray’)
title = "n = %s" % i
[Link](title)
[Link]()

ESSAI - 2017/2018 72 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

A.4.3 Un code plus intéracif

import numpy
import pylab
U, s, V = [Link](matrix(2,2,[1,1,1,-1]))
U, s, V

import pylab
import numpy

# Read in a png image as a *numpy array* : [Link]

# Convert to greyscale : [Link](...,2)
#A_image = [Link]([Link](DATA + ’[Link]’), 2)
#A_image = [Link]([Link](DATA + ’CapturFiles_11.png’), 2)
A_image = [Link]([Link](DATA + ’[Link]’), 2)

# Use numpy to compute the singular value decomposition of the image

u,s,v = [Link](A_image)
S = [Link](A_image.shape)
S[:len(s),:len(s)] = [Link](s)
n = A_image.shape[0]

# The SVD is a factorization of A_image as a product: U * S * V

# We verify that this is indeed the case.
# WARNING u * S is *NOT* matrix multiplication in numpy,
# it’s componentwise multiplication. Instead use "[Link]" to multiply.

@interact
def svd_image(i = ( "Eigenvalues(quality)",(20,(1..A_image.shape[0])))):

# Thus the following expression involves only the first i eigenvalues,

# and the first i rows of v and the first i columns of u.
# Thus we use i*m + i*n + i = i*(m+n+1) storage, compared to the
# original image that uses m*n storage.
A_approx = [Link]([Link](u[:,:i], S[:i,:i]), v[:i,:])
p = show(graphics_array([matrix_plot(A_image), matrix_plot(A_approx)]), axes=False

# Of course, nothing is ever exactly right with floating point matrices...

# so we check for closeness.
print [Link](A_approx, A_image, 1e-1, 1e-1)

Il s’agit d’abord de tester la notion de problème bien ou mal posé par la voie des matrices bien
ou mal conditionnées. Ensuite, il faudra implémenter la factorisation LU et la factorisation de
Cholesky et de comparer ses deux factorisations.

ESSAI - 2017/2018 73 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

A.5 Matrice bien ou mal conditionnée

Considérons les deux systèmes suivants : Ax = b avec respectivement
   
4.218613 6.327917 10.546530
A= ,b =
3.141592 4.712390 7.853982

et    
4.218611 6.32791 10.546530
A= ,b =
3.141594 4.712390 7.853980
1. Dessinez les deux droites associées au premier système. Vous remarquerez que le premier
système est formé de deux droites presque parallèles. Ainsi, résoudre le système revient
à chercher l’intersection de ces deux droites presque parallèles.
2. Résoudre ces deux systèmes et comparer les solutions. Que remarquez-vous ? Il est clair
que si on perturbe le deuxième système, l’intersection de ces deux nouvelles droites
presque parallèles sera fortement modifié.
3. Calculer le nombre de conditionnement ||A||.||A−1 || de chacun de ces système. Le condi-
tionnement étant loin de 1, ceci explique les différences des solutions. On parle de pro-
blème mal conditionné.

A.6 Factorisation LU
On considère la matrice  
4 8 12
A=  3 8 13 
2 9 18
1. Ecrire l’algorithme de factorisation LU.
2. Comparer le résultat avec celui de la fonction python qui permet une factorisation LU.

A.7 Factorisation de Cholesky

On considère une matrice A tridiagonale, symétrique et définie-positive d’ordre N . On cherche
à compléter un algorithme permettant de résoudre un système linéaire Ax = b par la méthode
de Cholesky. Notons ai les éléments de la diagonale de A et cj les éléments de la sur-diagonale
et sous-diagonale de A pour i = 1, . . . , N et j = 1 . . . , N − 1 :
 
a1 c 1
 c1 a2 c2 
 
 c2 . . . ... 
A=  
 · · · c N −2


 cN −2 aN −1 cN −1 
cN −1 aN

ESSAI - 2017/2018 74 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Soit R une matrice triangulaire inférieure telle que A = RRt . Alors R aura la forme suivante :
 
l1
 m1 l2 
 
 ... ... 
R=  
 · · · c N −2


 · · · lN −1 
mN −1 lN

1. Trouver le l et m en fonction de a et de c.
2. Compléter l’algorithme de décomposition de Cholesky suivant (l est mémorisé dans a et
m
 est mémorisé pdans c) :

 a(1) = a(1)
Pour i allant de 1 à N-1 Faire


 c(i) = ?/ ? √
a(i + 1) = ?−?∗?


3. La résolution de Ax = b nécessite la résolution du système linéaire triangulaire inférieur

Ry = b, puis du système triangulaire supérieur Rt x = y. Compléter l’algorithme de sorte
à obtenir la solution x du système linéaire Ax = b :

par une méthode de descente de Ry = b :

(a) 
 b(1) = b(1)/a(1)
Pour i allant de 1 à N-1 Faire
b(i + 1) = (b(i + 1)−?∗?)/?


par une méthode de remontée de Rt x = y :

(b) 
 b(N ) = b(N )/a(N )
Pour i allant de N-1 à 1 (avec un pas de -1) Faire
b(i) = (b(i)−?∗?)/?


Les éléments non nuls de la matrice R correspondent-ils aux éléments non nuls de la
matrice A ?

Application en Probabilité - Simulation de la loi N (m, Γ)). Ici m ∈ Rd et Γ est

une matrice de Md (R), symétrique et définie-positive. Soient C une matrice telle que
CC t = Γ et Z un vecteur aléatoire gaussien de loi N (0, Id ). Alors X = CZ + m a pour
loi N (m, Γ)).
Correction
SAGEmath: N =3
SAGEmath: a = [3.0 for i in range(N)]
SAGEmath: c = [-1.0 for i in range(N-1)]
SAGEmath: A = matrix(RR,N,N)
SAGEmath: for i in range(N):
... A[i,i]=a[i]
...
SAGEmath: for i in range(N-1):
... A[i,i+1]=c[i]

ESSAI - 2017/2018 75 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

...
SAGEmath: for i in range(1,N):
... A[i,i-1]=c[i-1]
...
SAGEmath: show(A)
 
3.00000000000000 −1.00000000000000 0.000000000000000
 −1.00000000000000 3.00000000000000 −1.00000000000000 
0.000000000000000 −1.00000000000000 3.00000000000000
show(A.cholesky_decomposition())
 
1.73205080756888 0.000000000000000 0.000000000000000
 −0.577350269189626 1.63299316185545 0.000000000000000 
0.000000000000000 −0.612372435695794 1.62018517460197
SAGEmath: def cholesky(a,c):
... a[0]=sqrt(a[0]).n()
... for i in range(N-1):
... c[i]=c[i]/a[i].n()
... a[i+1]=(sqrt(a[i+1]-c[i]^2)).n()
... A = matrix(CC,N,N)
... for i in range(N):
... A[i,i]=a[i]
... for i in range(1,N):
... A[i,i-1]=c[i-1]
... return(A)
SAGEmath: a = [3.0 for i in range(N)]
SAGEmath: c = [-1.0 for i in range(N-1)]
SAGEmath: show(cholesky(a,c))

 
1.73205080756888 0.000000000000000 0.000000000000000
 −0.577350269189626 1.63299316185545 0.000000000000000 
0.000000000000000 −0.612372435695794 1.62018517460197

A.8 Méthodes itératives

Il s’agit de résoudre un système d’équation linéaires par les méthodes développées précédemment
(LU, Cholesky, Jacobi, Gradient conjugué).
1. Ecrire les algorithmes de descentes puis de remontée appliqués à des matrices triangu-
laires (respectivement inférieures et supérieures)
2. Résoudre le système d’équations linéaires Ax = b avec
   
2 -1 1
A= ,b =
-1 2 1

3. Vérifiez que A est symétrique définie-positive. Calculez la factorisation de Cholesky de

A puis résoudre le système Ax = b.

ESSAI - 2017/2018 76 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

4. La méthode de Jacobi consiste à calculer la suite xk+1 = Jxk +D−1 b avec J = In +D−1 A
et D la matrice diagonale contenant les éléments  de la diagonale de A. Appliquez la
3

méthode de Jacobi, en partant de x0 = 2 . Vous remarquerez que x2 = x∗ la
2
solution du système Ax = b.
5. Résoudre le système Ax
 =b grâce à la factorisaton de Cholesky puis grâce à la méthode
0
de Jacobi, avec x0 =  0  :
0
   
4 2 1 4
A=  -1 2 0  ,b =  2 
2 1 4 9

ESSAI - 2017/2018 77 inestej@[Link]

Annexe B

Exerices de Révision et Applications

Exercice B.0.1. Citer deux méthodes d’interpolation polynomiale. Quelle est la différence entre
ces méthodes et la méthode des moindres carrés discrets ?

Exercice B.0.2. Soit f une fonction de R dans R telle que

f (x) = x3 + 3x2 − 1 pour x ∈ R.

1. Montrer que l’équation f (x) = 0 admet une unique solution sur [0, 1].
2. Nous désirons déterminer la solution de l’équation f (x) = 0 avec la méthode de Newton
sur [0, 1].
(a) Est-ce que la fonction f vérifie les conditions du théorème de convergence globale de
la méthode de Newton pour x0 = 1. Donner le théorème.
(b) Calculer la racine de f sur l’intervalle [0, 1] avec une précision de calculs égale à 10−6 .
3. Ecrire l’algorithme général de la méthode de Newton qui prend en entrée les points x0 ,
a < b et une fonction f et rend la racine de f (x) = 0 sur [a, b] ou bien un mesSAGEmath
d’erreur.

Exercice B.0.3 (10pt). Soit f une fonction de R dans R telle que

f (x) = x3 + 3x2 − 1 pour x ∈ R.

1. Montrer que l’équation f (x) = 0 admet une unique solution sur [0, 1].
2. Nous désirons déterminer la solution de l’équation f (x) = 0 avec la méthode de Newton
sur [0, 1]. En quoi consiste cette méthode ?
3. Ecrire l’algorithme général de la méthode de Newton qui prend en entrée les points x0 ,
a < b et une fonction f et rend la racine de f (x) = 0 sur [a, b] ou bien un mesSAGEmath
d’erreur.

78
Cours d’Analyse Numérique Ines Abdeljaoued Tej

4. Est-ce que la fonction f vérifie les conditions du théorème de convergence globale de la

méthode de Newton pour x0 = 1. Donner le théorème.
5. Calculer la racine de f sur l’intervalle [0, 1] avec une erreur de 10−6 .

Exercice B.0.4 (10pt). Soit A une matrice définie par :

 
1 2 -1 2
 2 3 -2 -1 
A=  2 2

1 3 
-1 -2 0 1

1. Donner la définition de la décomposition LU d’une matrice carré.

2. Calculer la factorisation LU de A.
3. En déduire le déterminant de A.
4. Le système linéaire  
1
 0 
Ax =  
 -2  où x ∈ R4
2
admet-il une solution ? est-elle unique ?
5. Résoudre le système en utilisant la factorisation LU de la matrice A.
Exercice B.0.5. 1. Soit f une fonction continue sur un intervalle [a, b] ⊂ R avec f (a).f (b) <
0. Expliquer géométriquement la méthode de Newton dans la recherche d’une solution
de l’équation f (x) = 0 sur [a, b].
2. Donner l’algorithme de la méthode de Newton pour la recherche d’une racine simple
d’une fonction de classe C 2 sur [a, b].

Exercice B.0.6. Nous considérons l’ensemble des points

x0 = −1, x1 = 0, x2 = 1, x3 = 2,
f0 = 1, f1 = −1, f2 = 3, f3 = 19.

1. Afin de déterminer une approximation du nuage de points (xi , fi ) pour i de 0 à 3, calculer

la droite de régression linéaire P0 par la méthode des moindres carrés discrets.

2. Citer deux méthodes permettant de calculer le polynôme de Lagrange.

3. Montrer que les polynômes d’interpolation de Lagrange vérifient Li (xi ) = 1 et Li (xj ) = 0
pour tout i et j 6= i de 0 à 3.
4. Calculer le polynôme de Lagrange P sur les points (x0 , f0 ), (x1 , f1 ), (x2 , f2 ) et (x3 , f3 ).
5. Quelle est la différence entre P et P0 .

6. Donner le théorème de convergence globale de la méthode de Newton pour une fonction

de classe C 2 .

ESSAI - 2017/2018 79 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

7. Calculer la racine de P (x) = 0 sur l’intervalle [0, 3] par la méthode de Newton avec
x0 = 1. La précision des calculs est à 10−6 près.
Exercice B.0.7. 1. Soit f une fonction continue sur un intervalle [a, b] ⊂ R avec f (a).f (b) <
0. Expliquer géométriquement la méthode de la sécante dans la recherche d’une solution
de l’équation f (x) = 0 sur [a, b].
2. Donner l’ordre de la méthode de Newton pour la recherche d’une racine simple d’une
fonction de classe C 2 .

Exercice B.0.8. Nous considérons l’ensemble des points

x0 = −1, x1 = 0, x2 = 1, x3 = 2,
f0 = −4, f1 = −1, f2 = −2, f3 = 5.

1. Afin de déterminer une approximation du nuage de points (xi , fi ) pour i de 0 à 3, calculer

la droite de régression linéaire P0 par la méthode des moindres carrés discrets.

2. Citer deux méthodes permettant de calculer le polynôme de Lagrange.

3. Montrer que les polynômes d’interpolation de Lagrange vérifient Li (xi ) = 1 et Li (xj ) = 0
pour tout i et j 6= i de 0 à 3.
4. Calculer le polynôme de Lagrange P sur les points (x0 , f0 ), (x1 , f1 ), (x2 , f2 ) et (x3 , f3 ).
5. Quelle est la différence entre P et P0 .

6. Vérifier que P (xi ) = fi pour i de 0 à 3 et montrer qu’il existe une unique racine de P
sur l’intervalle [0, 3].
7. Donner le théorème de convergence globale de la méthode de Newton pour une fonction
de classe C 2 .
8. Calculer la racine de P (x) = 0 sur l’intervalle [0, 3] par la méthode de Newton avec
x0 = 3. La précision des calculs est à 10−6 près.
9. Ecrire l’algorithme de Newton qui prend en entrée les points x0 , a < b et une fonction f
et rend la racine de f (x) = 0 sur [a, b] ou bien un mesSAGEmath d’erreur.
Exercice B.0.9. 1. Calculer par la méthode de la puissance, la valeur propre dominante à
−2
10 près de    
4 -1 -1 1
A=  -1 2 -1  avec x0 =  1 .
-1 -1 2 1
2. Ecrire l’algorithme de la méthode de la puissance inversée.

Exercice B.0.10. Etant donnés une matrice A triangulaire inférieure de taille (n, n) avec des
1 sur la diagonale et un vecteur b de taille n, calculer la complexité en nombre d’opérations
élémentaires de la méthode de descente appliquée au système Ax = b.

Exercice B.0.11. Soit n = 3. Nous considérons l’ensemble de points

ESSAI - 2017/2018 80 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

x0 = −1, x1 = 0, x2 = 1, x3 = 2,
f0 = 1, f1 = 0, f2 = 1, f3 = 4.

1. Afin de déterminer une approximation du nuage de points (xi , fi ) par la droite de ré-
gression linéaire, nous appliquerons la méthode des moindres carrés discrets :
(a) Soit A la matrice associée au système des moindres carrés :
 Pi=n 2 Pi=n 
i=0 xi i=0 xi
A= 
Pi=n
i=0 xi n+1

(b) Donner la décomposition LU de A

   Pi=n 
a0 i=0 xi fi
(c) Résoudre le système A  =  par une descente puis une remon-
Pi=n
a1 i=0 fi
tée.
2. Calculer le polynôme de Lagrange P sur les points (xi , fi ) pour i entre 0 et 3.
3. Quelle est la différence entre P et la droite de régression linéaire ?

4. Donner le théorème de convergence globale de la méthode de Newton pour une fonction

de classe C 2 . Peut-on appliquer ce théorème au calcul de la racine de P (x) = 0 sur
l’intervalle [0, 3] avec x0 = 1 ?

Exercice B.0.12. Etant donnés une matrice A triangulaire supérieure de taille (n, n) inversible
et un vecteur b de taille n, calculer la complexité en nombre d’opérations élémentaires de la
méthode de remontée appliquée au système Ax = b.

ESSAI - 2017/2018 81 inestej@[Link]

Annexe C

Tests, DS et Examens corrigés

C.1 2008-2009
Nom et Prénom : Groupe :
Ecole ou Institut de l’année dernière :

Test 1 : Interpolation et Approximation

Exercice : Soient (xi , fi ) pour i = 0..n des points réels ou complexes. Ecrire un algorithme
permettant de calculer le polynôme de Lagrange se basant sur les polynômes d’interpolation
Li (x) pour i = 0..n.
Déterminer la complexité de cet algorithme en fonction de n.

Correction :
Entrée : (xi , fi ) pour i = 0..n
Sortie : Le polynôme de Lagrange P(x).

P := 1 ;
Pour i allant de 0 à n Faire
L :=1 ;
Pour j allant de 0 à n Faire
Si i <> j alors
L := L + (x − xj )/(x − xi ) ;
Fin Si ;
Fin Pour ;
P :=P + fi + L ;
Fin Pour ;
Afficher(P ) ;
Fin.

Il y a o(4n2 ) opérations élémentaires. Il s’agit de o(2n2 ) additions/soustractions, o(n2 ) multi-

82
Cours d’Analyse Numérique Ines Abdeljaoued Tej

plications et o(n2 ) divisions.

Examen du module Analyse Numérique

Exercice 1 (6pt).
1. Donner le théorème de Gershgorin appliqué à une matrice A de taille n × n quelconque.
2. Expliquer la méthode de la puissance appliqée à une matrice A en donnant l’algorithme.
3. Calculer les deux premières itérations de la méthode de la puissance apliquée à :
!
 
√1
1 0 2
A= avec y0 = √1
.
0 2 2

Exercice 2 (4pt). Exercice : Soient (xi , fi ) pour i = 0..n des points réels ou complexes.
Ecrire un algorithme permettant de calculer le polynôme de Lagrange se basant sur les poly-
nômes d’interpolation Li (x) pour i = 0..n.
Déterminer la complexité de cet algorithme en fonction de n.

Problème (10pt). Soit n = 3. Nous considérons l’ensemble de points

x0 = −1, x1 = 0, x2 = 1, x3 = 2,
f0 = 1, f1 = 0, f2 = 1, f3 = 4.
1. Donner le théorème de convergence globale de la méthode de Newton pour une fonction
de classe C 2 . Peut-on appliquer ce théorème au calcul de la racine de P (x) = 0 sur
l’intervalle [0, 3] avec x0 = 1 ?
 Pi=n Pi=n 
i=0 x2i i=0 xi
2. Soit A =   . Donner la décomposition LU de A.
Pi=n
i=0 xi
n+1
   Pi=n 
a0 i=0 x i f i
3. Résoudre le système A  =  par une descente puis une remontée.
Pi=n
a1 i=0 fi
4. En déduire une approximation du nuage de points (xi , fi ) par la droite de régression
linéaire.
5. Quelle est la différence entre le polynôme de Lagrange et la droite de régression linéaire ?

Correction de l’Examen d’Analyse Numérique

Exercice 1 (6pt).

ESSAI - 2017/2018 83 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

1. Le théorème de Gershgorin appliqué à une matrice A de taille n × n quelconque donne

une borne sur le spectre de A : toute valeur propre appartient à l’un au moins des disques
de Gerschgorin : X
| aii − λ |≤ aij
i6=j

2. La méthode de la puissance appliqée à une matrice A consiste à choisir au hasard un

Ayn
vecteur y0 et à caluler la suite yn+1 = kAynk
.
     
0.71 0.45 0.24
3. Pour y0 = , nous avons y1 = et y2 = .
0.71 0.89 0.97

Exercice 2 (4pt). Entrée : (xi , fi ) pour i = 0..n

Sortie : Le polynôme de Lagrange P(x).

P := 1 ;
Pour i allant de 0 à n Faire
L :=1 ;
Pour j allant de 0 à n Faire
Si i <> j alors
L := L + (x − xj )/(x − xi ) ;
Fin Si ;
Fin Pour ;
P :=P + fi + L ;
Fin Pour ;
Afficher(P ) ;
Fin.
Il y a o(4n2 ) opérations élémentaires. Il s’agit de o(2n2 ) additions/soustractions, o(n2 ) multi-
plications et o(n2 ) divisions.

Problème (10pt). Soit n = 3. Nous considérons l’ensemble de points

x0 = −1, x1 = 0, x2 = 1, x3 = 2,
f0 = 1, f1 = 0, f2 = 1, f3 = 4.
1. Voir le cours pour le théorème de convergence globale de la méthode de Newton pour
une fonction de classe C 2 .
On ne peux appliquer ce thèorème au calcul de la racine de P (x) = x2 sur l’intervalle
00
[0, 3] : il faudra restreindre l’intervalle car P 0 (0) = 0 et P (0) = 0).
 
6 2
2. Soit A =  . La décomposition LU de A est triviale :
2 4
   
1 0 6 2
L= , U = .
1 10
3
1 0 3

ESSAI - 2017/2018 84 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

   
a0 8  
8
3. La résolution du système A  =  par une descente donne 10 puis la
a1 6 3
 
1
remontée donne .
1
4. Une approximation du nuage de points (xi , fi ) est donnée par la droite de régression
linéaire d’équation P0 (x) = x + 1.
5. Conclure sur la variation en terme de degré, de complexité des calculs et en terme de
différence entre interpolation et approximation.

Examen de contrôle du module Analyse Numérique

Problème 1 (10pt). Notons A une matrice symétique définie positive d’ordre n.

√
r1,1 = a1,1
pour i de 2 à n faire :
Pk=j−1
ai,j − rik rjk
ri,j = k=1
rj,j
pour j = 1..i − 1
et q
ai,i − k=i−1 2
P
ri,i = k=1 ri,k .
1. Ecrire l’algorithme de la méthode de Cholesky
2. Supposons que la fonction racine suppose au plus 9 opérations élémentaires. Déterminer
la complexité de la méthode de factorisation de Cholesky de A en fonction de n.
3. Donner les 2 premières itérations de cette méthode appliquée à la matrice
 
1 1
A=
1 2

Problème 2 (10pt). Soit n = 3. Nous considérons l’ensemble de points

x0 = −1, x1 = 0, x2 = 1, x3 = 2,
f0 = 1, f1 = 0, f2 = 1, f3 = 4.
1. Donner le théorème de convergence globale de la méthode de Newton pour une fonction
de classe C 2 . Peut-on appliquer ce théorème au calcul de la racine de P (x) = 0 sur
l’intervalle [0, 3] avec x0 = 1 ?
 Pi=n 2
Pi=n 
i=0 xi i=0 xi
2. Soit A =   . Donner la décomposition LU de A.
Pi=n
i=0 xi n+1

ESSAI - 2017/2018 85 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

   Pi=n 
a0 i=0 xi fi
3. Résoudre le système A  =  par une descente puis une remontée.
Pi=n
a1 i=0 fi
4. En déduire une approximation du nuage de points (xi , fi ) par la droite de régression
linéaire.
5. Quelle est la différence entre le polynôme de Lagrange et la droite de régression linéaire ?

Correction de l’Examen de contrôle du module Analyse Numérique

Problème 1 (10pt). Notons A une matrice symétique définie positive d’ordre n.

1. L’algorithme de la méthode de Cholesky est le suivant :
Entrée : ai,j
Sortie : ri,j
√
r1,1 = a1,1
Pour i allant de 2 à n Faire
ri,i = ai,i
Pour j allant de 1 à i − 1 Faire
ri,j = ai,j
Pour k allant de 1 à j − 1 Faire
ri,j = ri,j − ri,k rj,k
Fin Pour ;
ri,j
ri,j = rj,j
Fin Pour
Pour k allant de 1 à i − 1 Faire
ri,i = ri,i − ri,k ri,k
Fin Pour ;
√
ri,i = ri,i
Fin Pour ;
Fin.
2. Sachant que la fonction racine suppose au plus 9 opérations élémentaires, la complexité
de la méthode de factorisation de Cholesky de A en fonction de n est calculée comme
suit : P
Il y a ni=2 ( i−1
Pn n(n−1)(n−2)
+ n(n−1)
P
j=1 (j − 1)) + i=2 (i − 1) = 6 2
additions.
Pn Pi−1 Pn n(n−1)(n−2) n(n−1)
Il y a i=2P(n j=1 j) + i=2 (i − 1) = 6
+ 2 multiplications.
Et il y a i=2 9 = 9(n − 1). Donc 9(n − 1) opérations élémentaires permettant le calcul
des racines carrés.

Nous obtenons une complexité de l’ordre de o( 31 n3 + 9n).

3. La matrice R étant triangulaireinférieure,
 elle est calculée par la méthode de Cholesky
1 0
appliquée à A. Nous obtenons .
1 1

ESSAI - 2017/2018 86 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Problème 2 (10pt). Soit n = 3. Nous considérons l’ensemble de points

x0 = −1, x1 = 0, x2 = 1, x3 = 2,
f0 = 1, f1 = 0, f2 = 1, f3 = 4.
1. Voir le cours pour le théorème de convergence globale de la méthode de Newton pour
une fonction de classe C 2 .
On ne peux appliquer ce thèorème au calcul de la racine de P (x) = x2 sur l’intervalle
00
[0, 3] : il faudra restreindre l’intervalle car P 0 (0) = 0 et P (0) = 0).
 
6 2
2. Soit A =  . La décomposition P A = LU est triviale :
2 4
   
1 0 6 2
L= , U = .
1 10
3
1 0 3

Ici, la matrice de permutation P est égale à la matrice identité d’ordre 2.

   
a0 8  
8
3. La résolution du système A   =   par une descente donne 10 puis la
a 6 3
1
 
1
remontée donne .
1
4. Une approximation du nuage de points (xi , fi ) est donnée par la droite de régression
linéaire d’équation P0 (x) = x + 1.
5. La variation entre ces deux approches se notent en terme de degré : le polynôme d’inter-
polation ayant un degré égal au nombre de points-1 alors que dans l’approximation, le
degré en est indépendant. Dans le cas de l’approximation, le degré est largement plus pe-
tit que dans l’interpolation. La différence entre les deux approches est constatée dans la
complexité des calculs : le cas étudié en cours, donne une approximation plus rapidement
que l’interpolation de Lagrange par exemple.

C.2 2009-2010

Examen du module Analyse Numérique

Exercice 1 (6pt). Notons A = (aij ) une matrice symétrique définie positive d’ordre n. La
factorisation de la matrice A par la méthode de Cholesky est détaillée comme suit :
√
r1,1 = a1,1
pour i de 2 à n faire :
Pk=j−1
ai,j − rik rjk
ri,j = k=1
rj,j
pour j = 1..i − 1
et q
ai,i − k=i−1 2
P
ri,i = k=1 ri,k .

ESSAI - 2017/2018 87 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

1. Ecrire l’algorithme de la méthode de Cholesky qui prend en entrée A et qui calcule une
matrice R triangulaire vérifiant A = Rt R.
2. Supposons que la fonction racine nécessite au plus 9 opérations élémentaires. Déterminer
la complexité de la méthode de factorisation de Cholesky de A en fonction de n.
3. Donner la factorisation de Cholesky de la matrice
 
1 1
A=
1 2
Exercice 2 (4pt). Calculer par la méthode de la puissance la valeur propre et le vecteur propre
!
√1
 
1 0 2
dominants de la matrice A = avec X0 = √1
.
0 2 2

Problème (10pt). Soit n = 3. Nous considérons l’ensemble de points

x0 = −1, x1 = 0, x2 = 1, x3 = 2,
f0 = 0, f1 = −1, f2 = 0, f3 = 3.
1. Calculer le poynôme de Lagrange P sur les points d’interpolation xi , i = 0..3.
2. Donner le théorème de convergence de la méthode du point fixe pour une fonction de
classe C 1 sur un intervalle [a, b]. Peut-on appliquer ce théorème au calcul de la racine
de P (x) = 0 sur l’intervalle [−2, 0] avec g(x) = 12 P (x) + x et X0 = 0 ? Si oui, calculer la
racine de P sur [−2, 0].
 Pi=n Pi=n     Pi=n 
i=0 x2i i=0 xi a0 i=0 xi f i
3. Soit A =  . Résoudre le système A  = 
Pi=n Pi=n
i=0 xi n + 1 a1 i=0 fi
par une descente puis une remontée (suite à une factorisation LU).
4. En déduire une approximation du nuage de points (xi , fi ) par la droite de régression
linéaire p1 (x) = a0 x + a1 .

Examen de contrôle du module Analyse Numérique

Exercice (8pt).
1. Etant donné une matrice A = (aij ), λ une de ses valeurs propres associé au vecteur
propre v et i l’indice de la coordonné vi de v vérifiant | vi |= max1≤j≤n | vj |. Montrer
que
j=n
X
| λ − aii |≤ | aij | .
j=1,j6=i

2. Calculer la valeur propre et le vecteur propre dominants de A par la méthode de la

puissance :    
1 0 1
A= avec X0 = .
0 2 1

ESSAI - 2017/2018 88 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

 
0
3. Que se passe-t-il lorsqu’on considère X0 = ?
2

Problème (12pt). Soit n = 3. Nous considérons l’ensemble de points

x0 = −1, x1 = 0, x2 = 1, x3 = 2,
f0 = 0, f1 = −1, f2 = 0, f3 = 3.
1. Calculer le poynôme de Lagrange P sur les points d’interpolation xi , i = 0..3.
2. Donner le théorème de convergence de la méthode du point fixe pour une fonction de
classe C 1 sur un intervalle [a, b]. Peut-on appliquer ce théorème au calcul de la racine de
P (x) = 0 sur l’intervalle [0, 2] avec g(x) = 21 P (x) + x et X0 = 0 ?
 Pi=n Pi=n 
i=0 x2i i=0 xi
3. Soit A =   . Donner la décomposition LU de A.
Pi=n
i=0 xi
n+1
   Pi=n 
a0 i=0 xi fi
4. Résoudre le système A  =  par une descente puis une remontée.
Pi=n
a1 i=0 fi
5. En déduire une approximation du nuage de points (xi , fi ) par la droite de régression
linéaire p1 (x) = a0 x + a1 .
6. Quelle est la différence entre le polynôme de Lagrange et la droite de régression linéaire ?

C.3 2010-2011

Devoir Surveillé du module Analyse Numérique

Interpolation.
1. Soient f une fonction de classe C 2 sur [0, 1] et xdata = [x0 , x1 , . . . , xn ] une liste de n + 1
points deux à deux distincts de [0, 1]. Ecrire l’algorithme de la méthode des différences
divisées qui prend en entrée f et xdata et qui calcule le polynôme de Lagrange de f aux
points x0 , x1 , . . . , xn .

Indication
 : la formule des différences divisées est la suivante :
ci,i = f (xi ) pour i = 0, .., n
c −ci,j−1
ci,j = i+1,j xj −xi
pour 0 ≤ i < j ≤ n.

2. Répondre par Vrai ou Faux aux affirmations suivantes et justifier votre réponse :
(a) Il existe une infinité de polynômes de degré > n interpolant les points (xi , f (xi ))i=0,.,.n .
π
(b) Le polynôme de Lagrange de sin(x) sur les points 0, 2
et π est égal à P (x) = x3 +3x2 .

Intégration.

ESSAI - 2017/2018 89 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

R2
3. Calculer par la méthode des trapèzes une approximation T (h) de 0 x2 dx. Considérer
pour cela une subdivision de l’intervalle [0, 2] en 2 sous intervalles de même amplitude.
Pn−1 h
4. Soient h = b−a
n
et xi = a + ih. Sachant que T (h) = i=0 2
(f (xi ) + f (xi+1 )), vérifier que
n−1
h T (h) h X h
T( ) = + f (a + kh + ).
2 2 2 k=0 2

Connaissant T (h), donner le nombre d’opérations élémentaires nécessaires au calcul de

T ( h2 ). Nous supposons qu’une évaluation de f nécessite α multiplications et β additions.
R2
5. En déduire 0 x2 dx pour une subdivision de l’intervalle [0, 2] en 4 sous intervalles de
même amplitude.

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

6. Donner le théorème de convergence globale de la méthode de Newton pour une fonction
de classe C 2 sur un intervalle [a, b].
7. Peut-on appliquer ce théorème au calcul de la racine x∗ de f (x) = x3 − 1 sur l’intervalle
[0.5, 1.5] ? Si oui, déterminer x∗ à 10−1 près.

Analyse Numérique - Correction du DS

Interpolation.
1. Soient f une fonction de classe C 2 sur [0, 1] et xdata = [x0 , x1 , . . . , xn ] une liste de
n + 1 points deux à deux distincts de [0, 1]. L’algorithme suivant reprend la méthode des
différences divisées.
Entrée : Une fonction f et xdata == [x0 , x1 , . . . , xn ]
Sortie : Le polynôme de Lagrange de f aux points x0 , x1 , . . . , xn
# Méthode des différences divisées (formule de Newton)
def newton(f,xdata):
n = len(xdata)-1
c = matrix(QQ,n+1,n+1)
for i in range(n+1):
c[i,i] = f(xdata[i])
#print c
#print("-------")
for k in range(1,n+1):
for i in range(n+1-k):
c[i,i+k] = (c[i+1,i+k]-c[i,i+k-1])/(xdata[i+k]-xdata[i])
#print c
#print("-------")
P = c[0,0]
P += sum(c[0,k]*prod((x-xdata[j]) for j in range(k)) for k in range(1,n+1))
return(expand(P))
2. Répondre par Vrai ou Faux aux affirmations suivantes et justifier votre réponse :

ESSAI - 2017/2018 90 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

(a) Vrai. Il existe une infinité de polynômes de degré supérieur strictement à n interpo-
lant les points (xi , f (xi ))i=0,.,.n . En particulier, par deux points passe une infinité de
courbes, mais une seule est une droite.
(b) Faux. Le polynôme de Lagrange de sin(x) sur les points 0, π2 et π n’est pas égal à
P (x) = x3 + 3x2 . Ici, P est un polynôme de degré supérieur au nombre de points +
1 et P (xi ) 6= sin(xi ).

Intégration.
R2
3. Une approximation T (h) de 0 x2 dx par la méthode des trapèzes avec une subdivision
de l’intervalle [0, 2] en 2 sous intervalles de même amplitude est égale à T (1) = 3.
4. Soient h2 = b−a2n
et xi = a + i h2 pour i = 0, . . . , 2n. Il suffit de regrouper les indices i dans
T ( h2 ) en deux groupes : le premier constitué d’indices pairs (donnant T (h) 2
) et le second
h
Pn−1 h
groupe d’indices impairs (donnant 2 k=0 f (a + kh + 2 )). Ainsi, nous obtenons que
n−1
h T (h) h X h
T( ) = + f (a + kh + ). (C.3.1)
2 2 2 k=0 2

Connaissant T (h), Le nombre d’opérations élémentaires nécessaires au calcul de T ( h2 )

est égal à
n(β + 2) + 1 additions 0 soustractions
n(α + 1) multiplications 2 divisions
où α multiplications et β additions sont nécessaires pour une seule évaluation de f .
R2
5. Remplacer a = 0, b = 2, f (x) = x2 , n = 2 et h = 1 dans (C.3.1). On obtient 0 x2 dx
pour une subdivision de l’intervalle [0, 2] en 4 sous intervalles de même amplitude, i.e.
2.75 :

1
Figure C.1 – Intégration pour h = 1 Figure C.2 – Intégration pour h = 2

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

6. Le théorème de convergence globale de la méthode de Newton. Soit f une fonction de
classe C 2 sur un intervalle [a, b] vérifiant f (a).f (b) < 0. Si pour tout x ∈ [a, b], f 0 (x) 6= 0
et f ”(x) 6= 0 alors ∀x0 ∈ [a, b] tel que f (x0 ).f ”(x0 ) > 0, la suite xn+1 = xn − ff0(x n)
(xn )
converge vers l’nique racine de f sur [a, b].

ESSAI - 2017/2018 91 inestej@[Link]

Cours d’Analyse Numérique Ines Abdeljaoued Tej

Figure C.3 – Résolution de x3 − 1

7. La fonction f (x) = x3 − 1 vérifie les conditions du théorème sur l’intervalle [0.5, 1.5].
Elle est strictement monotone et elle change de signe. Elle admet donc une unique racine
obtenue, par la méthode de Newton, de la manière suivante : x∗ à 10−1 près. Partant
de x0 = 1.5, nous obtenons x1 = 1.1481, x2 = 1.0182 et x3 = 1.0003. La courbe rouge
représente f 0 et la courbe bleue représente f ”.

ESSAI - 2017/2018 92 inestej@[Link]