Programmation Analyse Num Erique: Licence 2 - Maths Info

Université de la Polynésie française
Programmation
Analyse numérique
Licence 2 - Maths Info
Benjamin Graille
Université Paris-Saclay
15 avril 2022
Table des matières
1 Introduction et généralités 5
1 Nécessité de l’analyse numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2 Rappels d’analyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1 Régularité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Formules de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.3 Convexité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2 Interpolation polynomiale 11
1 Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2 Exemples avec peu de points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1 Interpolation à un seul point . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Interpolation à deux points : exemple de la droite . . . . . . . . . . . . . . 13
2.3 Interpolation à trois points . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3 Polynôme interpolateur de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.1 Définition et propriétés . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.2 Expression dans la base canonique . . . . . . . . . . . . . . . . . . . . . . 18
3.3 Formule de Lagrange et poids barycentriques . . . . . . . . . . . . . . . . 20
3.4 Comportement asymptotique lorsque n tend vers l’infini . . . . . . . . . . 23
3 Intégration numérique 27
1 Formules de quadratures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
1.1 Définitions et premières propriétés . . . . . . . . . . . . . . . . . . . . . . 28
1.2 Formules de quadratures élémentaires et composées . . . . . . . . . . . . 31
2 Méthodes de quadrature classiques . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.1 Une liste des méthodes classiques . . . . . . . . . . . . . . . . . . . . . . . 33
2.2 Analyse des méthodes des rectangles . . . . . . . . . . . . . . . . . . . . . 35
2.3 Analyse de la méthode des trapèzes . . . . . . . . . . . . . . . . . . . . . 37
2.4 Analyse de la méthode du point milieu . . . . . . . . . . . . . . . . . . . . 38
2.5 Analyse de la méthode de Simpson . . . . . . . . . . . . . . . . . . . . . . 40
2.6 Contrôle de l’erreur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3 Les algorithmes des méthodes de quadrature classiques . . . . . . . . . . . . . . . 43
4 Résolution d’équations ordinaires 45

1 Exemples d’application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
1.1 Schémas numériques pour équations différentielles ordinaires . . . . . . . 45
1.2 Méthode de tir pour les problèmes aux limites du second ordre . . . . . . 47
2 Position correcte du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.1 Existence de solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.2 Notion de conditionnement . . . . . . . . . . . . . . . . . . . . . . . . . . 48
2.3 vitesse de convergence - ordre de convergence . . . . . . . . . . . . . . . . 50
3 Méthodes de type encadrement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.1 méthode de la dichotomie . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.2 Méthode de la sécante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4 Méthodes de type interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
⁄
4.1 méthode de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.2 méthode de la fausse position . . . . . . . . . . . . . . . . . . . . . . . . . 60
Bibliographie 63
1 Introduction et généralités
Dans ce cours, nous décrivons les premiers outils de l’analyse numérique et nous proposons une
implémentation en python des différents modèles numériques rencontrés.
L’analyse numérique est une discipline à l’interface des mathématiques et de l’informatique. Elle
s’intéresse tant aux fondements qu’à la mise en pratique des méthodes permettant de résoudre,
par des calculs purement numériques, des problèmes d’analyse mathématique.
Plus formellement, l’analyse numérique est l’étude des algorithmes permettant de résoudre numé-
riquement par discrétisation les problèmes de mathématiques continues (distinguées des mathé-
matiques discrètes). Cela signifie qu’elle s’occupe principalement de répondre de façon numérique
à des questions à variable réelle ou complexe comme l’algèbre linéaire numérique sur les champs
réels ou complexes, la recherche de solution numérique d’équations différentielles et d’autres
problèmes liés survenant dans les sciences physiques et l’ingénierie. Branche des mathématiques
appliquées, son développement est étroitement lié à celui des outils informatiques.
D’une manière assez caricaturale, nous pouvons définir l’analyse numérique comme l’ensemble
des mathématiques appliquées dont l’objet d’étude est centré sur des algorithmes informatiques.
Elle est constituée de trois blocs inter-dépendants :
modélisation ce bloc consiste à proposer un modèle mathématique qui devra représenter le plus
fidèlement possible le phénomène réel qui nous intéresse ;
calcul scientifique ce bloc consiste à proposer et à implémenter un algorithme permettant de
calculer une solution approchée du modèle mathématique ;
analyse ce bloc consiste à analyser mathématiquement les propriétés du modèle et de “sa” (ou
“ses”) solution(s).
Afin de mener à bien l’ensemble de ce programme, de nombreuses branches de mathématiques
appliquées ou d’analyse numérique ont vu le jour. Nous pouvons citer par exemple les équations
différentielles ordinaires (ODEs), les équations aux dérivées partielles (PDEs), l’optimisation,
l’analyse du signal, l’algèbre linéaire...
Dans ce chapitre, pour mieux comprendre ce qu’est l’analyse numérique, nous nous intéressons
à un exemple de problème de dynamique de population, c’est-à-dire que nous essayons de déter-
miner un ou plusieurs modèles permettant de décrire l’évolution d’une population d’individus
au cours du temps. Nous ne prenons pas en compte les effets spatiaux pour simplifier : une seule
variable est utilisée pour décrire l’évolution, il s’agit du temps. En particulier, cela signifie que
nos individus ne peuvent pas se déplacer dans l’espace (afin de se regrouper ou de se répartir
sur un territoire par exemple).
L’objectif est évidemment de décrire le ou les modèles mais également d’étudier leurs solutions.
L’écriture du problème mathématique est intéressante en elle-même mais sa ou ses solutions le
sont encore plus pour bien des applications. Nous serons alors amenés à introduire et utiliser de
nombreux concepts et outils mathématiques pour étudier les solutions d’un point de vue théo-
rique (existence, unicité, propriétés de régularité, de périodicité,...) mais aussi d’un point de vue
numérique. En effet, pour de nombreux modèles, il est impossible (c’est un théorème) d’écrire de
manière simple la solution à l’aide de sommes, produits, composées de fonctions, éventuellement
fonctions mathématiques classiques (cos, sin, exp, ln,. . .). Afin de réellement prédire l’évolution
de la population qui nous intéresse, ne restera alors que la simulation numérique, c’est-à-dire le
calcul approché par un ordinateur de la solution exacte.
⁄
1 Nécessité de l’analyse numérique
Supposons que nous souhaitions résoudre un problème général similaire aux modèles continus
de la dynamique des populations :
u0 (t) = f (t, u(t)), t > 0,

(1.1)
u(0) = u0 ,
∗ × R).
où u0 ∈ R et f : (t, x) 7→ f (t, x) ∈ C 1 (R+
Il n’est en général pas possible de donner une expression analytique de la solution t 7→ u(t) (cette
solution est bien définie de manière unique par le théorème de Cauchy-Lipschitz). Il peut alors
être intéressant de déterminer un algorithme permettant de trouver une solution approchée du
problème de Cauchy 1.1. Une méthode pour y parvenir est d’intégrer l’équation entre t = 0 et
t>0: Z t Z t
u(t) = u0 + u0 (s) ds = u0 + f (s, u(s)) ds.
0 0
Nous sommes alors ramener à déterminer une valeur approchée d’une intégrale.
De même, il n’existe pas nécessairement de primitive analytique de la fonction s 7→ f (s, u(s)),
surtout que la fonction u n’est pas connue... Pour construire des méthodes de quadrature, c’est-
à-dire des méthodes numériques d’approximation d’intégrales, il est possible de remplacer la
fonction à intégrer par un polynôme qui “ressemble” à la fonction à intégrer.
L’objectif de ce cours est de construire les outils d’analyse numérique permettant de trouver
des méthodes numériques (et de les analyser) approchant la solution du problème 1.1. Nous
définissons un pas de temps ∆t > 0 et nous construisons une suite (un )n∈N telle que un approche
u(tn ) avec tn = n∆t. L’idée est de remplacer la fonction f par un polynôme Pn et de poser
Z tn+1
n+1 n
u =u + Pn (s) ds.
tn
Nous pouvons alors récrire le problème sous la forme
un+1 = un + φ(un , un+1 ),
où φ est une fonction régulière. Nous devrons donc peut-être résoudre une équation pour dérminer
un+1 en fonction de un .
Nous devons donc
. proposer une méthode pour construire le plynôme Pn à partir de quelques valeurs (un et
éventuellement un+1 ) : nous étudierons la théorie de l’interpolation polynomiale ;
. proposer une méthode de quadrature pour approcher l’intégrale par celle d’un polynôme
interpolateur : nous étudierons la théorie de l’intégration numérique ;
. proposer des méthodes pour résoudre des équations ordinaires pour calculer un+1 en fonc-
tion de un : nous étudierons les méthodes de résolution approchée d’équations ordinaires.
2 Rappels d’analyse
Nous souhaitons à présent rappeler quelques résultats essentiels en analyse numérique autour
des formules de Taylor. En effet, l’analyse locale des fonctions passe le plus souvent par un
⁄
Chapitre 1. Introduction et généralités
développement asymptotique et les formules de Taylor sont très utiles dès que les fonctions sont
suffisamment régulières.
Tous les résultats sont présentés sans démonstration car ce sont des prérequis de ce cours.
2.1 Régularité
Nous commençons par rappeler quelques définitions concernant la régularité des fonctions réelles
à valeurs réelles, ainsi que des théorèmes fondamentaux en analyse.
Définition 1.1 – continuité

Soit f : I → R où I est un intervalle de R. Soit x ∈ I. La fonction f est dite continue en
x si
∀ε > 0 ∃η > 0 : y ∈ I ∩ ]x−η, x+η[ =⇒ |f (x) − f (y)| 6 ε.
De plus, la fonction f est dite continue sur I si elle est continue en tout point x de I.
Théorème 1.2 – Valeurs intermédiaires

Soit f : I → R une fonction continue. Soit (a, b) ∈ I 2 tel que f (a) < 0 < f (b). Alors il
existe c ∈]a, b[ tel que f (c) = 0.
Ce théorème fondamental en analyse peut aussi se résumer par la phrase : « l’image d’un segment
par une fonction continue est aussi un segment ». La démonstration peut se faire de manière
constructive à l’aide de l’algorithme de la dichotomie.
Définition 1.3 – dérivabilité

Soit f : I → R où I est un intervalle de R. La fonction f est dite dérivable en x si la
limite
f (x + h) − f (x)
lim
h→0 h
existe et est finie. Dans ce cas, nous appelerons dérivée en x cette limite et nous la noterons
f 0 (x). De plus, la fonction f est dite dérivable sur I si elle est dérivable en tout point x
de I.
Théorème 1.4 – Rolle

Soit f : [a, b] → R une fonction continue sur [a, b] et dérivable sur ]a, b[. Supposons que
f (a) = f (b). Alors il existe c ∈]a, b[ tel que f 0 (c) = 0.
Théorème 1.5 – Accroissements finis

Soit f : [a, b] → R une fonction continue sur [a, b] et dérivable sur ]a, b[. Alors il existe
c ∈]a, b[ tel que
f (b) − f (a) = f 0 (c)(b − a).
⁄
Définition 1.6
Soit I un intervalle ouvert de R. Nous définissons par récurrence les espaces C k (I), k ∈ N,
n o
C 0 (I) = f : I → R continue sur I ,
n o
C k (I) = f : I → R dérivable sur I telle que f 0 ∈ C k−1 (I) , k > 0.
2.2 Formules de Taylor
Les formules de Taylor sont toujours démontrées à partir du théorème des accroissements finis
et sont basées sur la remarque évidente suivante
Z 1
f (x + h) − f (x) = h f 0 (x + ht) dt.
0
Tous les théorèmes proposés ont la même base mais les hypothèses sont légèrement différentes
(sur la régularité de la fonction) et il est en général utile de réfléchir au choix de le “bonne”
formule selon ce que l’on veut démontrer.
Théorème 1.7 – Taylor avec reste intégral

Soit f : I → R ∈ C n+1 (I). Alors pour tout x ∈ I et h ∈ R tel que x+h ∈ I
n 1
hk
Z
X 1 (n+1)
f (x + h) = f (k) (x) + hn+1 f (x + hs)(1 − s)n ds. (1.2)
k! 0 n!
k=0
Théorème 1.8 – Taylor – Lagrange

Soit f : I → R ∈ C n (I) dérivable jusqu’à l’ordre n+1. Alors pour tout x ∈ I et h ∈ R tel
que x+h ∈ I il existe ξ ∈]0, 1[ tel que
n
X hk hn+1 (n+1)
f (x + h) = f (k) (x) + f (x + hξ). (1.3)
k! (n + 1)!
k=0
Théorème 1.9 – Taylor – Young

Soit f : I → R ∈ C n−1 (I) dérivable jusqu’à l’ordre n. Alors pour tout x ∈ I et h ∈ R tel
que x+h ∈ I
n
X hk
f (x + h) = f (k) (x) + Rn (h), (1.4)
k!
k=0
où la fonction Rn est négligeable devant la fonction h 7→ hn au voisinage de 0. C’est-à-dire
Rn (h)
lim = 0.
h→0,h6=0 hn
2.3 Convexité
Nous rappelons les définitions d’une fonction convexe, strictement convexe, concave et stricte-
ment concave. Ces fonctions sont très importantes en analyse et au moins localement pour les
⁄
Chapitre 1. Introduction et généralités
fonctions régulières, il s’agit de cas générique.
Définition 1.10 – fonction convexe

Soit f : [a, b] → R. La fonction f sera dite convexe si
∀(x, y) ∈ [a, b]2 , x 6= y,

∀λ ∈]0, 1[, f λx + (1 − λ)y 6 λf (x) + (1 − λ)f (y).
La fonction sera dite strictement convexe si l’inégalité est stricte.
Finalement une fonction f sera dite concave si −f est convexe.
Proposition 1.11 – convexité pour les fonctions régulières

Soit f : [a, b] → R de classe C 2 . La fonction f est convexe si, et seulement si, f 00 est
positive.
Démonstration. Nous allons prouver plus exactement l’équivalence entre être convexe et avoir une dé-
rivée croissante pour une fonction f de classe C 1 .
Supposons que f est convexe. Soit x < y. Nous comparons f 0 (x) et f 0 (y) en utilisant la corde entre ces
deux points. Nous pouvons récrire la convexité en prenant
h
λx + (1 − λ)y = x + h =⇒ (1 − λ)(y − x) = h =⇒ λ = 1 − .
y−x
Nous en déduisons que
f (x + h) − f (x) f (y) − f (x)
6 .
h y−x
Nous récrivons ensuite la convexité en prenant
h
λx + (1 − λ)y = y − h =⇒ λ(y − x) = h =⇒ λ = .
y−x
f (y) − f (x) f (y) − f (y − h)
6 .
y−x h
En passant à la limite h → 0, nous obtenons que la fonction f 0 est croissante, donc f 00 est positive.
Réciproquement, si f 00 est positive, la fonction f 0 est croissante. Soit (x, y) ∈ [a, b]2 et soit λ ∈ [0, 1].
Pour fixer les idées, supposons que x 6 y. Posons z = λx + (1 − λ)y. Nous avons donc x 6 z 6 y.
Z y Z 1
f (y) − f (z) = f 0 (t) dt = λ(y − x) f 0 (sy + (1 − s)z) ds,
z 0
Z y Z 1
f (y) − f (x) = f 0 (t) dt = (y − x) f 0 (sy + (1 − s)x) ds.
x 0
Comme f 0 est croissante,
∀s ∈ [0, 1] sy + (1 − s)z > sy + (1 − s)x =⇒ f 0 (sy + (1 − s)z) > f 0 (sy + (1 − s)x).
Nous concluons
f (y) − f (z) > λf (y) − λf (x) ⇐⇒ f (z) 6 λf (x) + (1 − λ)f (y),
ce qui termine la preuve en reprenant l’expression de z = λx + (1 − λ)y. !
Nous pouvons également remarquer que si une fonction est de classe C 2 avec une dérivée seconde
strictement positive, alors la fonction est strictement convexe. La réciproque n’est pas vraie car
une fonction strictement convexe peut avoir une dérivée seconde qui s’annule (par exemple en
un point comme la fonction x 7→ x4 ). Cependant, il est possible de démontrer que la dérivée
⁄
seconde d’une fonction strictement convexe ne peut pas s’annuler sur un intervalle de longueur
non nulle.
⁄
2 Interpolation polynomiale
En analyse numérique, les polynômes de Lagrange, du nom de Joseph-Louis
Lagrange, permettent d’interpoler une série de points par un polynôme qui
passe exactement par ces points appelés aussi nœuds. Cette technique d’in-
terpolation polynomiale a été découverte par Edward Waring en 1779 et
redécouverte plus tard par Leonhard Euler en 1783.
Edward Waring (1736–1798) était un mathématicien britannique. Il est en-
tré au Magdalene College de Cambridge en tant que sizar, et est devenu
Senior wrangler en 1757. Waring a écrit un certain nombre d’articles dans
les Philosophical Transactions of the Royal Society, traitant de la résolution
d’équations algébriques, de théorie des nombres, de séries, de l’approxima- Joseph Louis de Lagrange
tion des racines, de l’interpolation, de la géométrie des sections coniques et
de la dynamique.
Leonhard Euler (1707–1783) est un mathématicien et physicien suisse, qui
passa la plus grande partie de sa vie dans l’Empire russe et en Allemagne.
Il était notamment membre de l’Académie royale des sciences de Prusse à
Berlin. Euler fit d’importantes découvertes dans des domaines aussi variés
que le calcul infinitésimal et la théorie des graphes. Il introduisit également
une grande partie de la terminologie et de la notation des mathématiques
modernes, en particulier pour l’analyse mathématique, comme la notion de
fonction mathématique. Il est aussi connu pour ses travaux en mécanique,
en dynamique des fluides, en optique et en astronomie ou en géométrie du Edward Waring
triangle.
Joseph Louis de Lagrange (1736–1813) né de parents français descendants de
Descartes, est un mathématicien, mécanicien et astronome sarde, naturalisé
français. En mathématiques, il fonde le calcul des variations, avec Euler, et
la théorie des formes quadratiques. En physique, en précisant le principe de
moindre action, avec le calcul des variations, vers 1756, il invente la fonction
de Lagrange, qui vérifie les équations de Lagrange, puis développe la mé-
canique analytique, vers 1788, pour laquelle il introduit les multiplicateurs
de Lagrange. Il entreprend aussi des recherches importantes sur le problème
des trois corps en astronomie, un de ses résultats étant la mise en évidence Leonhard Euler
des points de libration (dits points de Lagrange) (1772).
1 Interpolation
L’interpolation d’une fonction de R dans R est un problème ancien qui a de très nombreuses
applications que ce soit d’un point de vue purement théorique comme l’intégration numérique
(et la résolution d’équations différentielles par conséquent), la construction d’algorithme de ré-
solution de f (x) = 0 ou d’un point de vue industriel comme la visualisation de données, la
représentation 3D d’objet. Essentiellement, l’objectif est de remplacer la fonction compliquée
(ou le nuage de points) qui nous est donnée par une fonction beaucoup plus simple : un po-
lynôme. Malheureusement deux concepts différents apparaissent rapidement : l’interpolation et
l’approximation.
Si nous imaginons un nuage de points (correspondant par exemple à un ensemble discret de
valeurs (xi , f (xi )) pour une certaine fonction f : R → R), il est possible :
. de l’interpoler, c’est-à-dire de faire passer un polynôme par tous les points ;
. de l’approcher (par exemple au sens des moindres carrés), c’est-à-dire de trouver un poly-
nôme de bas degré dont la différence pour une certaine norme est petite.
A la figure Fig. 2.1, nous trouvons des illustrations de l’interpolation et de l’approximation par
⁄
un polynôme de degré 2 lorsque le nombre de points N varie. La fonction f associée à ces nuages
de points est une perturbation aléatoire de la fonction x → 4x(1 − x).
Interpolation ou approximation Interpolation ou approximation

1.2 1.2
1.0 1.0
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0.0 points 0.0 points

interpolation interpolation
0.2 approximation 0.2 approximation
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

1.2 1.2
1.0 1.0
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

1.2 1.2
1.0 1.0
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
Figure 2.1 – Interpolation et approximation pour N ∈ {4, 6, 10, 25, 50, 100}.
Nous pouvons remarquer (sans savoir pour le moment comment l’interpolée a été calculée) que
le polynôme interpolateur passe par les bons points lorsque N est petit (moins que 10 dans
notre exemple) mais qu’il est mal calculé lorsque N devient grand : bien que par définition il
devrait passer par tous les points d’interpolation, le calcul retourne un polynôme qui ne passe
pas par ces points. Ainsi, il faut bien comprendre qu’interpoler ne signifie pas approcher et que
lorsque le nombre de points est élevé, l’interpolation perd de son intérêt. Cependant, cette notion
d’interpolation est très utile en analyse numérique et en calcul scientifique pour approcher (oui
cette fois c’est vrai) localement une fonction à l’aide de quelques points d’interpolation. Nous
nous en servirons pour proposer des méthodes de calcul d’intégrales approchées, des formules de
résolution d’équations différentielles, des méthodes de résolution d’équations non linéaires...
⁄
Chapitre 2. Interpolation polynomiale
Dans ce chapitre, nous nous intéresserons essentiellement à l’interpolation par des polynômes
interpolateurs de Lagrange, mais nous verrons également que ce n’est pas la seule. Nous commen-
çons par nous intéresser à quelques exemples simples avant de généraliser les résultats obtenus.
2 Exemples avec peu de points
Pour bien comprendre la notion de polynôme interpolateur, nous allons commencer par des
exemples pour N petit, où N est le nombre de points par lequel nous imposons que le polynôme
doit passer.
2.1 Interpolation à un seul point
Prenons dans cette partie N = 1, c’est-à-dire que nous cherchons à faire passer un polynôme
par un point de coordonnées (x1 , y1 ). Il est évident qu’il en existe une infinité si le degré du
polynôme n’est pas limité. Nous pouvons donc chercher un ou le polynôme de degré le plus petit
possible passant par ce point. Une illustration est donnée à la figure 2.2.
Figure 2.2 – Exemple de polynôme interpolateur de Lagrange lorsque N = 1.
Dans ce cas, nous cherchons le polynôme de degré au plus 0 qui passe par le point de coordonnées
(x1 , y1 ). C’est donc le polynôme constant P = y1 . Il est évidemment unique.
2.2 Interpolation à deux points : exemple de la droite
Nous nous intéressons à présent à l’exemple le plus simple dont les applications sont malgré tout
innombrables : par deux points passe une unique droite. Ce résultat ou plutôt cette définition
selon Archimède doit évidemment être suivi d’un calcul analytique afin de déterminer l’équation
de cette droite. Dans cette section, nous prenons donc N = 2.
⁄
2.2.1 Un peu de géométrie euclidienne
Théorème 2.1 – Polynôme interpolateur de degré 1

Soit (x1 , y1 ) et (x2 , y2 ) deux points distincts de R2 . Il existe une unique droite D passant
par ces deux points a :
(x, y) ∈ D ⇐⇒ (x − x1 )(y2 − y1 ) − (y − y1 )(x2 − x1 ) = 0.
Si de plus x1 6= x2 , il existe un unique polynôme P ∈ R1 [X] tel que
(X − x1 )y2 − (X − x2 )y1
D = {(x, P (x)), x ∈ R} avec P = .
x2 − x1
a. On distinguera bien D qui est un object géométrique c’est-à-dire une partie du plan, de P le polynôme
(ou parfois fonction polynomiale) dont le graphe est égal à D.
Démonstration. Nommons les points afin de faire de la géométrie euclidienne ! Notons M = (x, y),
M1 = (x1 , y1 ), M2 = (x2 , y2 ), les trois points de R2 qui nous intéressent. Le point M est sur la droite
# » # »
D = (M1 M2 ) si, et seulement si, les vecteurs M1 M et M1 M2 doivent être colinéaires, ce qui peut se
récrire sous la forme
x − x1 x2 − x1
= 0.
y − y1 y2 − y1
Le calcul de ce déterminant donne directement l’expression de la droite D. Puis, dans le cas x1 6= x2 , il

est possible de paramétrer cette droite par l’abscisse x en divisant par x2 − x1 . Cela induit l’expression
du polynôme P . !
Remarquez l’expression proposée pour le polynôme interpolateur P qui est symétrique par rap-
port aux deux points. Nous pouvons en proposer deux autres : la première est la décomposition
de ce polynôme dans la base canonique (1, X)
y2 − y1 x2 y1 − x1 y2
P =X + ,
x2 − x1 x2 − x1
la seconde dans une base “adaptée” aux points d’interpolation (1, X − x1 )
y2
P = (X − x1 ) + y1 .
x2 − x1
La question que nous devons nous poser à présent est la suivante : comment calculer le polynôme
P (qui est seulement de degré 1 dans notre exemple) et surtout comment calculer les valeurs de
P (x) pour de nombreuses valeurs de x. En effet, il est parfois nécessaire de connaitre l’expression
exacte du polynôme (de chacun de ses coefficients dans une certaine base), parfois plutôt de
calculer l’évaluation de ce polynôme un grand nombre de fois.
La réponse à cette question ne peut pas être l’expression trouvée dans ce théorème puisque la
démonstration géométrique pour la trouver ne fonctionnera plus lorsque le nombre de points
d’interpolation sera supérieur à 2.
2.2.2 Une première méthode de calcul
La première méthode qui vient à l’esprit est de décomposer le polynôme recherché dans la base
canonique et de résoudre les équations associées. Nous cherchons donc à déterminer deux réels
⁄
a0 et a1 , ces deux réels étant les coefficients dans la base canonique du polynôme P , c’est-à-dire
tels que P = a0 + a1 X. Nous écrivons ensuite les relations qui doivent être satisfaites :
(
a0 + a1 x1 = y1

1 x1 a0 y1
⇐⇒ = .
a0 + a1 x2 = y2 1 x 2 a 1 y2
Comme x1 6= x2 , la matrice est inversible et nous pouvons laisser le soin à Python et son
package numpy.linalg de résoudre le problème : les coefficients a0 et a1 existent donc bien et
sont uniques.
Allons un peu plus loin dans l’analyse de cette méthode. Nous posons

1 x1
A= .
1 x2
Nous avons donc

1 x2 −x1
A−1 = .
x2 − x1 −1 1
Ainsi, si x1 et x2 sont proches, le calcul de la matrice A−1 et plus précisément la résolution

du système linéaire conduira à des erreurs d’arrondis. D’une manière générale, ce phénomène
s’amplifiera lorsque le nombre de points d’interpolation augmente : c’est ce que nous voyons à
la Fig. 2.1 où le polynôme interpolateur “numérique” ne passe plus par les points d’interpolation
lorsque le nombre de points est trop grand.
2.2.3 Expression du polynôme dans une base adaptée
Nous proposons une seconde méthode pour évaluer le polynôme interpolateur. L’idée est de le
représenter dans une autre base que la base canonique : une base adaptée aux points d’interpo-
lation. Nous cherchons donc à déterminer deux réels a0 et a1 tels que P = a0 + a1 (X − x1 ). Le
calcul se fait alors plus simplement :
y2 − y1
P (x1 ) = y1 =⇒ a0 = y1 , P (x2 ) = y2 =⇒ a1 = .
x2 − x1
Par ailleurs, nous remarquons que les calculs du polynôme interpolateur passant par les points
(x1 , y1 ) et (x2 , y2 ), noté P1,2 ici, peuvent être réutilisés si l’on souhaite ajouter un point d’inter-
polation.
2.3 Interpolation à trois points
Prenons dans cette partie N = 3. Une illustration est donnée à la figure 2.3.
Dans ce cas, nous cherchons le polynôme de degré au plus 2 qui passe par les trois points de
coordonnées (x1 , y1 ), (x2 , y2 ) et (x3 , y3 ). C’est donc une parabole que l’on peut chercher sous la
forme P = aX 2 + bX + c. Les calculs commencent à être plus compliqués et il est nécessaire
d’avoir une stratégie pour trouver les coefficients a, b et c. Ces trois coefficients sont solutions
du système linéaire suivant  2
P (x1 ) = ax1 + bx1 + c = y1 ,

P (x2 ) = ax22 + bx2 + c = y2 ,

P (x3 ) = ax23 + bx3 + c = y3 .

⁄
Figure 2.3 – Exemple de polynôme interpolateur de Lagrange lorsque N = 3.
Un calcul un peu fastidieux permet de déterminer la valeur des trois coefficients a, b et c :

y1 (x2 − x3 ) + y2 (x3 − x1 ) + y3 (x1 − x2 )
a=− ,
(x2 − x3 )(x3 − x1 )(x1 − x2 )
y1 (x22 − x23 ) + y2 (x23 − x21 ) + y3 (x21 − x22 )
b= ,
(x2 − x3 )(x3 − x1 )(x1 − x2 )
y1 x2 x3 (x2 − x3 ) + y2 x3 x1 (x3 − x1 ) + y3 x1 x2 (x1 − x2 )
c=− .
(x2 − x3 )(x3 − x1 )(x1 − x2 )
Plusieurs remarques peuvent alors être faites :

. les calculs deviennent rapidement lourds et il faut trouver une méthode générale
Q ;
. il semble que les coefficients obtenus pour N points ont un dénominateur i<j (xi − xj )
(le même que celui des polynômes de Lagrange).
Exercice 2.2 – Polynôme interpolateur de degré 2

Déterminez l’expression du polynôme interpolateur passant par les points (x1 , y1 ), (x2 , y2 )
et (x3 , y3 ) dans la base adaptée (1, (X − x1 ), (X − x1 )(X − x2 )).
Solution. Nous écrivons P = a0 + a1 (X − x1 ) + a2 (X − x1 )(X − x2 ). Nous avons en évaluant en x1

P (x1 ) = a0 + a1 (x1 − x1 ) + a2 (x1 − x1 )(x1 − x2 ) = y1 =⇒ a0 = y1 .
Puis en évaluant en x2 :
y2 − y1
P (x2 ) = y1 + a1 (x2 − x1 ) + a2 (x2 − x1 )(x2 − x2 ) = y2 =⇒ a1 = .
x2 − x1
Enfin en évaluant en x3 :
y2 − y1
P (x3 ) = y1 + (x3 − x1 ) + a2 (x3 − x1 )(x3 − x2 ) = y3
x2 − x1
=⇒
y3 − y1 y2 − y1
y3 − y1 − xy22 −y1
(x − x ) −
−x1 3 1 x − x1 x2 − x1
a2 = = 3
(x3 − x1 )(x3 − x2 ) x3 − x2
⁄
3 Polynôme interpolateur de Lagrange
Dans cette section, nous étudions les propriétés du polynôme interpolateur de Lagrange et nous
décrivons quelques algorithmes permettant de le calculer.
3.1 Définition et propriétés
Etant donnés N points du plan (x1 , y1 ), . . . , (xN , yN ), avec N ∈ N∗ , nous souhaitons construire
un polynôme de degré minimal qui passe par tous ces points.
Le premier résultat fondamental est que le degré minimal est N −1, qu’il est toujours possible
de construire un tel polynôme et qu’enfin ce polynôme est unique.
Théorème 2.3 – Isomorphisme d’espaces vectoriels

Soit (x1 , . . . , xN ) ∈ RN tels que xi 6= xj si i 6= j. L’application
Φ : RN −1 [X] −→ RN
P 7−→ (P (x1 ), . . . , P (xN ))
est bijective (c’est un isomorphisme).
N
Démonstration. Pour montrer que l’application linéaire Φ : RN −1 [X] → R définie par Φ(P ) =
P (x1 ), . . . , P (xN ) est un isomorphisme d’espaces vectoriels, il suffit de vérifier qu’elle est injective,
puisque les espaces sont de même dimension N . Supposons donc que P ∈ RN −1 [X] tel que Φ(P ) = 0.
Nous avons donc un polynôme de degré au plus N −1 qui a N racines, c’est le polynôme nul. !
Définition 2.4 – Polynôme interpolateur de Lagrange

Etant donnés des réels x1 , . . . , xN deux à deux distincts et des réels quelconques y1 , . . . , yN ,
l’unique polynôme P qui vérifie
P ∈ RN −1 [X] et P (xk ) = yk pour 1 6 k 6 N
est appelé le polynôme interpolateur de Lagrange aux points (x1 , y1 ), . . . , (xN , yN ).
Un cas particulièrement intéressant lorsque l’on fait de l’analyse est le cas où le nuage de points
(x1 , y1 ), . . . , (xN , yN ) est sur le graphe d’une fonction régulière f . C’est-à-dire que f (xi ) = yi ,
1 6 i 6 N . Dans ce cas, on dira que le polynôme P est le polynôme interpolateur de la fonction
f aux points x1 , . . . , xN .
Définition 2.5 – Polynôme interpolateur de Lagrange

Si f est une fonction continue sur un intervalle [a, b] et à valeurs dans R et si x1 , . . . , xN
sont N points distincts de [a, b], alors l’unique polynôme P ∈ RN −1 [X] tel que P (xi ) =
f (xi ) pour 1 6 i 6 N , est appelé polynôme d’interpolation de Lagrange de f aux points
x1 , . . . , xN .
Il est possible d’estimer l’écart entre une fonction f et son polynôme interpolateur de Lagrange
lorsque la fonction est régulière. Cette estimation repose sur le lemme suivant.
⁄
Lemme 2.6 – Rolle généralisé

Soient a, b deux réels tels que a < b et g : [a, b] → R une fonction continue. Si g s’annule en
m (m > 2) points t1 , . . . , tm deux à deux distincts de [a, b] et si g est m − 1 fois dérivable
sur ]a, b[, alors il existe t ∈]t1 , tm [ tel que g (m−1) (t) = 0.
Démonstration. La preuve se fait par récurrence sur m > 2. Pour m = 2, c’est exactement le théorème
de Rolle (Th. 1.4). Supposons le résultat du lemme vrai pour m > 2 et prenons g une fonction m fois
dérivable qui s’annule en t1 , . . . , tm+1 points distincts. D’après le théorème de Rolle 1.4, la fonction g 0
s’annule sur chaque intervalle ]ti , ti+1 [, 1 6 i 6 m, en un point noté t̃i . La fonction g 0 est donc m−1
fois dérivable et s’annule m fois, par hypothèse de récurrence, il existe t ∈]t̃1 , t̃m [⊂]t1 , tm+1 [ tel que
g (m) (t) = 0. !
Théorème 2.7 – Théorème de l’erreur

Soient a, b deux réels tels que a < b et f : [a, b] → R une fonction continue. Etant donné
N points x1 , x2 , . . . , xN deux à deux distincts de [a, b], on pose
ωN (X) = (X − x1 ) . . . (X − xN )
et on note P ∈ RN −1 [X] le polynôme d’interpolation de Lagrange de f aux points

x1 , . . . , xN .
Si f est N fois dérivable sur ]a, b[, alors pour tout x ∈ [a, b] il existe t ∈]a, b[ tel que
f (N ) (t)
f (x) − P (x) = ωN (x) .
N!
Démonstration. Pour x fixé dans [a, b]\{x1 , . . . , xN }, on applique le lemme de Rolle généralisé (Lm. 2.6)
à la fonction auxiliaire
f (x) − P (x)
F (t) = f (t) − P (t) − ωN (t).
ωN (x)
La fonction F est N fois dérivable et s’annule aux N +1 points x, x1 , . . . , xN . D’après le lemme de Rolle
généralisé (Lm. 2.6), il existe t ∈]a, b[ tel que F (N ) (t) = 0. Calculons la dérivée N ième de F .
f (x) − P (x) f (N ) (t)

F (N ) (t) = f (N ) (t) − 0 − N! =⇒ f (x) − P (x) = ωN (x) .
ωN (x) N!
Par ailleurs, si x ∈ {x1 , . . . , xN }, nous avons f (x)−P (x) = 0 et ωN (x) = 0, ce qui termine la preuve. !
Corollaire 2.8
Si f (N ) est uniformément bornée par M sur ]a, b[, alors pour tout x ∈ [a, b]
M M
|f (x) − P (x)| 6 |ωN (x)| 6 (b − a)N .
N! N!
3.2 Expression dans la base canonique
Soit N ∈ N, N > 1. Etant donné N réels deux à deux distincts x1 , . . . , xN et N réels quelconques
y1 , . . . , yN , on note P le polynôme d’interpolation de Lagrange aux points (x1 , y1 ), . . . , (xN , yN ).
On souhaite évaluer P (x) pour x ∈ R\{x1 , . . . , xN } quelconque.
La première idée pour évaluer le polynôme interpolateur est de résoudre un système linéaire
satisfait par les coordonnées de ce polynôme dans la base canonique. Cette méthode conduit à
⁄
résoudre un système de Vandermonde. Si nous notons
N
X −1
P = ai X i ,
i=0
l’expression du polynôme interpolateur dans la base canonique, nous avons
N
X −1
P (xk ) = ai xik = yk
i=0
ce qui se récrit sous forme matricielle V (x1 , . . . , xN )a = y avec
x01 . . . x1N −1
     
a0 y1
V (x1 , . . . , xN ) =  ... ..  , a =  ...  , y =  ...  .
    
. 
N −1
x0N . . . xN aN −1 yN
Cette méthode a deux défauts : la matrice est pleine (donc la résolution du système linéaire
est coûteuse) et elle est mal conditionnée (le conditionnement croı̂t comme l’exponentielle de
N ). Nous envisagerons donc d’autres méthodes de calcul. Cependant, lorsque l’expression du
polynôme interpolateur est donnée dans la base canonique, il est bon d’avoir un algorithme
rapide d’évaluation de ce polynôme en un point x. L’algorithme de Horner permet de n’effectuer
que N additions et N multiplications pour réaliser ce calcul.
Algorithme 2.1 – Horner

Si x est un réel et si Q est le polynôme défini par
Q(X) = a0 X n + a1 X n−1 + . . . + an−1 X + an ,
la suite finie q0 , q1 , . . . , qn définie par q0 = a0 et qk = qk−1 x + ak pour k = 1, . . . , n vérifie

qn = Q(x).
def Horner (p , x ):
y = 0
for ak in p :
y *= x
y += ak
return y
h i
Démonstration. Il suffit d’écrire Q(X) = an + X an−1 + X an−2 + X . . . a1 + X(a0 ) . !
La classe poly1d du package numpy de python est une classe pour les polynômes d’une variable.
Un polynôme est stocké sous la forme d’un tableau numpy contenant les coefficients dans la base
canonique dans l’ordre décroissant. Il est possible d’évaluer ce polynôme en un point (ou même
en plusieurs points directement) en utilisant l’algorithme de Horner. Par exemple
import numpy as np
p = np . poly1d ([1 , 2 , 3]) # polynomial X ^2 + 2 X + 3
p (2) # p (2) = 11
Ainsi, nous pouvons proposer un premier algorithme de calcul de l’interpolation aux points
(xi , yi ), 1 6 i 6 N .
⁄
Algorithme 2.2 – InterpVdM
import numpy as np
def InterpVdM (x , y , xx ):
"""
Interpolation with the Vandermonde method
"""
x = np . asanyarray (x , dtype = ’ float ’)
y = np . asanyarray (y , dtype = ’ float ’)
xx = np . asanyarray ( xx , dtype = ’ float ’)
n = x . size
if y . size != n :
print ( ’ Error in InterpVdm : x and y do not have the same size ’)
A = np . ones (( n , n ))
for k in range (1 , n ):
A [: , k ] = A [: , k -1] * x
p = np . linalg . solve (A , y )
yy = Horner ( p [:: -1] , xx )
return yy if xx . size != 1 else np . asscalar ( yy )
La fonction InterpVdM prend trois arguments de type numpy array : x, y et xx et calcule

le polynôme interpolateur aux points (xi , yi )16i6N puis l’évalue aux points xx en utilisant
l’algorithme de Horner.
3.3 Formule de Lagrange et poids barycentriques
Etant donnés N réels (x1 , . . . , xN ) deux à deux distincts, l’application
Φ : RN −1 [X] −→ RN
P 7−→ (P (x1 ), . . . , P (xN ))
est un isomorphisme d’espaces vectoriels. Nous pouvons alors définir l’image de la base canonique
de RN par l’application Φ−1 et nous la noterons (L1 , . . . , LN ). A la figure 2.4, vous trouverez les
premiers polynômes de Lagrange pour des points équi-répartis.
Figure 2.4 – Polynôme de Lagrange pour des points équi-répartis.
⁄
Définition 2.9 – Polynômes de Lagrange

Le i-ème polynôme de Lagrange, noté Li ∈ RN −1 [X], pour 1 6 i 6 N est défini par
Y X − xj
Li = , 1 6 i 6 N.
xi − xj
j6=i
Il est tel que (

1 si i = j,
Li (xj ) =
0 sinon.
En utilisant cette base adaptée aux points (x1 , . . . , xn ), nous pouvons donner une nouvelle ex-
pression du polynôme interpolateur de Lagrange. Nous avons de manière évidente
n
X
P (X) = yi Li (X).
i=1
Une récriture astucieuse de cette dernière expression conduit à la formule de Lagrange, très
efficace pour évaluer ce polynôme. A nouveau on pose
ωN (X) = (X − x1 ) . . . (X − xN ).
Proposition 2.10 – Formule de Lagrange

On a une formule dite formule de Lagrange pour P (x)
N
P yi
N (x−xi )ωN0 (x )
X ωN (x) i=1
i
P (x) = yi 0 (x ) ou encore P (x) =
(x − xi )ωN i PN
1
i=1 0 (x )
(x−xi )ωN i
i=1
valables lorsque x est un réel distinct de x1 , . . . , xN .
Démonstration. Les polynômes L1 , . . . , LN vérifient
Li (xi ) = 1 et Li (xj ) = 0 pour tout j ∈ {1, . . . , N }\{i}.

N
P
Ces polynômes forment une base de RN −1 [X] et vérifient Li (X) = 1.
i=1
Pour x ∈ R\{x1 , . . . , xN }, les valeurs Li (x) sont données par
Y x − xj ωN (x) Y
0
Li (x) = = 0 (x ) car ωN (xi ) = (xi − xj ).
xi − xj (x − xi )ωN i
j6=i j6=i
N
P
Dans la base L1 , . . . , LN , le polynôme P s’écrit P (X) = yi Li (X). On obtient donc pour P (x) la
i=1
première expression de la proposition valable lorsque x est un réel distinct de x1 , . . . , xN .
N
P
Pour la seconde, nous utilisons que Li (X) = 1. Nous avons donc
i=1
N
X ωN (x) 1
0 (x ) = 1
(x − xi )ωN
=⇒ ωN (x) = N
,
i=1 i P 1
0 (x )
(x−xi )ωN i
i=1
ce qui termine la preuve. !
⁄
Cette formule de Lagrange permet d’évaluer une valeur P (x) en un point x en O(2N 2 ) itérations
pour le premier x, puis O(2N ) itérations pour les suivants. En effet, le calcul de ωN 0 (x ) =
Q i
j6=i (x i − x j ) pour chaque valeur de i prend 2N − 2 opérations. On peut d’ailleurs stocker ces
valeurs si l’on doit calculer une autre valeur de P (x) puisque cette valeur ne dépend pas de x.
Pour compléter le calcul de P (x), il ne reste plus qu’à faire 2 multiplications pour chaque valeur
de i et sommer les différents termes.
Cependant, pour un traitement vectoriel (et en python, c’est le plus efficace) de cette formule,
il est nécessaire de traiter séparément les points d’interpolation puisque la formule contient des
divisions par 0 dans ce cas.
Algorithme 2.3 – InterpLagrange
import numpy as np
import itertools as it
def InterpLagrange (x , y , xx ):
"""
Interpolation using the Lagrange formula
"""
x = np . asanyarray (x , dtype = ’ float ’)
y = np . asanyarray (y , dtype = ’ float ’)
xx = np . asanyarray ( xx , dtype = ’ float ’)
n = x . size
if y . size != n :
print ( " Error in InterpLagrange : x and y don ’t have the same size " )
w = np . ones (( n ,))
for i in range ( n ):
for j in it . chain ( range ( i ) , range ( i +1 , n )):
w [ i ] *= ( x [ i ] - x [ j ])
weight = 1./ w
if xx . size == 1:
N , D = 0. , 0.
if xx == x [ i ]:
return y [ i ]
dxi = weight [ i ] / ( xx - x [ i ])
N += y [ i ] * dxi
D += dxi
else :
N , D = np . zeros ( xx . shape ) , np . zeros ( xx . shape )
ind = []
ind . append ( np . where ( xx == x [ i ])) # indices where P ( xi ) = yi
indloc = np . where ( xx != x [ i ])
dxi = weight [ i ] / ( xx [ indloc ] - x [ i ]) # avoid divide by 0
N [ indloc ] += y [ i ] * dxi
D [ indloc ] += dxi
for i in range ( n ): # fix yi in xi
N [ ind [ i ]] , D [ ind [ i ]] = y [ i ] , 1.
return N / D
La fonction InterpLagrange prend trois arguments de type numpy array : x, y et xx et

calcule le polynôme interpolateur aux points (xi , yi )16i6N puis l’évalue aux points xx en
utilisant la formule de Lagrange.
⁄
3.4 Comportement asymptotique lorsque n tend vers l’infini
On peut montrer que si f est une fonction définie par une série entière de rayon de convergence
infini (par exemple f (x) = exp(x) ou f (x) = cos(x)), le polynôme d’interpolation de Lagrange
de f aux points x1 , x2 , . . . , xN de l’intervalle [a, b] converge uniformément vers f sur [a, b] lorsque
N tend vers l’infini, quelque soit la répartition des points d’interpolation.
Lorsque les points x1 , . . . , xN sont équirépartis sur [a, b] (i.e. pour tout j ∈ [[1, N ]], xj = a +
(j − 1)/h où h = (b − a)/(N − 1)), on pourrait s’attendre à ce que le polynôme d’interpolation
de Lagrange d’une fonction f définie sur [a, b] converge uniformément vers f sur [a, b] quand N
tend vers l’infini, au moins dans le cas d’une fonction assez régulière. Or dans chacun des cas
suivants (et ce ne sont que des exemples)
f : [0, 1] −→ R f : [−1, 1] −→ R f : [−1, 1] −→ R

√
x 7−→ x x 7−→ |x| x 7−→ 4x21+1
si x1 , . . . , xN sont équirépartis, on a max |P (t) − f (t)| → +∞ lorsque N tend vers l’infini. Ce

comportement est appelé phénomène de Runge.
Un moyen de limiter ce phénomène est d’utiliser d’autres points d’interpolation. Si f est lip-
schitzienne sur [−1, 1] ou plus généralement si f est höldérienne sur [−1, 1] (i.e. il existe α ∈]0, 1]
et C > 0 tels que |f (x) − f (y)| 6 C|x − y|α pour tout x, y dans [−1, 1]) et si on prend
pour points d’interpolation les zéros x1 , . . . , xN du N -ième polynôme de Tchebychev, alors
max |P (t) − f (t)| → 0 lorsque N tend vers l’infini.
Définition 2.11 – Polynômes de Tchebychev

Les polynômes de Tchebychev sont définis par récurrence
T0 = 1, T1 = X, Tn = 2XTn−1 − Tn−2 , n > 2.
Propriété 2.12
Le nième polynôme de Tchebychev Tn vérifie les assertions suivantes :
1. Tn est de degré exactement n et son coefficient de plus haut degré vaut 2n−1 pour
n > 1.
2. Tn est scindé à racines simples.
(2j − 1)π
Tn (x) = 0 ⇐⇒ x ∈ {x1 , . . . , xn }, xj = cos , 1 6 j 6 n.
2n
3. |Tn (x)| 6 1 pour tout x ∈ [−1, 1]. De plus
kπ
|Tn (x)| = 1 ⇐⇒ x ∈ {x00 , . . . , x0n }, x0k = cos , 0 6 k 6 n.
n
Démonstration. Nous démontrons le premier point par récurrence. Nous définissons (pn ) la propriété
suivante : « Tn est de degré n et son coefficient de plus haut degré vaut 2n−1 pour n > 1 ». Les propriétés
(p0 ) et (p1 ) sont vraies. Supposons que les propriétés (pk ) sont vraies pour 0 6 k 6 n et montrons que
(pn+1 ) est vraie. Nous avons Tn+1 = 2XTn − Tn−1 . Ainsi (pn+1 ) est vraie.
L’essentiel de la suite de la preuve consiste à remarquer que, pour x ∈ [−1, 1], la fonction polynomiale
Tn (x) coı̈ncide avec une fonction trigonométrique cos(n arccos(x)). Cette propriété est vraie pour n = 0
et n = 1 puisque cos(0) = 1 et cos(arccos(x)) = x si x ∈ [−1, 1]. Puis nous utilisons la relation de
⁄
récurrence :
cos((n + 1) arccos(x)) = cos(n arccos(x)) cos(arccos(x)) − sin(n arccos(x)) sin(arccos(x)),

cos((n − 1) arccos(x)) = cos(n arccos(x)) cos(arccos(x)) + sin(n arccos(x)) sin(arccos(x)),
cos((n + 1) arccos(x)) = 2 cos(n arccos(x)) cos(arccos(x)) − cos((n − 1) arccos(x)),
= 2x cos(n arccos(x)) − cos((n − 1) arccos(x)).
Ainsi, la suite de fonctions x →

7 cos(n arccos(x)) vérifie la même relation de récurrence que Tn (x), les
deux fonctions coı̈ncident sur [−1, 1]. Ainsi, pour n > 1 et pour x ∈ [−1, 1],
π
Tn (x) = 0 ⇐⇒ n arccos(x) = 2 mod π,
π
⇐⇒ arccos(x) = mod nπ ,
2n
n o
⇐⇒ x ∈ cos (2j−1)π

2n , j ∈ [[1, n]] .
En utilisant les propriétés de la fonction cos, nous obtenons immédiatement que |Tn (x)| 6 1 pour tout
x ∈ [−1, 1]. De plus
|Tn (x)| = 1 ⇐⇒ n arccos(x) = 0 mod π,

⇐⇒ arccos(x) = 0 mod nπ ,
⇐⇒ x ∈ cos kπ

n , k ∈ [[0, n]] ,
Proposition 2.13
Si Qn est un polynôme de degré n et de même coefficient de plus haut degré que celui de
Tn , alors
max |Qn (x)| > max |Tn (x)| = 1.

x∈[−1,1] x∈[−1,1]
Démonstration. Soit Q un polynôme de degré n > 1 dont le coefficient de degré n vaut 2n−1 (le cas
n = 0 est trivial). Supposons que
max |Q(x)| < max |Tn (x)| = 1.

x∈[−1,1] x∈[−1,1]
Nous avons Tn (x0k ) = (−1)k pour tout k ∈ [[0, n]]. Le polynôme Qn − Tn admet donc au moins une racine
dans l’intervalle ]x0k−1 , x0k [ pour chaque k ∈ [[1, n]]. Mais c’est impossible car Qn − Tn est un polynôme
non nul de degré 6 n − 1. !
Corollaire 2.14
Si ξ1 , . . . , ξN sont N points deux à deux distincts de [−1, 1], on a
N N
Y Y 1 1
max (x − ξj ) > max (x − xj ) = max |TN (x)| = .
x∈[−1,1] x∈[−1,1] x∈[−1,1] 2N −1 2N −1
j=1 j=1
Revenons à l’étude des polynômes interpolateurs de Lagrange. Considérons f une fonction suf-
fisamment régulière et définissons P son polynôme interpolateur en N points de [−1, 1] notés
ξ1 , . . . , ξN . D’après le Corollaire 2.8, l’erreur d’interpolation vérifie
N
kf (N ) k∞ kf (N ) k∞ Y
|f (x) − P (x)| 6 |ωN (x)| 6 kωN k∞ , ωN (x) = (x − ξj ),
N! N!
j=1
où la norme infinie est à prendre sur l’intervalle [−1, 1]. Ainsi, le choix des points ξi qui fournit
⁄
une valeur de kωN k∞ la plus petite possible (on ne peut pas descendre sous 2−N +1 ) est celui
des zéros du polynôme de Tchebychev.
Evidemment, le cas particulier de l’intervalle [−1, 1] est facilement généralisable. En effet, si
l’intervalle d’interpolation vaut [a, b], les points d’interpolation associés sont
a+b b−a (2j − 1)π

xj = + cos , 1 6 j 6 N.
2 2 2N
points équidistants points de Tchebychev
a| | | | | | | | | | |b a|| | | | | | | | ||b

b−a a+b b−a 2i − 1
xi,N = (i − 1) xi,N = + cos π 16i6N
N −1 2 2 2N
On observe que, si la convergence semble bien fonctionner au centre de l’intervalle dans le cas
des points équirépartis, de fortes oscillations de plus en plus violentes apparaissent au bord du
domaine : il n’y a pas de convergence uniforme de PN vers f dans ce cas. En revanche, avec les
points de Tchebychev, la convergence semble bien uniforme. L’évolution de la norme infinie de
l’erreur lorsque N augmente est donnée à la figure 2.5.
⁄
Figure 2.5 – Illustration du phénomène de Runge pour x 7→ 1/(1 + x2 ).
⁄
3 Intégration numérique
Étant donnée une fonction f : [a, b] → R, le calcul de son intégrale
Z b
I(f ) = f (x) dx
a
est un problème qui est souvent compliqué. Dans le cas où l’expression d’une primitive F de f
est connue, le problème est réglé puisque I(f ) = F (b) − F (a). Malheureusement, cela n’est pas
toujours possible :
. F peut ne pas avoir d’expression analytique, par exemple dans les cas suivants
√ 1
x 7→ exp(−x2 ), x 7→ 1 + 2 cos x, x 7→ ;
log x
. f peut n’être connue que par ses valeurs f (x) en certains points x (issues de données
expérimentales par exemple) ;
. f peut n’être connue que comme valeur retournée par un algorithme coûteux permettant
de la calculer pour tout x.
Dans tous ces cas-là, il est nécessaire de calculer une approximation numérique Iapp de I(f ), ce
à quoi on va s’intéresser dans tout ce chapitre.
1 Formules de quadratures
Rb
Comment calculer une valeur approchée de I(f ) = a f (x) dx ? Un résultat théorique donne
une piste : on peut recourir à des sommes de Riemann. Si Σ est une subdivision a = x0 < x1 <
x2 < . . . xn = b pointée en chaque sous-intervalle par xj+1/2 ∈ [xj , xj+1 ] pour tout i ∈ [[0, n − 1]],
on note
n−1
X
S(f, Σ) = f (xj+1/2 )(xj+1 − xj )
j=0
la somme de Riemann de f associée à Σ, et δ(Σ) = max |xj+1 − xj | le pas maximal de la

06j6n−1
subdivision.
Théorème 3.1 – Sommes de Riemann

Soit f ∈ C 0 ([a, b]). Étant donnée une suite de subdivisions Σn de [a, b] dont le pas maximal
δ(Σn ) tend vers 0 lorsque n → ∞, on a
Z b
S(f, Σn ) −−−−−→ f (x) dx.
n→+∞ a
Remarquons tout d’abord que le calcul de In (f ) = S(f, Σn ) n’utilise que des évaluations de f en
certains points xi+1/2 et des opérations arithmétiques standards (sommes et produits) : on peut
donc effectuer ce calcul “à la main” ou sur machine.
Néanmoins, ce théorème ne quantifie pas à quelle vitesse ces sommes de Riemann convergent
vers l’intégrale de f . En réalité, si l’on prend une subdivision Σn régulièrement espacée en n
⁄
sous-intervalles, ce qui correspond à xj = a + jh où h = (b − a)/n et xj+1/2 = xj par exemple, et

si f est de classe C 1 , on ne peut s’attendre en général qu’à une convergence en O(1/n), au sens
où
C
I(f ) − In (f ) 6
n
pour une certaine constante C > 0. C’est-à-dire que pour obtenir une précision à 10−6 de
I(f ) il faudrait évaluer f (x) de l’ordre de 1 million de fois, et effectuer autant d’opérations
arithmétiques : c’est très coûteux et surtout, il risque d’y avoir une forte accumulation d’erreurs
d’arrondis. La question est donc la suivante : peut-on trouver d’autres procédés plus efficaces
pour approcher l’intégrale ?
1.1 Définitions et premières propriétés
On va s’intéresser à des approximations d’une intégrale par des formules faisant intervenir uni-
quement des combinaisons linéaires de valeurs prises par f en certains points. Nous introduisons
une nouvelle notation pour l’intégrale que nous souhaitons calculer afin d’éviter les confusions
dans la suite de ce cours. Dans ce paragraphe, nous nous intéressons donc à évaluer
Z β
I= ϕ(ξ) dξ,
α
pour ϕ : [α, β] → R une fonction continue.
Définition 3.2 – Formule de quadrature

Étant donnés N points ξ1 , . . . , ξN dans un intervalle [α, β] et N poids ω1 , . . . , ωN ∈ R asso-
ciés à chaque point, on appelle formule de quadrature associée aux (ξi )16i6N et (ωi )16i6N
l’application linéaire Ie : C 0 ([α, β]) → R définie par
N
X
I(ϕ)
e = ωi ϕ(ξi ).
i=1
Remarque 3.3
En général, une formule de quadrature est donnée sous la forme
Z β N
X
ϕ(ξ) dξ ≈ ωi ϕ(ξi ),
α i=1
points ξ1 ... ξN
ou est décrite par un tableau .
poids ω1 ... ωN
On sait qu’au voisinage d’un point, une fonction de classe C p−1 s’approche par un polynôme de
degré p − 1. Ainsi, si une formule de quadrature approche correctement l’intégrale de certains
polynômes (par exemple si elle lui est égale), elle devrait mieux approcher l’intégrale de ϕ. Cette
remarque motive la définition suivante d’ordre, et on montrera plus loin dans le cours qu’elle
mesure bien la qualité d’approximation de la méthode.
⁄
Chapitre 3. Intégration numérique
Définition 3.4 – Ordre d’une méthode

On dit qu’une formule de quadrature Ie est d’ordre p si elle est exacte pour tout polynôme
de degré inférieur ou égal à p − 1 (et pas plus), c’est-à-dire si
∀P ∈ Rp−1 [X] I(P

e ) = I(P ),
∃Q ∈ Rp [X] : I(Q)
e 6= I(Q).
Remarquons que, comme les applications Ie et I sont linéaires, il est équivalent de dire : la formule
de quadrature Ie est d’ordre p si
e k ) = I(X k ),
I(X 06k 6p−1 e p ) 6= I(X p ).
et I(X
On peut alors envisager une autre manière d’approcher l’intégrale d’une fonction ϕ en utilisant
des polynômes : remplaçons ϕ par son polynôme interpolateur de Lagrange Pϕ associé aux N
points deux à deux distincts ξ1 , . . . , ξN ∈ [a, b], c’est-à-dire
Z β Z β
ϕ(ξ) dξ ≈ Pϕ (ξ) dξ = I(ϕ). e
α α
On rappelle que Pϕ s’écrit dans la base de Lagrange composée des polynômes

Y X − ξj
Li = , 16i6N
ξi − ξj
j6=i
sous la forme
N
X
Pϕ = ϕ(ξi )Li .
i=1
De sorte que
Z N
βX N
X Z β N
X
I(ϕ)
e = ϕ(ξi )Li (ξ) dξ = ϕ(ξi ) Li (ξ) dξ = ωi ϕ(ξi ),
α i=1 i=1 α i=1
où l’on a posé

Z β
ωi = Li (ξ) dξ 1 6 i 6 N.
α
Cette méthode est donc un cas particulier de formule de quadrature, associée aux points ξi et
aux poids ωi , qu’on qualifie de type interpolation à N points.
Définition 3.5 – Formule de quadrature de type interpolation

Etant donnés N points distincts de [α, β] notés ξ1 , . . . , ξN , la formule de quadrature de
type interpolation associée aux points (ξi )16i6N est définie par
N
X Z β Y ξ − ξj
I(ϕ)
e = ωi ϕ(ξi ), avec ωi = dξ 1 6 i 6 N.
α j6=i ξi − ξj
i=1
Proposition 3.6
Une formule de quadrature à N points est d’ordre au moins N si, et seulement si, elle est
de type interpolation à N points.
⁄
Démonstration. Supposons que Ie est une formule de quadrature de type interpolation aux points
ξ1 , . . . , ξN
Z β
I(ϕ)
e = Pϕ (ξ) dξ,
α
où Pϕ est le polynôme d’interpolation de Lagrange associé à ϕ aux points ξ1 , . . . , ξN . Si ϕ est un

polynôme de degré inférieur on égal à N − 1, ϕ = Pϕ et donc I(ϕ)
e = I(ϕ), ce qui conclut.
Réciproquement soit Ie est une formule de quadrature à N points (distincts)
N
X
I(ϕ)
e = ωi ϕ(ξi )
i=1
qui est exacte sur RN −1 [X]. Considérons la base de Lagrange formée des Li ∈ RN −1 [X], 1 6 i 6 N ,
e j ) = β Lj (ξ) dξ de sorte que
R
associée aux (ξi )16i6N . On sait que ωj = I(L α
N
X N Z
X β Z N
βX Z β
I(ϕ)
e = ωi ϕ(ξi ) = Li (ξ)ϕ(ξi ) dξ = ϕ(ξi )Li (ξ) dξ = Pϕ (ξ) dξ,
i=1 i=1 α α i=1 α
et la formule est de type interpolation à N points. !
Proposition 3.7 – Erreur des méthodes de type interpolation

Si Ie est de type interpolation à N point sur [α, β] et ϕ est de classe C N alors on a
l’estimation
(β − α)N +1 (N )
I(ϕ) − I(ϕ)
e 6 kϕ k∞ .
(N + 1)!
Démonstration. On utilise l’estimation vue au Corollaire 2.8

kϕ(N ) k∞
ϕ(ξ) − P (ξ) 6 |ωN (ξ)|
N!
où ωn (ξ) = (ξ − ξ1 ) . . . (ξ − ξN ). En utilisant l’inégalité |ωN (ξ)| 6 (ξ − α)N , pour α 6 ξ 6 β, on obtient
β β
kϕ(N ) k∞ (β − α)N +1 (N )
Z Z
I(ϕ) − I(ϕ)
e = ϕ(ξ) − P (ξ) dξ 6 |ωN (ξ)| dξ 6 kϕ k∞ ,
α α N! (N + 1)!
Ainsi, la proposition précédente pourrait nous faire penser que, pour mieux approcher l’intégrale
d’une fonction, il est suffisant d’augmenter le nombre de points N dans la formule de quadrature.
Il n’en est rien comme nous l’avons vu au chapitre 2 consacré à l’interpolation : l’estimation
dépendant de la dérivée N ième ne garantie pas la convergence. L’étude du phénomène de Runge
a même été plus loin dans l’analyse : la différence entre la fonction et le polynôme interpolateur
peut tendre vers l’infini pour certaines fonctions. Nous devons donc changer de stratégie. La
figure 3.1 illustre la divergence du calcul de l’intégrale approchée en utilisant directement le
polynôme interpolateur de Lagrange (avec des points équirépartis ou les points de Tchebychev).
La figure a été construite en calculant l’intégrale du polynôme interpolateur pour la fonction
x 7→ 1/(1 + x2 ) sur [0, 5].
⁄
Figure 3.1 – Illustration de la non convergence de la méthode d’intégration par interpolation.
1.2 Formules de quadratures élémentaires et composées
Soit f : [a, b] → R une fonction continue. Essayons d’approcher

Z b
I(f ) = f (x) dx.
a
Reprenons notre exemple précédent de somme de Riemann In (f ) = S(f, Σn ) où Σn est la

subdivision régulière donnée par h = (b − a)/n, xj = a + jh, 0 6 j 6 n et ξj = xj pour
0 6 j 6 N − 1. On a une somme de n termes de la forme (xj+1 − xj )f (xj ), or par la relation de
Chasles, I(f ) s’écrit aussi comme une somme de n termes
n−1
X Z xj+1
I(f ) = f (x) dx,
j=0 xj
Rx
de sorte qu’approcher I(f ) par In (f ) consiste à remplacer xjj+1 f (x) dx par (xj+1 − xj )f (xj ).
Rx
Or lorsque f est positive, xjj+1 f (x) dx représente l’aire sous le graphe de f entre les abscisses
xj et xj+1 , et (xj+1 − xj )f (xj ) représente l’aire du rectangle de même base et de hauteur f (xj ),
c’est-à-dire la valeur de la fonction à gauche de l’intervalle. On a donc découpé l’intervalle [a, b]
en petits intervalles sur lesquels on a approché l’aire sous le graphe de f par une formule de
quadrature élémentaire, appelée ici méthode des rectangles à gauche
Z xj+1
f (x) dx ≈ (xj+1 − xj )f (xj ),
xj
puis on a sommé pour obtenir une formule de quadrature composée

Z b N
X −1
f (x) dx ≈ (xj+1 − xj )f (xj ).
a j=0
⁄
Nous allons à présent définir des méthodes de quadrature élémentaires (qui seront de type
interpolation vues à la section précédente) plus précises que la méthode des rectangles à gauche
afin de rester dans l’esprit des sommes de Riemann mais en améliorant la vitesse de convergence.
Remarquons tout d’abord qu’il est suffisant de définir les formules de quadratures élémentaires
sur l’intervalle [−1, 1]. En effet, le changement de variables affine
φj : [−1, 1] −→ [xj , xj+1 ]

(3.1)
t 7−→ 1−t 1+t
2 xj + 2 xj+1
permet de ramener le calcul de l’intégrale sur [xj , xj+1 ] à un calcul sur [−1, 1] selon la formule
Z xj+1
xj+1 − xj 1
Z
f (x) dx = f (φj (t)) dt,
xj 2 −1
de sorte qu’il suffira d’appliquer les formules de quadrature élémentaires sur l’intervalle [−1, 1]
aux fonctions ϕj = f ◦ φj .
Définition 3.8 – Formule de quadrature élémentaire – composite

On appelle formule de quadrature élémentaire une formule de quadrature Ie sur [−1, 1]
N
X
Ie (ϕ) = ωi ϕ(ξi ) pour ϕ ∈ C 0 ([−1, 1]).
i=1
Etant donnée une subdivision de l’intervalle [a, b] en n sous-intervalles de la forme
a = x0 < x1 < . . . < xn = b,
la formule de quadrature élémentaire Ie induit une formule de quadrature composite Ic

sur [a, b] qui s’écrit
n−1 n−1 N
X xj+1 − xj X xj+1 − xj X
Ic (f ) = Ie (f ◦ φj ) = ωi f (xi,j ), (3.2)
2 2
j=0 j=0 i=1
avec
1 − ξi 1 + ξi
xi,j = φj (ξi ) = xj + xj+1 , 1 6 i 6 N, 0 6 j 6 n − 1.
2 2
2 Méthodes de quadrature classiques
Dans cette section, quelques-unes des méthodes de quadratique parmi les plus classiques seront
présentées, puis analysées afin de déterminer une estimation de l’erreur commise. Commençons
par appliquer la méthode générale des méthodes composites dans les cas les plus simples.
Précisons les notations. Nous noterons
Z 1 N
X
I0 (ϕ) = ϕ(ξ) dξ, Ie (ϕ) = ωi ϕ(ξi ),
−1 i=1
n−1 n−1
(3.3)
b
xj+1 − xj xj+1 − xj
Z X X
I(f ) = f (x) dx = I0 (f ◦ φj ), Ic (f ) = Ie (f ◦ φj ),
a 2 2
j=0 j=0
⁄
Nom points d’interpolation formule composite

n−1
b−a X
rectangle à droite 1 f (xj+1 )
n j=0
n−1
b−a X
rectangle à gauche −1 f (xj )
n j=0
n−1
b − a X xj+1 + xj
point milieu 0 f
n j=0 2
n−1
b−a X
trapèze −1, 1 f (xj + 1) + f (xj )
2n j=0
n−1
b−a X x
j+1 + xj

Simpson −1, 0, 1 f (xj + 1) + 4f + f (xj )
6n j=0 2
Table 3.1 – Les méthodes de quadratures classiques : méthodes composées à partir des méthodes
de type interpolation.
1−t 1+t
où φj : [−1, 1] → [xj , xj+1 ] t 7→ 2 xj + 2 xj+1 est définie en (3.1).
2.1 Une liste des méthodes classiques
Nous considérons ici différentes méthodes de quadrature élémentaires de type interpolation et

nous précisons la formule de quadrature composite obtenue. Nous choisissons pour simplifier les
formules une subdivision régulière de l’intervalle [a, b]. C’est-à-dire que nous avons xi = a + ih
avec h = (b − a)/n. Le tableau 3.1 récapitule les formules classiques de quadrature.
Détaillons à présent comment sont obtenues ces formules. Pour la méthode des rectangles à droite,
nous avons un seul point d’interpolation : le point ξ1 = 1. Ainsi, la fonction ϕ est remplacée par
son polynôme interpolateur de Lagrange au point ξ1 = 1. Nous pouvons écrire
Z 1
Pϕ (ξ) = ϕ(1) =⇒ Ie (ϕ) = ϕ(1) dξ = 2ϕ(1).
−1
Nous en déduisons à l’aide de la formule (3.2)
n−1
b−aX
Ic (f ) = f (xj+1 ). (Rectangles à droite)
n
j=0
La formule des rectangles à gauche s’obtient exactement de la même manière en choisissant

comme polynôme interpolateur le polynôme interpolateur de Lagrange au point ξ1 = −1 :
Z 1
Pϕ (ξ) = ϕ(−1) =⇒ Ie (ϕ) = ϕ(−1) dξ = 2ϕ(−1).
−1
Nous en déduisons à l’aide de la formule (3.2)
n−1
b−aX
Ic (f ) = f (xj ). (Rectangles à gauche)
n
j=0
⁄
La formule du point milieu est obtenue quant à elle en choisissant le polynôme interpolateur de
Lagrange au point ξ1 = 0 :
Z 1
Pϕ (ξ) = ϕ(0) =⇒ Ie (ϕ) = ϕ(0) dξ = 2ϕ(0).
−1
Dans ce cas, nous avons x1,j = (xj + xj+1 )/2 et
n−1
b − a X xj+1 + xj
Ic (f ) = f . (Points milieux)
n 2
j=0
La formule des trapèzes est obtenue en choisissant le polynôme interpolateur de Lagrange aux
deux points ξ1 = −1 et ξ2 = 1. Dans ce cas, le polynôme est de degré au plus 1, ainsi
Z 1
1−ξ 1+ξ
Pϕ (ξ) = ϕ(−1) + ϕ(1) =⇒ Ie (ϕ) = Pϕ (ξ) dξ = ϕ(−1) + ϕ(1).
2 2 −1
Dans ce cas, nous avons x1,j = xj et x2,j = xj+1 et
n−1
b − a Xh i
Ic (f ) = f (xj ) + f (xj+1 ) . (Trapèzes)
2n
j=0
Finalement, la méthode de Simpson est obtenue en choisissant le polynôme interpolateur de

Lagrange aux trois points ξ1 = −1, ξ2 = 0 et ξ3 = 1. Ce polynôme de degré au plus 2 s’écrit
ξ(ξ − 1) ξ(1 + ξ)
Pϕ (ξ) = ϕ(−1) + (1 + ξ)(1 − ξ)ϕ(0) + ϕ(1).
2 2
Nous avons donc (en utilisant que les contributions des parties impaires sont nulles)
Z 1 3 1
ξ(ξ − 1) ξ 1
ω1 = dξ = = ,
−1 2 6 −1 3
Z 1 1
ξ3

2 4
ω2 = (1 + ξ)(1 − ξ) dξ = 1 − =2− = ,
−1 3 −1 3 3
Z 1 3 1
ξ(ξ + 1) ξ 1
ω3 = dξ = = .
−1 2 6 −1 3
Nous en déduisons alors la formule composée
n−1
b − a Xh x
j+1 + xj
i
Ic (f ) = f (xj + 1) + 4f + f (xj ) . (Simpson)
6n 2
j=0
La figure 3.2 illustre la convergence de ces 5 méthodes d’intégration numérique. Les courbes ont
été obtenues en calculant l’erreur du calcul approché avec N points (N variant entre 1 et 100)
pour calculer Z 5
1
I= 2
dx.
0 1+x
En échelle logarithmique, l’erreur est asymptotiquement proche d’une droite de pente entière (1
pour les méthodes des rectangles, 2 pour les méthodes des points milieux et des trapèzes et 4
pour la méthode de Simpson). Les sections suivantes vont permettre de comprendre pourquoi
ces pentes sont observées.
⁄
Figure 3.2 – Illustration de la convergence des méthodes d’intégration numérique.
2.2 Analyse des méthodes des rectangles
Considérons la formule des rectangles à gauche
Ie (ϕ) = 2ϕ(−1)
pour ϕ ∈ C 0 ([−1, 1]). On cherche à calculer l’ordre de la méthode et quantifier l’erreur commise
par rapport à la valeur exacte de l’intégrale. Une illustration des méthodes des rectangles à
gauche et à droite se trouve à la figure 3.3.
Figure 3.3 – Illustration de la méthode du rectangle.
Proposition 3.9
La méthode des rectangles à gauche est d’ordre 1.
⁄
Démonstration. On sait déjà que c’est une méthode d’interpolation à 1 point donc elle est d’ordre au
moins 1. Pour montrer qu’elle n’est pas d’ordre plus grand il suffit d’exhiber un polynôme
R1 P de degré
1 tel que Ie (P ) 6= I0 (P ). Le choix P = X convient : Ie (P ) = −2 tandis que I0 (P ) = −1 x dx = 0. !
Proposition 3.10 – Erreur de la méthode des rectangles à gauche

Si ϕ est de classe C 1 sur [−1, 1], il existe c ∈] − 1, 1[ tel que l’erreur de la méthode de
quadrature élémentaire s’écrit
Ee (ϕ) = I0 (ϕ) − Ie (ϕ) = 2ϕ0 (c).
Si f est de classe C 1 sur [a, b] l’erreur de quadrature de la méthode composite associée à

une subdivision régulière de pas h est majorée par
b−a
|Ec (f )| = |I(f ) − Ic (f )| 6 hkf 0 k∞ .
2
Démonstration. Soit φ une primitive de ϕ sur l’intervalle [−1, 1]. Effectuons un développement de Taylor
avec reste de Lagrange de φ au point −1 :
∃c ∈] − 1, 1[: φ(1) = φ(−1) + 2φ0 (−1) + 2φ00 (c).
or I0 (ϕ) = φ(1) − φ(−1) et φ0 = ϕ, φ00 = ϕ0 donc

Z 1
I0 (ϕ) − Ie (ϕ) = ϕ(t) dt − 2ϕ(−1) = 2ϕ0 (c).
−1
Pour f ∈ C 1 ([a, b], R), nous avons grâce aux formules (3.3)
n−1
X xj+1 − xj
I(f ) − Ic (f ) = I0 (f ◦ φj ) − Ie (f ◦ φj )
j=0
2
n−1
X n−1
X
= (xj+1 − xj )(f ◦ φj )0 (cj ) = (xj+1 − xj )f 0 (φj (cj ))φ0j (cj ),
j=0 j=0
où cj ∈] − 1, 1[, 0 6 j 6 n−1. Comme φ0j (cj ) = (xj+1 − xj )/2, nous obtenons
n−1
X (xj+1 − xj )2 0 b−a
I(f ) − Ic (f ) 6 kf k∞ = hkf 0 k∞ ,
j=0
2 2
La formule des rectangles à droite se traite de la même manière.

Remarquons à ce niveau que nous pouvons faire le lien entre l’ordre de la méthode de quadrature
élémentaire (la propriété d’être une formule exacte pour les polynôme d’un certain degré) et
l’ordre de convergence de la méthode de quadrature. Nous pouvons en effet démontrer une
estimation de l’erreur en O(h) sans utiliser le fait que la méthode est de type interpolation. Il
est suffisant de remarquer que la méthode des rectangles (à gauche ou à droite, cela ne change
rien) est une formule de quadrature d’ordre 1, c’est-à-dire que Ie (P ) = I0 (P ) dès que P ∈ R0 [X].
Or nous pouvons écrire
∀ξ ∈ [−1, 1] ∃η ∈ [−1, 1] : ϕ(ξ) = ϕ(0) + ϕ0 (η)ξ.
En notant P le polynôme obtenu en tronquant le développement de Taylor à l’ordre le plus bas
(c’est-à-dire P = ϕ(0)), nous avons
∀ξ ∈ [−1, 1] ϕ(ξ) − P (ξ) 6 kϕ0 k∞ |ξ|.
⁄
Nous en déduisons par linéarité de I0 et de Ie que
I0 (ϕ) − Ie (ϕ) = I0 (ϕ) − I0 (P ) + Ie (P ) − Ie (ϕ) ,

Z 1 Z 1
I0 (ϕ) − I0 (P ) 6 ϕ(ξ) − P (ξ) dξ 6 kϕ0 k∞ |ξ| dξ = kϕ0 k∞ ,
−1 −1
Ie (P ) − Ie (ϕ) = 2 P (−1) − ϕ(−1) 6 2kϕ0 k∞ .
En reprenant la preuve de la proposition précédente, c’est cette inégalité qui permet ensuite
d’obtenir la majoration |I(f ) − Ic (f )| 6 O(h). Ainsi, si la formule de quadrature est d’ordre 1,
l’erreur de la méthode composée associée est en O(h).
2.3 Analyse de la méthode des trapèzes
Considérons la formule des trapèzes
Ie (ϕ) = ϕ(−1) + ϕ(1)
pour ϕ ∈ C 0 ([−1, 1]). Une illustration de la méthode des trapèzes se trouve à la figure 3.4. On
cherche à calculer l’ordre de la méthode et à quantifier l’erreur commise par rapport à la valeur
exacte de l’intégrale.
Figure 3.4 – Illustration de la méthode du point milieu à gauche et du trapèze à droite.
Proposition 3.11
La méthode des trapèzes est d’ordre 2.
Démonstration. On sait déjà que c’est une méthode d’interpolation à 2 points donc elle est d’ordre au
moins 2. Pour montrer qu’elle R 1n’est pas d’ordre plus grand, nous choisissons le polynôme P = X 2 . Nous
avons Ie (P ) = 2 et I0 (P ) = −1 x2 dx = 2/3. !
⁄
Proposition 3.12 – Erreur de la méthode des trapèzes

Si ϕ est de classe C 2 sur [−1, 1], l’erreur de la méthode de quadrature élémentaire est
majorée par
2
|Ee (ϕ)| = |I0 (ϕ) − Ie (ϕ)| 6 kϕ00 k∞ .
3
b−a
|Ec (f )| = |I(f ) − Ic (f )| 6 h2 kf 00 k∞ .
12
Démonstration. Comme la méthode des trapèzes est une méthode de type interpolation à deux points,
nous appliquons la proposition 3.7 : pour ϕ ∈ C 2 ([−1, 1]) nous avons
23 00 4
|Ee (ϕ)| 6 kϕ k∞ = kϕ00 k∞ .
3! 3
Pour améliorer la constante, il est possible de reprendre la preuve de la proposition (3.7) pour écrire
1
x3 i1
Z
1 00 1 00 h 2
|Ee (ϕ)| 6 kϕ k∞ (1 − x)(1 + x) dx 6 kϕ k∞ x − = kϕ00 k∞ .
2! −1 2 3 −1 3
Pour f ∈ C 2 ([a, b], R), nous avons grâce aux formules (3.3)
n−1
X xj+1 − xj
|I(f ) − Ic (f )| 6 I0 (f ◦ φj ) − Ie (f ◦ φj )
j=0
2
n−1
1X
6 (xj+1 − xj )k(f ◦ φj )00 k∞ .
3 j=0
Comme φ0j (x) = (xj+1 − xj )/2, nous avons
h 0 h2 00
(f ◦ φj )0 (ξ) = f ◦ φj (ξ), (f ◦ φj )00 (ξ) = f ◦ φj (ξ).
2 4
Nous obtenons donc
n−1
1 X (xj+1 − xj )3 0 b−a
I(f ) − Ic (f ) 6 kf k∞ = h2 kf 0 k∞ ,
3 j=0 4 12
2.4 Analyse de la méthode du point milieu
Considérons la formule du point milieu
Ie (ϕ) = 2ϕ(0)
pour ϕ ∈ C 0 ([−1, 1]). Une illustration de la méthode des points milieux se trouve à la figure 3.4.
Proposition 3.13
La méthode du point milieu est d’ordre 2.
⁄
Démonstration. Comme c’est une méthode d’interpolation à 1 point, la méthode est au moins d’ordre
1. Testons la méthode sur la base canonique des polynômes :
Z 1
P = 1, I0 (P ) = 1 dx = 2, Ie (P ) = 2P (0) = 2,
−1
Z 1
P = X, I0 (P ) = x dx = 0, Ie (P ) = 2P (0) = 0,
−1
Z 1
2
P = X 2, I0 (P ) = x2 dx = , Ie (P ) = 2P (0) = 0.
−1 3
La méthode est donc d’ordre exactement 2. !
Proposition 3.14 – Erreur de la méthode du point milieu

Si ϕ est de classe C 2 sur [−1, 1], il existe c ∈] − 1, 1[ tel que l’erreur de la méthode de
quadrature élémentaire s’écrit
1
Ee (ϕ) = I0 (ϕ) − Ie (ϕ) = ϕ00 (c).
3
b−a
|Ec (f )| = |I(f ) − Ic (f )| 6 h2 kf 00 k∞ .
24
Rx
Démonstration. On pose φ(x) = −1
ϕ(t) dt et on fait un développement de Taylor-Lagrange à l’ordre
3 en 0
1 1
φ(−1) = φ(0) − φ0 (0) + φ00 (0) − φ(3) (c1 ),
2 6
0 1 00 1
φ(1) = φ(0) + φ (0) + φ (0) + φ(3) (c2 ),
2 6
R1
pour certains c1 , c2 ∈] − 1, 1[. Or −1 ϕ(t) dt = φ(1) − φ(−1) et φ0 = ϕ donc en faisant la différence :
Z 1
1
ϕ(t) dt = 2ϕ(0) + (ϕ00 (c1 ) + ϕ00 (c2 ))).
−1 6
La fonction ϕ00 étant continue, d’après le théorème des valeurs intermédiaires, il existe c ∈]a, b[ tel que
(ϕ00 (c1 ) + ϕ00 (c2 ))/2 = ϕ00 (c), d’où
Z 1
1
ϕ(t) dt − 2ϕ(0) = ϕ00 (c).
−1 3
En utilisant le changement de variable φj sur l’intervalle [xj , xj+1 ], on en déduit que si f est de classe
C 2 sur [a, b],
n−1
X xj+1 − xj
I(f ) − Ic (f ) = I0 (f ◦ φj ) − Ie (f ◦ φj )
j=0
2
n−1 n−1
X (xj+1 − xj )3
X xj+1 − xj
= (f ◦ φj )00 (cj ) = f (2) (φj (cj )).
j=0
6 j=0
24
Donc
b−a
|I(f ) − Ic (f )| 6 h2 kf 00 k∞
24
⁄
2.5 Analyse de la méthode de Simpson
Considérons la formule de Simpson
1 4 1
Ie (ϕ) = ϕ(−1) + ϕ(0) + ϕ(1)
3 3 3
pour ϕ ∈ C 0 ([−1, 1]). Une illustration de la méthode de Simpson se trouve à la figure 3.5.
Figure 3.5 – Illustration de la méthode de Simpson.
Proposition 3.15
La méthode de Simpson est d’ordre 4.
Démonstration. Comme c’est une méthode d’interpolation à 3 point, la méthode est au moins d’ordre
3. Testons la méthode sur la base canonique des polynômes :
Z 1
1 4 1
P = 1, I0 (P ) = 1 dx = 2, Ie (P ) = P (−1) + P (0) + P (1) = 2,
−1 3 3 3
Z 1
1 4 1
P = X, I0 (P ) = x dx = 0, Ie (P ) = P (−1) + P (0) + P (1) = 0,
−1 3 3 3
Z 1
2 1 4 1 2
P = X 2, I0 (P ) = x2 dx = , Ie (P ) = P (−1) + P (0) + P (1) = ,
−1 3 3 3 3 3
Z 1
1 4 1
P = X 3, I0 (P ) = x3 dx = 0, Ie (P ) = P (−1) + P (0) + P (1) = 0,
−1 3 3 3
Z 1
2 1 4 1 2
P = X 4, I0 (P ) = x4 dx = , Ie (P ) = P (−1) + P (0) + P (1) = .
−1 5 3 3 3 3
La méthode est donc d’ordre exactement 4. !
⁄
Proposition 3.16 – Erreur de la méthode de Simpson

Si ϕ est de classe C 4 sur [−1, 1], l’erreur de la méthode de quadrature élémentaire est
majorée par
2
|Ee (ϕ)| = |I0 (ϕ) − Ie (ϕ)| 6 kϕ(4) k∞ .
45
b−a
|Ec (f )| = |I(f ) − Ic (f )| 6 h4 kf (4) k∞ .
720
Rx
Démonstration. On pose φ(x) = −1
ϕ(t) dt et on fait un développement de Taylor-Lagrange à l’ordre
5 en 0
1 1 (3) 1 (4) 1 (5)
φ(−1) = φ(0) − φ0 (0) + φ00 (0) − φ (0) + φ (0) − φ (c1 ),
2 6 24 120
1 1 (3) 1 (4) 1 (5)
φ(1) = φ(0) + φ (0) + φ00 (0) +
0
φ (0) + φ (0) + φ (c2 ),
2 6 24 120
pour certains c1 , c2 ∈] − 1, 1[. Ainsi
1 1
I0 (ϕ) = φ(1) − φ(−1) = 2ϕ(0) + ϕ00 (0) + ϕ(4) (c1 ) + ϕ(4) (c2 ) .

3 120
Nous effectuons alors un développement de Taylor-Lagrange à l’ordre 4 de la fonction ϕ en 0 également
1 1 (3) 1 (4)
ϕ(−1) = ϕ(0) − ϕ0 (0) + ϕ00 (0) − ϕ (0) + ϕ (c3 ),
2 6 24
1 1 (3) 1 (4)
ϕ(1) = ϕ(0) + ϕ (0) + ϕ00 (0) +
0
ϕ (0) + ϕ (c4 ),
2 6 24
pour certains c3 , c4 ∈] − 1, 1[. Ainsi
1 1 (4)
Ie (ϕ) = 2ϕ(0) + ϕ00 (0) + ϕ (c3 ) + ϕ(4) (c4 ) .

3 72
2 2 (4) 2
|Ee (ϕ)| = |I0 (ϕ) − Ie (ϕ)| 6 + kϕ k∞ = kϕ(4) k∞ .
120 72 45
En utilisant le changement de variable φj sur l’intervalle [xj , xj+1 ], on en déduit que si f est de classe
C 4 sur [a, b],
n−1
X xj+1 − xj
I(f ) − Ic (f ) = I0 (f ◦ φj ) − Ie (f ◦ φj )
j=0
2
n−1 n−1
X (xj+1 − xj )5
X xj+1 − xj
= (f ◦ φj )(4) (cj ) = f (2) (φj (cj )).
j=0
45 j=0
720
Donc
b−a
|I(f ) − Ic (f )| 6 h4 kf (4) k∞
720
⁄
2.6 Contrôle de l’erreur
En utilisant les théorèmes de majoration de l’erreur de quadrature, il est possible de construire

un maillage de l’intervalle [a, b] afin de garantir que l’erreur soit inférieure à un certain seuil .
Supposons que nous aillons à disposition une formule de quadrature d’ordre k, c’est-à-dire qu’elle
vérifie (avec les notations utilisées dans cette section)
|Ee (ϕ)| 6 Ckϕ(k) k∞
pour une certaine constante C. Si l’on considère une subdivision de l’intervalle [a, b] non néces-
sairement régulière x0 = a < x1 < . . . < xn−1 < xn = b, l’erreur globale de la méthode est alors
majorée par
n−1
X (xj+1 − xj )k
|I(f ) − Ic (f )| 6 C(b − a) sup |f (k) (t)|.
2k t∈[xj ,xj+1 ]
j=0
Il est ainsi possible de choisir une subdivision telle que chacun des termes de la somme soit
inférieur à un certain /n, ce qui implique que l’erreur globale sera inférieure à ce .
2 2
La figure 3.6 permet de visualiser pour la fonction x 7→ e−x − e−25x que l’on intègre sur l’in-
tervalle [−4, 4] des maillages construits pour que chaque méthode produise une erreur inférieure
à 10−3 .
Figure 3.6 – Illustration de maillages adaptés pour les différentes méthodes : à gauche = 10−3 ,
à droite = 10−4 .
Dans le tableau 3.2, nous avons affiché l’erreur de quadrature qui est effectivement inférieure à
10−3 et le nombre de points du maillage. On observe que les méthodes d’ordre 1 (les méthodes
des rectangles) nécessitent beaucoup plus de points que les méthodes d’ordre plus élevé : plus
de 17000 alors que la méthode de Simpson d’ordre 4 seulement 41. Lorsque le seuil est diminué
(dans le tableau 3.3, = 10−4 ), la conclusion est encore plus flagrante. Grossièrement, pour
diviser le seuil par 10 avec une méthode d’ordre 1, il faut environ 10 fois plus de points. Pour
une méthode d’ordre plus élevé, la croissance est beaucoup plus lente.
La conclusion est que (lorsque la fonction que l’on intègre est régulière), il est préférable d’utiliser
des méthodes d’ordre élevé. Même si chaque calcul coûte plus cher (l’évaluation de la méthode
de Simpson nécessite environ 2 fois plus de calcul que les autres méthodes), comme elle nécessite
un maillage beaucoup plus grossier, elle est avantageuse.
⁄
Nom de la méthode erreur nombre de points dans le maillage

rectangle à droite 3.323E − 05 17693
rectangle à gauche 8.314E − 06 17693
point milieu 2.198E − 04 250
trapèze 2.245E − 04 189
Simpson 3.233E − 05 41
Table 3.2 – Illustration de l’adaptation du maillage pour un seuil = 10−3
Nom de la méthode erreur nombre de points dans le maillage

rectangle à droite 1.982E − 06 203638
rectangle à gauche 1.268E − 06 203638
point milieu 2.759E − 05 831
trapèze 2.895E − 05 569
Simpson 2.407E − 06 69
Table 3.3 – Illustration de l’adaptation du maillage pour un seuil = 10−4
3 Les algorithmes des méthodes de quadrature classiques
Dans cette section, nous proposons une version python des algorithmes qui ont été étudiés dans
ce chapitre. Ils n’ont pas la prétention d’être parfait... et sont limités à des maillages uniformes,
c’est-à-dire que les points sont définis par xi = a + ih avec h = (b − a)/n.
Toutes les fonctions proposées prennent quatre arguments : la fonction f à intégrer, deux réels a
et b qui définissent le domaine d’intégration et un entier N qui détermine le nombre de découpage
pour la méthode composée.
Algorithme 3.1 – Rectangles à gauche
def rectangles_gauche (f , a , b , N ):
x , h = np . linspace (a , b , N +1 , retstep = True )
xl = x [: -1]
return h * np . sum ( f ( xl ))
Algorithme 3.2 – Rectangles à droite
def rectangles_droite (f , a , b , N ):
xr = x [1:]
return h * np . sum ( f ( xr ))
Algorithme 3.3 – Trapèzes
def trapezes (f , a , b , N ):
xl , xr = x [: -1] , x [1:]
return .5* h *( np . sum ( f ( xl )) + np . sum ( f ( xr )))
⁄
Algorithme 3.4 – Points milieux

def points_milieux (f , a , b , N ):
xm = .5*( x [: -1]+ x [1:])
return h * np . sum ( f ( xm ))
Algorithme 3.5 – Simpson
def simpson (f , a , b , N ):
xl , xr = x [: -1] , x [1:]
xm = .5*( xl + xr )
return h /6*(
np . sum ( f ( xl )) + 4* np . sum ( f ( xm )) + np . sum ( f ( xr ))
)
⁄
4 Résolution d’équations ordinaires
Ce chapitre est consacré à la recherche de solution approchée d’équations non linéaires scalaires.
Soit I un intervalle de R et soit f : I → R continue. On considère le problème suivant
Chercher x? ∈ I tel que f (x? ) = 0. (4.1)
Définition 4.1
Toute solution x? ∈ I de (4.1) est dite racine ou zéro de f dans I.
Si f est de classe C r sur I avec r > 1 et x? une racine de f dans I, alors
0
. x? est dite racine simple si f (x? ) 6= 0,
. x? est dite racine de multiplicité p < r, si f (k) (x? ) = 0 pour 0 6 k < p et
f (p) (x? ) 6= 0.
1 Exemples d’application
Nous présentons deux exemples classiques pour lesquels il est nécessaire de résoudre de telles
équations.
1.1 Schémas numériques pour équations différentielles ordinaires
L’une des utilisations importantes est la résolution numérique des équations différentielles or-
dinaires. La construction des schémas numériques conduit parfois à des méthodes implicites,
c’est-à-dire nécessitant la résolution d’équations ordinaires. Ces méthodes peuvent présenter
des avantages numériques même si elles peuvent parfois sembler plus compliquées. Pour s’en
convaincre, considérons le problème simple
u0 (t) = −u(t) t ∈]0, T [, avec u(0) = u0 ,
dont la solution est donnée par u(t) = u0 e−t .

En considérant un pas de temps constant ∆t > 0, le schémas d’Euler explicite appliqué à la
résolution de cette équation conduit à la construction de la suite ue = (uen )n∈N définie par
uen+1 = (1 − ∆t)uen , n > 0, (4.2)
alors que le schéma d’Euler implicite conduirait à la suite ui = (uin )n∈N
uin+1 = uin − ∆tuin+1 , n > 0.
Pour uin donné, la détermination de uin+1 se ramène ici à la résolution de l’équation x = uin −x∆t
dont la solution est évidente car l’équation à résoudre est linéaires... Nous avons donc
1
uin+1 = ui , n > 0. (4.3)
1 + ∆t n
⁄
Figure 4.1 – Comparaison des méthodes explicites et implicites
La figure 4.1 montre les premiers termes de la suite obtenues par les méthodes explicites et
implicites avec différentes valeurs du pas de temps ∆t.
Pour la méthode explicite, on montre que, si l’on prend u0 > 0,
. si 0 < ∆t 6 1, la suite est décroissante et converge vers 0, ce qui est le bon comportement ;
. si 1 < ∆t < 2, la suite prend alternativement des valeurs positives et négatives mais
converge vers 0 ;
. si ∆ > 2, la suite diverge en prenant alternativement des valeurs positives et négatives de
plus en plus grandes en valeur absolue.
Pour la méthode implicite, on montre que, si l’on prend u0 > 0, la suite reste positive et tend vers
0 quelque soit la valeur du pas de temps ∆t. Ainsi, le comportement de la solution approchée reste
toujours cohérent avec celui de la solution exacte. La méthode implicite semble plus intéressante.
Cependant il n’est pas toujours possible de résoudre directement le problème posé par les schémas
implicites. En particulier, lorsque la fonction n’est pas linéaire, la résolution n’admet pas toujours
une seule solution et les solutions n’ont pas toujours d’expression analytique. Voici un exemple
⁄
Chapitre 4. Résolution d’équations ordinaires
illustrant ce propos.
u0 (t) = eu(t) t ∈]0, T [, avec u(0) = u0 . (4.4)
Ici le schéma d’Euler implicite s’écrit
un+1 = un + ∆teun+1 , n > 0.
Pour n > 0, cette équation ne pourra être résoluble que de manière approchée (numériquement).
Sachant que cette résolution devra être faite à chaque itération, il faudra qu’elle soit efficace,
c’est-à-dire rapide (peu couteuse en temps de calcul et taille mémoire) et précise (afin de ne pas
accroı̂tre l’erreur locale de la méthode à un pas implicite).
1.2 Méthode de tir pour les problèmes aux limites du second ordre
Dans un problème de tir à canon, on est souvent amené à déterminer l’angle vertical d’orientation
du canon afin d’atteindre une cible précise.
En une dimension d’espace, le problème peut se mettre sous la forme (4.5) où f est une certaine
fonction donnée de manière implicite sous la forme suivante.
Etant donnés deux réels x0 et xT donnés, on cherche une solution du problème
 00
x (t) = f (t, x(t)), t ∈]0, T [,

x(0) = x0 , (4.5)

x(T ) = xT .

On obtient ici un problème aux limites et non un problème de Cauchy pour une équation
différentielle ! Si l’on dispose d’un bon solveur d’EDOs, on peut en faire usage. On peut en effet
transformer ce problème en un problème de Cauchy (en une edo), en introduisant une inconnue
nécessaire pour définir une condition initiale.
En effet, si pour v ∈ R on considère le problème
 00
z (t) = f (t, z(t)), t ∈]0, T [,

z(0) = x0 , (4.6)
 0

z (0) = v,
on a alors un problème de Cauchy pour une EDO. Si l’on pose alors g(t, v) la solution de cette
EDO, en l’évaluant à t = T on définit une fonction g(T, v) de la seule variable v.
Le problème de détermination de la vitesse initiale revient alors à résoudre l’équation d’inconnue
v : G(v) = 0, avec G(v) = g(T, v) − xT .
Il apparaı̂t donc que dans certaines équations non linéaires, l’expression de la fonction dont
on cherche une racine peut être compliquée : ici une évaluation de la fonction G nécessite une
résolution d’équation différentielle. Il est donc nécessaire de chercher des méthodes numériques
de résolution d’équations non-linéaires qui soient efficaces et si possible sans recours au calcul
des dérivées.
2 Position correcte du problème
Avant de nous lancer dans la construction des schémas numériques pour le problème 4.1 d’équa-
tions non-linéaires, il est important de répondre aux questions suivantes.
⁄
1. Le problème admet-il une solution ?

2. Si oui cette solution est-elle unique ?
3. Cette solution dépend-t-elle continûment des données ? Autrement dit est-elle stable vis-
à-vis des données du problème ?
4. La solution du problème si elle existe a-t-elle une régularité particulière ?
La réponse à ces questions est ce que l’on appelle vérification du côté bien posé du problème.
Elle est nécessaire lorsqu’on souhaite résoudre numérique (de manière approchée) tout problème
donné. Si le point 3 est négligé, les conséquences peuvent être dramatiques au niveau numérique.
2.1 Existence de solution
Puisqu’on est dans un cadre scalaire (c’est-à-dire f : R 7→ R), l’existence de solution est une
simple conséquence du théorème de Rolle (Thm. 1.4).
Proposition 4.2 – cas où l’équation est sous la forme f (x) = 0

Soit I = [a, b] et f : I → R continue telle que f (a)f (b) < 0, alors il existe x? ∈ I tel que
f (x? ) = 0.
Proposition 4.3 – cas où l’équation est sous la forme x = g(x)

Soit I = [a, b] et g : I → R continue telle que g(I) ⊂ I, alors il existe x? ∈ I tel que
g(x? ) = x? .
Démonstration. On applique le résultat précédent à f (x) = g(x) − x !
Remarque 4.4 – accessibilité de l’intervalle I

Les résultats que nous venons de présenter, font l’hypothèse que l’intervalle I est connu.
Mais ce n’est pas toujours le cas en pratique. Ainsi, I est lui même inconnu et est appelé
domaine d’attraction de la solution dans le cas d’un problème de point fixe.
Lorsque l’équation est sous la forme f (x) = 0, un procédé constructif permet aussi de définir
l’intervalle I : on se fixe un point a et un pas ∆. Puis on se déplace selon ce pas de part et
d’autre de a, jusqu’au premier changement de signe de f . On définit alors l’intervalle I. Pour ∆
suffisamment petit, cette procédure approche directement la racine de l’équation et on appelle
la schéma ainsi construit méthode de recherche incrémentale.
On retient alors à ce stade que certaines preuves de la position correcte des équations non-
linéaires cachent en elles mêmes des schémas numériques, certes basiques, mais exploitables
pour la résolution de ces équations.
2.2 Notion de conditionnement
Intéressons nous à présent à la stabilité. C’est-à-dire à l’effet, sur la solution, des perturbations
sur les données.
⁄
m x̄? |x̄? − x? | conditionnement

1 1.000000010000000 1.000E-08 1.000E+00
2 0.999900000000000 1.000E-04 1.000E+04
3 1.002154434690032 2.154E-03 2.154E+05
4 0.990000000000000 1.000E-02 1.000E+06
5 1.025118864315096 2.512E-02 2.512E+06
Table 4.1 – Influence du conditionnement sur la recherche de zéro
Définition 4.5 – Conditionnement

On appelle nombre de conditionnement d’un algorithme, le facteur d’amplification (ou de
réduction) de l’erreur d’évaluation d’un algorithme.
Si ce facteur s’intéresse aux erreurs relatives on parle de conditionnement relatif ou sim-
plement conditionnement. Si par contre ce facteur ne s’intéresse qu’aux erreurs absolues,
on parle de conditionnement absolu.
Le problème sera dit bien conditionné si le conditionnement est petit, c’est-à-dire proche
de 1. Il sera dit mal conditionné si ce nombre est très grand par rapport à 1.
Proposition 4.6 – conditionnement absolu

Supposons f de classe C p , p > 0 et soient x? et x̄? solutions de
f (x? ) = 0, f (x̄? ) + η(x̄? ) = 0,
où η est une perturbation (erreur de mesure) de f supposée bornée au voisinage de x? .

Si x? est une racine de multiplicité m < p de f alors
1/m
η(x̄? )
|x̄? − x? | . m! , (4.7)
f (m) (x? )
Démonstration. Effectuons le développement limité de la fonction f autour du point x? avec un reste

de Lagrange. Nous avons
1 ?
f (x̄? ) = f (x? ) + (x̄ − x? )m f (m) (ξ),
m!
où ξ est dans l’intervalle [x̄? , x? ] ou [x? , x̄? ]. Nous en déduisons que
1 ?
−η(x̄? ) = (x̄ − x? )m f (m) (ξ).
m!
Par continuité, nous obtenons la majoration du théorème
1/m
η(x̄? )
|x̄? − x? | . m! ,
f (m) (x? )
ce qui termine la preuve. Notons que, par définition de la multiplicité, le terme f (m) (x? ) n’est pas
nul. !
−1/m
Ainsi, le conditionnement absolu est donné par κ(f, x∗ ) = m! f (m) (x∗ ) |η(x̃)|1/m−1 . Ce qui
montre que la résolution des équations non-linéaires à racines multiples sera un problème moins
bien conditionné (et même mal conditionné) que celui de la recherche de racines simples.
Prenons pour exemple la recherche du zéro des fonctions fm : x 7→ (x − 1)m en supposant que
l’erreur η est constante (η(x) = −10−8 la valeur de η est choisie négative pour que l’équation ait
toujours des solutions). Le tableau 4.1 récapitule les résultats obtenus pour différentes valeurs
de la multiplicité m.
⁄
Ainsi, si l’on détermine le zéro avec 8 chiffres significatifs (précision du calcul de 10−8 ), le
conditionnement pour m = 1 permet de conserver la même erreur sur la solution, tandis que
pour m = 4 par exemple, l’erreur sur la solution est de 10−2 , c’est-à-dire que l’erreur est amplifiée
d’un facteur 106 . Il faut avoir conscience de ce phénomène d’amplification des erreurs lorsque
l’on fait du calcul scientifique.
2.3 vitesse de convergence - ordre de convergence
Comme nous le verrons, la recherche de solution approchée conduira à la construction de suites

qui convergent vers la solution du problème. Il est donc nécessaire de quantifier cette convergence.
Pour cela un rappel sur l’ordre de convergence des suites est nécessaire.
Définition 4.7 – ordre de convergence

Soit (xn )n∈N une suite qui converge vers une limite x? . La convergence de la suite (xn )
vers x? est dite d’ordre au moins p, p ∈ [1, +∞[, si
∃n0 ∈ N : ∀n > n0 |xn+1 − x? | 6 C|xn − x? |p ,
où C est une constante strictement positive et C < 1 si p = 1. La convergence est dite
d’ordre exactement p si on a de plus
∃n0 ∈ N : ∀n > n0 c|xn − x? |p 6 |xn+1 − x? | 6 C|xn − x? |p ,
où c est une autre constante.
Proposition 4.8 – critère d’ordre de convergence

Soit (xn )n∈N une suite qui converge vers une limite x? sans jamais être égale à x? . Si la
limite
|xn+1 − x? |
lim =K
n→∞ |xn − x? |p
existe pour un certain p ∈ [1, +∞[, alors

. si p = 1 et 0 < K < 1, la suite converge linéairement (exactement d’ordre 1),
. si p > 1 et 0 < K, la suite converge à l’ordre exactement p.
La limite K est alors appelé la constante asymptotique de l’erreur.
Pour bien comprendre cette notion d’ordre de convergence, supposons que nous avons une suite
(xn )n∈N qui converge à l’ordre p vers x? et que de plus nous avons
|xn+1 − x? |
lim = K > 0.
n→∞ |xn − x? |p
Le nombre dn = − log10 |xn − x? | mesure le nombre de chiffres (en base 10) de xn et de x? qui
coı̈ncident (à une constante additive près ne dépendant pas de n). On peut écrire
dn+1 ≈ pdn − log10 (K).
Ce qui montre qu’à une constante additive près, le nombre de chiffres exactes est multiplié par
p à chaque itération. En effet on a pour p > 1
log10 (K) log10 (K)
dn+1 + = p dn + .
1−p 1−p
On peut donc estimer le nombre d’itérations nécessaires pour gagner un chiffre exacte. En par-
ticulier, lorsque la convergence est linéaire, comme ce sera le cas pour la plupart des méthodes
⁄
que nous verrons, même celles d’ordre élevée lorsqu’elles seront confrontées à des situations
particulières comme celles de recherche des racines multiples.
Proposition 4.9 – Utilisation de la constante asymptotique d’erreur

Soit (xn )n∈N une suite qui converge linéairement avec une constante asymptotique d’erreur
égale à K. Alors le nombre d’itérations nécessaires pour gagner un chiffre exacte est le
plus petit entier supérieur à −1/ log10 (K).
Démonstration. la formule de récurrence du nombre de chiffres exactes par itération est donnée ici par
dn+1 = dn − log10 (K). Ainsi après m itérations à partir de l’itération n, on a dn+m = dn − m log10 (K).
Par conséquent, dn+m = dn + 1 si et seulement si m = −1/ log10 (K) (arrondi à l’entier supérieur). !
Quelques fois il est suffisant d’estimer l’ordre de convergence au moyen de comparaisons :
Définition 4.10
Soient (xn )n∈N , (yn )n∈N deux suites qui convergent respectivement vers x? , y ? . On dit que
−x?
(xn ) converge plus vite que (yn ) si limn→∞ xynn −y ? = 0.
Ainsi, si la suite (xn ) converge à l’ordre q, si la suite (yn ) converge à l’ordre p et si (xn ) converge
plus rapidement que (yn ), alors q > p. Par conséquent si l’on ne dispose que de l’ordre p de (yn ),
on pourra dire pour (xn ) convergeant plus vite que (yn ) qu’elle converge à l’ordre au moins p.
3 Méthodes de type encadrement
Nous commençons par présenter deux méthodes de recherche de 0 qui ne sont pas très rapides
mais qui ont l’intérêt d’être très robustes. En particulier, elles fonctionnent même pour les racines
multiples.
3.1 méthode de la dichotomie
Soit f : [a, b] → R continue telle que f (a)f (b) < 0. D’après le théorème de Rolle, f admet une
racine dans ]a, b[. Cette racine est certainement dans l’une des moitiés de l’intervalle [a, b]. C’est-
à dire si l’on pose c = a+b
2 . L’un des intervalles [a, c] ou [c, b] nous placera dans la configuration
de départ. On prend celui là et on rejette l’autre. On répète le processus, ce qui génère une
suite d’intervalles ([an , bn ]) emboı̂tés, dont la longueur tends vers 0 et qui vérifie pour tout
n, f (an )f (bn ) < 0. La méthode de dichotomie suit donc la procédure suivante :
. on pose [a0 , b0 ] = [a, b]
. [an , bn ] étant connu, on pose c = an +b
2
n
et on teste le signe de f (an )f (c)
. Si cette valeur est strictement négative, la racine de f est dans [an , c] on pose alors
[an+1 , bn+1 ] = [an , c]
. Si cette valeur est strictement positive, alors f (c)f (bn ) < 0, on pose alors [an+1 , bn+1 ] =
[c, bn ]
. Si elle est nulle alors c est la racine.
La figure 4.2 illustre le comportement de la méthode de la dichotomie.
⁄
Figure 4.2 – Illustration de la méthode de la dichotomie
Proposition 4.11 – Propriétés de la méthode de la dichotomie

Soit f continue sur [a, b] telle que f (a)f (b) < 0 et qui possède un unique zéro dans ]a, b[
noté x? . Soit (an ), (bn ) les suites générées par la méthode de dichotomie, alors
. x? ∈ [an , bn ] pour tout n ;
. [an+1 , bn+1 ] ⊂ [an , bn ] pour tout n ;
. bn − an = (b − a)/2n ;
. |an − x? | 6 (b − a)/2n et |bn − x? | 6 (b − a)/2n ;
. les suites (an ) et (bn ) converge vers x? .
Démonstration. La preuve est triviale et laissée en exercice au lecteur. !
Nous remarquons que, sous réserve que les hypothèses de la proposition précédente sont sa-
tisfaites, la méthode de la dichotomie converge toujours et il est possible d’estimer le nombre
d’itérations nécessaires pour déterminer le zéro x? avec une précision choisie. En effet, à l’étape
n, on est asuré que x? se trouve dans l’intervalle [an , bn ] de longueur (b − a)/2n . Ainsi, pour
fixé, le choix de
b − a b − a
n = E log2 + 1 = E log2
2
permet d’assurer que (an + bn )/2 est proche de x? à près.
⁄
Algorithme 4.1 – Dichotomie

def dichotomie (f , a , b , epsilon , verbose = False ):
an , bn = (a , b ) if a < b else (b , a )
fan , fbn = f ( an ) , f ( bn )
la , lb , lc = [ an ] , [ bn ] , []
if fan * fbn > 0:
raise ValueError ( f " Probleme dans dichotomie : { a } , { b } " )
while bn - an > epsilon :
cn = .5*( an + bn )
fcn = f ( cn )
if fan * fcn <= 0:
bn , fbn = cn , fcn
if fbn * fcn <= 0:
an , fan = cn , fcn
la . append ( an )
lb . append ( bn )
lc . append ( cn )
cn = .5*( an + bn )
lc . append ( cn )
if verbose :
return cn , np . asanyarray ( la ) , np . asanyarray ( lb ) , np . asanyarray ( lc )
return cn
La fonction dichotomie prend cinq arguments : la fonction f dont on cherche un zéro,

deux réels a et b qui définissent l’intervalle de recherche, un réel epsilon qui est un petit
paramètre utilisé pour arrêter l’algorithme lorsque la précision est atteinte et un booléen
verbose qui permet de modifier la sortie (seulement la solution trouvée ou bien toutes
les valeurs intermédiaires calculées).
Résumé de la méthode de dichotomie :

. avantages :
. méthode simple à implémenter,
. convergence certaine
. inconvénients :
. hypothèse de départ contraignante (encadrement du zéro),
. convergence lente.
3.2 Méthode de la sécante
On fera attention au vocabulaire : la méthode de la sécante est appelée method of false position
or regula falsi dans les publications anglophones tandis que la méthode de la fausse position est
appelée secant method.
Cette méthode suit le même principe que la méthode de dichotomie avec pour seule différence
qu’au lieu de prendre cn comme le milieu de an et de bn , on le prendra de manière équivalente
comme la racine du polynôme interpolateur de Lagrange associé au noeuds an et bn , ou comme
l’abscisse du point d’intersection avec l’axe des abscisses de la droite passant par (an , f (an )), et
(bn , f (bn )). C’est-à-dire
an f (bn ) − bn f (an )
cn = .
f (bn ) − f (an )
Comme f (an ) et f (bn ) sont toujours de signe différent (par définition de l’algorithme), le déno-
minateur ne s’annule jamais et la suite est bien définie.
⁄
On construit donc les trois suites (an ), (bn ) et (cn ) par a0 = a, b0 = b, et pour n ∈ N,
( (
an f (bn ) − bn f (an ) an si f (an )f (cn ) < 0, cn si f (an )f (cn ) < 0,
cn = , an+1 = bn+1 =
f (bn ) − f (an ) cn si f (an )f (cn ) > 0, bn si f (an )f (cn ) > 0,
(4.8)
Si f (cn ) = 0 pour un certain n ∈ N, on ne continue pas l’algorithme et les suites sont finies.
La figure 4.3 illustre le comportement de la méthode de la sécante.
Figure 4.3 – Illustration de la méthode de la sécante
Proposition 4.12 – Propriétés de la méthode de la sécante

Soit f continue sur [a, b] telle que f (a)f (b) < 0 et qui possède un unique zéro dans ]a, b[
noté x? . Soit (an ), (bn ) et (cn ) les suites générées par la méthode de la sécante, alors
. x? ∈ [an , bn ] pour tout n ;
. [an+1 , bn+1 ] ⊂ [an , bn ] pour tout n ;
. la suite (cn ) converge vers x? ;
Démonstration. Pour fixer les idées, on suppose comme à la figure 4.3 que f (a) < 0 et f (b) > 0.
Supposons que les trois suites sont infinies, c’est-à-dire que f (cn ) 6= 0 pour tout n ∈ N. La suite (an )
est croissante majorée, elle converge donc vers α, la suite (bn ) est décroissante minorée, elle converge
donc vers β. Par continuité, on a f (α) 6 0 et f (β) > 0. Deux cas se présentent alors.
⁄
Premier cas : f (α) = f (β). On a donc f (α) = f (β) = 0. Cela prouve que α = β = x? et que la suite
(cn ) tend vers x? .
Deuxième cas : f (α) 6= f (β). Par continuité de f , on en déduit que la suite (cn ) converge vers γ avec
αf (β) − βf (α)
γ= .
f (β) − f (α)
Or à chaque étape cn = an ou cn = bn . Supposons que cn = an pour une infinité de n ∈ N. Alors γ = α.

On en déduit que (α − β)f (α) = 0. Comme f (α) 6= f (β), on a α 6= β. Donc f (α) = 0, et donc α = x? ,
ce qui prouve bien que (cn ) tend vers x? . De même si cn = bn pour une infinité de n ∈ N, on montre
que β = x? et on obtient la même conclusion. !
On a donc prouvé que la suite (cn ) converge vers x? mais on ne sait rien des deux suites (an )
et (bn ). On sait seulement qu’une des deux converge vers x? . On peut également montrer que,
si la fonction est strictement convexe ou strictement concave, l’une des deux suites (an ) ou (bn )
reste constante.
Proposition 4.13 – Convergence de la méthode de la sécante

Soit f ∈ C 2 ([a, b]) telle que f (a)f (b) < 0 et f 00 n’a aucune racine dans l’intervalle [a, b]
(la fonction est strictement convexe ou strictement concave). Soit (an ), (bn ) générées par
la méthode de fausse position. Deux cas sont alors possibles :
. la suite (bn ) est constante, alors (an ) converge linéairement vers la racine x∗ de f
et on a
x? − an+1 ?
0 ? x −b
K1 = lim = 1 + f (x ) ;
n→∞ x? − an f (b)
. la suite (an ) est constante, alors (bn ) converge linéairement vers la racine x? de f
et on a
x? − bn+1 ?
0 ? x −a
K1 = lim = 1 + f (x ) .
n→∞ x? − bn f (a)
Démonstration. Pour fixer les idées, on suppose comme à la figure 4.3 que f (a) < 0 et f (b) > 0 et que
f 00 (x) > 0 pour x ∈ [a, b]. Les autres cas se traitent de manière identique. Comme la fonction f est
strictement convexe, la courbe est toujours sous ses sécantes, c’est-à-dire que f (cn ) < 0, la suite (bn )
est donc constante.
Nous avons donc une seule suite à étudier, la suite (an ) définie par
xf (b) − bf (x)
a0 = a, an+1 = φ(an ) avec φ(x) = .
f (b) − f (x)
Comme f (x? ) = 0, nous avons φ(x? ) = x? . Nous avons donc
an+1 − x? = φ(an ) − φ(x? ) = φ0 (ξ)(an − x? ),
où ξ est dans l’intervalle ]an , x? [. D’après la proposition précédente, la suite (an ) converge vers x? car
elle est égale à la suite (cn ). Nous avons donc
x? − an+1
lim = φ0 (x? ).
n→∞ x? − an
Dérivons la fonction φ. Nous avons

(f (b) − bf 0 (x))(f (b) − f (x)) + f 0 (x)(xf (b) − bf (x))
φ0 (x) =
(f (b) − f (x))2
f (b)(f (b) − f (x) − (b − x)f 0 (x))
= .
(f (b) − f (x))2
Cette expression se simplifie pour x = x? en
f (b) − (b − x? )f 0 (x? )
φ0 (x? ) = ,
f (b)
⁄
Cette proposition ne permet pas facilement d’estimer à l’avance le nombre d’itérations nécessaires
pour que la méthode de la sécante donne un résultat cn proche de x? à près pour > 0 fixé.
En effet, la limite K1 ne peut pas être calculée puisqu’elle dépend de x? inconnu.
Dans le cas de la figure 4.3 obtenue avec f (x) = e2x −4, nous avons x? = ln(2). Ainsi, pour b = 1,
K1 ' 0.2757. La convergence est donc plus rapide que celle obtenue pour la dichotomie. Pour
b = 1.6, K1 ' 0.6467. La convergence est donc plus lente que celle obtenue pour la dichotomie.
Algorithme 4.2 – Sécante

def secante (f , a , b , epsilon , verbose = False ):
an , bn = (a , b ) if a < b else (b , a )
fan , fbn , fcn = f ( an ) , f ( bn ) , 2* epsilon
la , lb , lc = [ an ] , [ bn ] , []
if fan * fbn > 0:
raise ValueError ( f " Probleme dans secante : { a } , { b } " )
while abs ( fcn ) > epsilon :
cn = ( an * fbn - bn * fan )/( fbn - fan )
fcn = f ( cn )
if fan * fcn <= 0:
bn , fbn = cn , fcn
if fbn * fcn <= 0:
an , fan = cn , fcn
la . append ( an )
lb . append ( bn )
lc . append ( cn )
cn = ( an * fbn - bn * fan )/( fbn - fan )
lc . append ( cn )
if verbose :
return cn , np . asanyarray ( la ) , np . asanyarray ( lb ) , np . asanyarray ( lc )
return cn
La fonction secante prend cinq arguments : la fonction f dont on cherche un zéro, deux
réels a et b qui définissent l’intervalle de recherche, un réel epsilon qui est un petit
paramètre utilisé pour arrêter l’algorithme lorsque la précision est atteinte et un booléen
verbose qui permet de modifier la sortie (seulement la solution trouvée ou bien toutes
les valeurs intermédiaires calculées).
Résumé de la méthode de la sécante :

. avantages :
. méthode simple à implémenter,
. convergence certaine et souvent meilleure que pour la dichotomie ;
. inconvénients :
. hypothèse de départ contraignante (encadrement du zéro),
. complexité un peu supérieure à celle de la dichotomie,
. convergence lente et très lente vers les racines multiples.
4 Méthodes de type interpolation
Nous présentons à présent deux méthodes de recherche de 0 qui sont plus efficaces. Le principe
est de remplacer la fonction f par un polynôme interpolateur de bas degré qui “ressemble” à la
fonction f lorsque l’on s’approche du 0 recherché.
⁄
4.1 méthode de Newton
La méthode de Newton est certainement une des méthodes les plus utilisées pour chercher les
zéros d’une fonction. Elle consiste à remplacer la fonction f dont on cherche le zéro par sa
tangente (il faut donc que la fonction f soit dérivable).
Etant donné un point x0 , nous calculons une nouvelle approximation x1 en cherchant le zéro de
la tangente à f au point x0 . Cette tangente est la droite Tf (x0 ) d’équation
Tf (x0 ) : y = f (x0 ) + (x − x0 )f 0 (x0 ).
Cette droite coupe l’axe des abscices en x1 = x0 − f (x0 )/f 0 (x0 ). Cela permet de construire la
suite (xn ) tant que la dérivée f 0 (xn ) ne s’annule pas :
f (xn )
xn+1 = xn − , n > 0.
f 0 (xn )
Lorsque la fonction f a un domaine de définition qui n’est pas R, il est possible qu’un terme xn
sorte du domaine de définition. La construction de la suite échoue, la méthode ne donne pas de
résultat intéressant... La conclusion est qu’il est nécessaire de bien choisir le point de départ x0
de la suite pour assurer que la suite est constructible et qu’elle converge bien vers x? .
La figure 4.4 illustre le comportement de la méthode de Newton dans un cas favorable.
Proposition 4.14 – Convergence de la méthode de Newton (racine simple)

Soit f une fonction de classe C 1 sur [a, b] qui possède un unique zéro dans ]a, b[ noté x? et
qui vérifie f 0 (x? ) 6= 0. Alors il existe h > 0 tel que, si x0 ∈ [x? − h, x? + h], alors la suite
(xn ) définie par
f (xn )
xn+1 = xn − , n > 0,
f 0 (xn )
est constructible et converge vers x? .

Si de plus f est de classe C 2 , la convergence est quadratique : nous avons en effet
xn+1 − x? f 00 (x? )
lim = .
n→∞ (xn − x? )2 2f 0 (x? )
Démonstration. Supposons pour simplifier que f 0 (x? ) = α > 0. L’autre cas se démondre de la même
façon. Comme la fonction f 0 est continue autour de x? ,
∀c < 1 < c, ∃h > 0 : |x − x? | 6 h =⇒ αc 6 f 0 (x) 6 αc.
Le choix des constantes sera précisé plus loin. Nous avons donc
f (xn )
xn+1 − x? = xn − x? − .
f 0 (xn )
Or, un développement limité de f au point xn donne
0 = f (x? ) = f (xn ) + (x? − xn )f 0 (ξ),
⁄
Figure 4.4 – Illustration de la méthode de Newton
pour ξ dans l’intervalle entre x? et xn . Nous obtenons donc

f 0 (ξ)
xn+1 − x? = (xn − x? ) 1 − 0
f (xn )
Or )
αc 6 f 0 (ξ) 6 αc c−c f 0 (ξ) c−c f 0 (ξ) c−c
=⇒ 61− 0 6 =⇒ 1 − 0 6 .
αc 6 f 0 (xn ) 6 αc c f (xn ) c f (xn ) c
Prenons par exemple c = 9/10 et c = 11/10, nous concluons que

2
|xn+1 − x? | 6 |xn − x? |.
9
Ainsi, quitte à diminuer h, on peut supposer que tous les termes de la suite sont dans l’intervalle
[x? − h, x? + h] et une récurrence immédiate donne |xn − x? | 6 (2/9)n h, ce qui permet de conclure à la
convergence de la suite.
Si de plus f ∈ C 2 ([a, b]) alors le développement limité autour du point xn peut être prolongé en
0 = f (x? ) = f (xn ) + (x? − xn )f 0 (xn ) + 21 (x? − xn )2 f 00 (ξ),
pour ξ dans l’intervalle entre x? et xn . On en déduit que
1 f 00 (ξ)
xn+1 − x? = (xn − x? )2 0 ,
2 f (xn )
ce qui permet de conclure après passage à la limite. !
⁄
La propriété précédente assure la convergence de la méthode seulement lorsque le point de départ

x0 est “suffisamment” proche du zéro recherché (avec des hypothèses de régularité). Il n’est pas
possible d’améliorer dans le cas général ce résultat comme on peut le voir grâce à l’illustration
proposée à la figure 4.5 :
. lorsque x0 est à gauche de la racine x? , la concavité de la fonction assure que la suite (xn )
converge vers x? ;
. lorsque x0 est proche de 0 ou de 1, la suite converge vers un 2-cycle, c’est-à-dire que les
valeurs oscillent entre les deux valeurs ;
. lorsque x0 est dans la partie convexe (x0 = 2 par exemple), le comportement n’est pas
prévisible : sur le dessin, la suite commence par osciller autour d’un 2-cycle puis converge
brutalement vers x? ...
Figure 4.5 – Illustration de différents comportement de la méthode de Newton
Algorithme 4.3 – Newton

def newton (f , df , x0 , epsilon , verbose = False ):
x = x0
fx , dfx = f ( x ) , df ( x )
lx = [ x ]
n , Nmax = 0 , 100
while abs ( fx ) > epsilon * min (1 , abs ( dfx )) and n < Nmax :
n += 1
x -= fx / dfx
fx , dfx = f ( x ) , df ( x )
lx . append ( x )
if verbose :
return x , np . asanyarray ( lx )
return x
La fonction newton prend cinq arguments : la fonction f dont on cherche un zéro, sa

dérivée df, un réel x0 qui sera le premier terme de la suite, un réel epsilon qui est un
petit paramètre utilisé pour arrêter l’algorithme lorsque la précision est atteinte et un
booléen verbose qui permet de modifier la sortie (seulement la solution trouvée ou bien
toutes les valeurs intermédiaires calculées).
Résumé de la méthode de Newton :

. avantages :
. méthode extrêmement rapide (toujours plus rapide que la sécante) ;
⁄
. hypothèse moins contraignante car il n’est pas nécessaire d’avoir un encadrement du

zéro ;
. inconvénients :
. évaluation obligatoire de la dérivée ;
. convergence non assurée et lente pour les racines multiples.
4.2 méthode de la fausse position
L’idée de la méthode de la fausse position est de remplacer f 0 (xn ) dans la méthode de Newton
par autre chose car il est rare de connaı̂tre en pratique une expression analytique de f 0 et on ne
peut en avoir qu’une approximation. Cela peut arriver par exemple lorsque la fonction f est le
résultat d’un calcul informatique ou bien une reconstruction à partir de données expérimentales.
On remplace f 0 (xn ) dans la méthode de Newton par
f (xn ) − f (xn−1 )
.
xn − xn−1
La méthode s’écrit ainsi : on choisit x0 = a et x1 = b ou l’inverse et on pose pour n ∈ N
xn − xn−1
xn+1 = xn − f (xn ) . (4.9)
f (xn ) − f (xn−1 )
C’est-à-dire que lorsque xn et xn−1 sont construits, on construit xn+1 comme l’intersection de
l’axe des abscisses et de la droite reliant les points f (xn ) et f (xn−1 ).
La figure 4.6 illustre le comportement de la méthode de la fausse position dans un cas favorable.
Il n’y a pas toujours constructiblité de cette méthode. Sans hypothèses supplémentaires, il se

peut que l’on sorte de l’intervalle de définition de la fonction ou bien que le dénominateur soit
nul. On a toutefois un résultat de convergence locale.
Proposition 4.15 – Convergence de la méthode de la fausse position

On suppose que f est de classe C 2 sur [a, b], qu’il existe un unique zéro de f dans [a, b], noté
x? , et que f 0 (x? ) 6= 0. Il existe alors > 0 tel que pour tout x0 , x1 tels que |x0 − x? | <
et |x1 − x? | < , la suite (xn ) converge vers x? et√l’erreur |xn − x? | est majorée par une
suite qui converge vers 0 à l’ordre au moins q = ( 5 + 1)/2.
Démonstration. La démonstration utilise des outils élémentaires d’analyse mais est un peu technique.
Le lecteur intéressé pourra la trouver dans [5]. !
⁄
Figure 4.6 – Illustration de la méthode de la fausse position
Algorithme 4.4 – Fausse position
def fausse_position (f , x0 , x1 , epsilon , verbose = False ):

x , xold = x1 , x0
fx , fxold = f ( x ) , f ( xold )
dfx = ( fx - fxold )/( x - xold )
lx = [ xold , x ]
n , Nmax = 0 , 100
while abs ( fx ) > epsilon * min (1 , abs ( dfx )) and n < Nmax :
n += 1
xold = x
x -= fx / dfx
fxold , fx = fx , f ( x )
dfx = ( fx - fxold )/( x - xold )
lx . append ( x )
if verbose :
return x , np . asanyarray ( lx )
return x
La fonction fausse_position prend cinq arguments : la fonction f dont on cherche un

zéro, deux réels x0 et x1 qui seront les deux premiers termes de la suite, un réel epsilon
qui est un petit paramètre utilisé pour arrêter l’algorithme lorsque la précision est atteinte
et un booléen verbose qui permet de modifier la sortie (seulement la solution trouvée ou
bien toutes les valeurs intermédiaires calculées).
⁄
Résumé de la méthode de la fausse position :

. avantages :
. méthode extrêmement rapide (presque aussi rapide que la méthode de Newton) ;
. hypothèse moins contraignante car il n’est pas nécessaire d’avoir un encadrement du
zéro ;
. pas d’évaluation de la dérivée ;
. inconvénients :
. le calcul de la dérivée approchée peut entrainer des erreurs d’arrondis importantes ;
. convergence non assurée et lente pour les racines multiples.
⁄
Bibliographie
[1] G. Allaire. Analyse numérique et optimisation. Edition de l’Ecole Polytechnique, 2002.
[2] H. Brezis. Analyse Fonctionnelle Théorie et Applications. Masson, 1983.
[3] P. Ciarlet. Introduction à l’analyse numérique matricielle et à l’optimisation. Masson, 1994.
[4] M. Crouzeix and A. L. Mignot. Analyse numérique des équations différentielles. Masson,
1992.
[5] J. Demailly. Analyse Numérique et Équations différentielles. EDP Sciences, 2006.
[6] G. Duvaut. Mécanique des milieux continus. Masson, 1990.
[7] D. Euvrard. Résolution numérique des EDP. Masson, 1990.
[8] L. Evans. Partial Differential Equations, volume 19 of Graduate Studies in Mathematics.

American Mathematical Society, 2002.
[9] E. Godlewski and P.-A. Raviart. Hyperbolic systems of conservation laws. Ellipses, 1991.
[10] L. Hörmander. Lectures on Nonlinear Hyperbolic Differential Equations. Springer, 1997.
[11] F. Hubert and J. Hubbard. Calcul scientifique de la théorie à la pratique. Vuibert, 2006.
[12] B. Lucquin. Équations aux dérivées partielles et leurs approximations. Ellipses, 2004.
[13] A. L. Pourhiet. Résolution numérique des EDP. Impr. du sud, 1988.
[14] P.-A. Raviart and J.-M. Thomas. Introduction à l’analyse numérique des équations aux
dérivées partielles. Masson, 1983.
[15] M. Renardy and R. Rogers. An introduction to partial differential equations. Springer, 1993.
[16] R. Richtmyer and K. Morton. Difference methods for initial value problems. Wiley-
Interscience, New-York, 1967.
[17] L. Schwartz. Méthodes mathématiques pour les sciences physiques. Hermann, 1965.
[18] L. Schwartz. Théroie des distributions. Hermann, 1966.
[19] D. Serre. Matrices : Theory and Applications. Springer, 2002.
[20] D. Serre. Systems of conservation laws 1. Cambridge University Press, 2003.
⁄

Programmation Analyse Num Erique: Licence 2 - Maths Info

Transféré par

Droits d'auteur :

Formats disponibles

Programmation Analyse Num Erique: Licence 2 - Maths Info

Transféré par

Informations du document

Titre original

Copyright

Formats disponibles

Partager ce document

Partager ou intégrer le document

Options de partage

Avez-vous trouvé ce document utile ?

Ce contenu est-il inapproprié ?

Droits d'auteur :

Formats disponibles

Programmation Analyse Num Erique: Licence 2 - Maths Info

Transféré par

Droits d'auteur :

Formats disponibles

Université de la Polynésie française

Licence 2 - Maths Info

4 Résolution d’équations ordinaires 45

1 Nécessité de l’analyse numérique

u0 (t) = f (t, u(t)), t > 0,

Nous pouvons alors récrire le problème sous la forme

un+1 = un + φ(un , un+1 ),

Définition 1.1 – continuité

∀ε > 0 ∃η > 0 : y ∈ I ∩ ]x−η, x+η[ =⇒ |f (x) − f (y)| 6 ε.

Théorème 1.2 – Valeurs intermédiaires

Définition 1.3 – dérivabilité

Théorème 1.4 – Rolle

Théorème 1.5 – Accroissements finis

f (b) − f (a) = f 0 (c)(b − a).

2.2 Formules de Taylor

Théorème 1.7 – Taylor avec reste intégral

Théorème 1.8 – Taylor – Lagrange

Théorème 1.9 – Taylor – Young

où la fonction Rn est négligeable devant la fonction h 7→ hn au voisinage de 0. C’est-à-dire

fonctions régulières, il s’agit de cas générique.

Définition 1.10 – fonction convexe

∀(x, y) ∈ [a, b]2 , x 6= y,

La fonction sera dite strictement convexe si l’inégalité est stricte.

Finalement une fonction f sera dite concave si −f est convexe.

Proposition 1.11 – convexité pour les fonctions régulières

Comme f 0 est croissante,

∀s ∈ [0, 1] sy + (1 − s)z > sy + (1 − s)x =⇒ f 0 (sy + (1 − s)z) > f 0 (sy + (1 − s)x).

f (y) − f (z) > λf (y) − λf (x) ⇐⇒ f (z) 6 λf (x) + (1 − λ)f (y),

ce qui termine la preuve en reprenant l’expression de z = λx + (1 − λ)y. !

Interpolation ou approximation Interpolation ou approximation

0.0 points 0.0 points

Interpolation ou approximation Interpolation ou approximation

0.0 points 0.0 points

Interpolation ou approximation Interpolation ou approximation

0.0 points 0.0 points

2 Exemples avec peu de points

2.1 Interpolation à un seul point

Figure 2.2 – Exemple de polynôme interpolateur de Lagrange lorsque N = 1.

2.2 Interpolation à deux points : exemple de la droite

2.2.1 Un peu de géométrie euclidienne

Théorème 2.1 – Polynôme interpolateur de degré 1

(x, y) ∈ D ⇐⇒ (x − x1 )(y2 − y1 ) − (y − y1 )(x2 − x1 ) = 0.

Si de plus x1 6= x2 , il existe un unique polynôme P ∈ R1 [X] tel que

Le calcul de ce déterminant donne directement l’expression de la droite D. Puis, dans le cas x1 6= x2 , il

2.2.2 Une première méthode de calcul

Nous avons donc

Ainsi, si x1 et x2 sont proches, le calcul de la matrice A−1 et plus précisément la résolution

2.2.3 Expression du polynôme dans une base adaptée

2.3 Interpolation à trois points

Figure 2.3 – Exemple de polynôme interpolateur de Lagrange lorsque N = 3.

Un calcul un peu fastidieux permet de déterminer la valeur des trois coefficients a, b et c :

Plusieurs remarques peuvent alors être faites :

Exercice 2.2 – Polynôme interpolateur de degré 2

Solution. Nous écrivons P = a0 + a1 (X − x1 ) + a2 (X − x1 )(X − x2 ). Nous avons en évaluant en x1

3 Polynôme interpolateur de Lagrange

3.1 Définition et propriétés

Théorème 2.3 – Isomorphisme d’espaces vectoriels

est bijective (c’est un isomorphisme).

Définition 2.4 – Polynôme interpolateur de Lagrange

a+b b−a (2j − 1)π