0% ont trouvé ce document utile (0 vote)

24 vues66 pages

Initiatio A Lanalyse Numerique

Cet ouvrage présente une initiation à l'analyse numérique, abordant des concepts fondamentaux tels que l'approximation de solutions d'équations, l'interpolation, et le calcul numérique des intégrales. Les méthodes classiques comme la dichotomie, les approximations successives, et la méthode de Newton sont expliquées en détail. Le texte est destiné aux étudiants de premier cycle universitaire et ne traite pas des méthodes avancées ou des équations aux dérivées partielles.

Transféré par

daniela6.belak9

Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.

Formats disponibles

Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd

0% ont trouvé ce document utile (0 vote)

24 vues66 pages

Initiatio A Lanalyse Numerique

Transféré par

daniela6.belak9

Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.

Formats disponibles

Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd

Initiation à l’Analyse Numérique

Robert Rolland
C.N.R.S. Institut de Mathématiques de Luminy, case 930, F13288
Marseille cedex 9, France.
E-mail address: [email protected]
Table des matières

Préface vii
Chapitre 1. Approximations de solutions d’équations 1
1.1. Introduction 1
Chapitre 2. Les outils classiques d’approximation 7
2.1. Introduction 7
2.2. La formule de Taylor 9
2.3. La formule d’Euler-Maclaurin 9
Chapitre 3. Interpolation des fonctions 15
3.1. Introduction 15
3.2. Interpolation de Lagrange 15
3.3. Le problème général de l’interpolation 25
3.4. Quelques exemples importants 25
Chapitre 4. Calcul numérique des intégrales définies 33
4.1. Introduction 33
4.2. Mise en œuvre de méthodes interpolatoires 33
4.3. Présentation générale des quadratures élémentaires - Ordre d’une
méthode 45
4.4. Accélération de convergence - Méthode de Romberg 45
Chapitre 5. Analyse numérique des équations différentielles 49
5.1. Introduction 49
5.2. Généralités sur les méthodes 50
Annexe A. Représentation des nombres 57
A.1. Introduction 57
A.2. Les nombres réels 57
A.3. Les entiers 59

v
Préface

Cet ouvrage est une initiation à l’analyse numérique. Les notions qui y sont
exposées et développées concernent les objets et méthodes de base de ce vaste
domaine des mathématiques. Le niveau est celui des deux ou trois premières années
d’Université. Nous n’abordons pas les méthodes techniquement très élaborées. Nous
n’abordons pas non plus l’immense domaine des équations aux dérivées partielles
et équations fonctionnelles générales faisant intervenir les résultats généraux sur les
espaces fonctionnels.
Nous essayons de montrer comment l’algèbre linéaire d’une part et quelques
outils classiques de l’analyse d’autre part, permettent de poser clairement les
problèmes d’approximation, d’interpolation, de calculs d’intégrales et de solu-
tions d’équations différentielles. Nous donnons les méthodes numériques classiques
concernant ces questions.

Robert Rolland

vii
CHAPITRE 1

Approximations de solutions d’équations

1.1. Introduction
Nous voulons résoudre numériquement l’équation :
(1.1) f (x) = 0
où f est une fonction suffisamment régulière. Nous supposerons f au moins
continue, parfois de classe C 1 ou plus. Il existe de nombreuses méthodes sophis-
tiquées pour attaquer ce problème. Nous donnons ici deux idées importantes à
la base de la plupart de ces méthodes : d’une part la dichotomie, d’autre part
l’approximation de Newton. Ces deux idées seront illustrées par deux méthodes
correspondantes.
1.1.1. Dichotomie. Nous supposons que nous avons isolé le zéro c que nous
voulons calculer, c’est-à-dire que nous avons trouvé un intervalle [a,b] tel que :

c ∈]a, b[,
∀x ∈ [a, b] \ {c}, f (x) 6= 0
Nous supposerons en outre que f change de signe en c, donc que :
f (a)f (b) < 0.
La méthode de calcul de c par dichotomie est très simple : on coupe l’intervalle sur
lequel on travaille en deux et on teste sur lequel des deux sous-intervalles obtenus
se trouve c. On réitère le procédé à partir du sous-intervalle déterminé. Après n
itérations on localise le nombre c sur un intervalle de longueur (b − a)/2n . Voici les
détails :
On fixe ǫ > 0,

A := a ;
B := b ;
E := ǫ ;
tant que B − A > E faire
C := (A + B)/2 ;
si f (C) = 0 alors
retourner C ;
sortir ;
finsi ;
si f (A)f (C) < 0 alors
B := C ;
sinon
A := C ;
1
2 1. APPROXIMATIONS DE SOLUTIONS D’ÉQUATIONS

finsi ;
fintq ;
retourner C ;

1.1.2. Approximations successives, point fixe. Avant de présenter la

méthode de Newton, disons quelques mots sur la méthode générale des approxi-
mations successives. Soit f une fonction de classe C 1 définie sur un intervalle [a, b]
à valeurs dans l’intervalle [a, b]. Soit u0 ∈ [a, b]. On définie alors par récurrence la
suite de terme général ui en posant pour tout i ≥ 1 ui = f (ui−1 ). On étudie la
convergence de cette suite. Nous allons montrer que sous certaines conditions on
peut affirmer qu’il existe un point fixe pour f , c’est-à-dire une solution de l’équation
f (x) = x (cf. figure 1).

x
a x0 b

Fig. 1. Théorème du point fixe

Théorème 1.1 (Théorème du point fixe). Sous les hypothèses précédentes, si

on suppose que :
sup |f ′ (x)| = k < 1,
x∈[a,b]

alors la fonction f (x) admet un point fixe x0 et un seul sur l’intervalle [a, b]. La
suite (ui )i définie précédemment converge vers le point fixe x0 .
Démonstration. D’après le théorème des accroissements finis :
|u2 − u1 | = |f (u1 ) − f (u0 )| ≤ k|u1 − u0 |,
|u3 − u2 | = |f (u2 ) − f (u1 )| ≤ k|u2 − u1 | ≤ k 2 |u1 − u0 |,
et en itérant ce calcul :
|un − un−1 | ≤ k n−1 |u1 − u0 |.
1.1. INTRODUCTION 3

On en déduit que si n ≥ p :
k p (1 − k n−p )
|un − up | ≤ |u1 − u0 |.
1−k
Donc la suite (un )n est une suite de Cauchy, elle converge vers un point x0 . Ce
point est un point fixe, en effet, un+1 = f (un ), comme la suite de terme général un
converge vers x0 et que la fonction f est continue, on conclut que f (un ) converge
vers f (x0 ) et donc que x0 = f (x0 ). On peut voir que le point fixe est unique en
constatant que si s est un point fixe alors :
|un − s| = |f (un−1 ) − f (s)| ≤ k|un−1 − s|,
et en itérant ce calcul :
|un − s| ≤ k n |u0 − s|.
Ceci prouve en effet que s est la limite de la suite de terme général un , et donc que
s = x0 .
Exemple 1.2 (Construction d’une table trigonométrique). Les propriétés géo-
métriques des lignes trigonométriques permettent de calculer facilement sin(π/6),
sin(π/4), sin(π/3), sin(π/5) ainsi évidemment que les cosinus des mêmes angles.
En utilisant la formule qui donne le sinus d’une différence on trouve sin(π/30),
en utilisant la formule qui donne le sinus de l’arc moitié on obtient sin(π/60). On
a donc une table donnant les sinus (et les cosinus) de 3◦ en 3◦ . Il faudrait donc
calculer sin(π/180) pour avoir une table trigonométrique donnant les lignes de tous
les angles en degré entier. On utilise alors la formule :
sin(3x) = 3 sin(x) − 4 sin3 (x).
On cherche alors à résoudre cette équation connaissant la valeur a = sin(3x).
L’équation se ramène à :
1
sin(x) = (4 sin3 (x) + a).
3
Donc si on pose :
1
f (u) = (4u3 + a),
3
on est amené à chercher un point fixe de f . Ici u est proche de 0 (sin(π/180)), donc la
dérivée est dans un voisinage du point fixe bien plus petite que 1 en valeur absolue,
et le théorème précédent s’applique. On peut prendre par exemple u0 = a/3.
1.1.3. Méthode de Newton. Il s’agit ici de remplacer la fonction f dont on
cherche un zéro par sa tangente en un point voisin du zéro cherché (cf. figure 2).
Ainsi à partir d’un point u0 proche de la solution x0 de l’équation f (x) = 0
(qu’on supposera unique, tout au moins dans un intervalle adapté), on construit
la tangente à la courbe y = f (x) au point (u0 , f (u0 ). Cette tangente coupe l’axe
des abscisses au point u1 . On itère cette construction à partir de la valeur u1 pour
obtenir le point u2 . Le calcul de l’équation de la tangente au point d’abscisse x et
de son intersection avec l’axe des x montre que si on introduit :
f (x)
F (x) = x − ,
f ′ (x)
alors on peut écrire :
u1 = F (u0 ), u2 = F (u1 ), · · · , un = F (un−1 ), · · ·
4 1. APPROXIMATIONS DE SOLUTIONS D’ÉQUATIONS

x
u1 u0

Fig. 2. Méthode de Newton

On vérifie imédiatement que x0 est point fixe de la fonction F (x). De plus si on

suppose la fonction f de classe C 2 par exemple et si on suppose que la dérivée f ′ de f
ne s’annule pas, alors F (x) est dérivable et sa dérivée est nulle au point x0 . Il y aura
donc un voisinage fermé de ce point fixe, pas nécessairement simple à déterminer
effectivement, dans lequel la théorie du paragraphe précédent s’applique. De plus,
comme la dérivée de F sera proche de zéro, on peut espérer une bonne convergence
de la suite (un )n , d’autant plus rapide qu’on va se rapprocher du point fixe. Nous
n’en dirons pas plus dans le cas général, nous allons développer deux exemples.
1.1.3.1. Exemple 1. Considérons la fonction :
f (x) = x2 − 2.
En nous restreignant à x ≥ 0, nous √ allons déterminer la solution positive de
l’équation x2 = 2, c’est à dire x = 2. Introduisons donc, comme nous l’indique la
méthode de Newton, la fonction :

x2 − 2 x 1 1 2
F (x) = x − =x− + = x+ .
2x 2 x 2 x
On voit que par exemple F ([1, 2]) ⊂ [1, 2] et sur cet intervalle la dérivée de F (x)
est en valeur absolue majorée par 1/2. Donc le théorème du point fixe s’applique
sur cet intervalle. Nous avons successivement :

1 2 √ √ 1 2 √ √
un + − 2 = F (un ) − F ( 2) = un + ( 2)2 − 2un 2 ,
2 un 2un

1 2 √ 1 √ 2
un + − 2= un − 2 ,
2 un 2un
1.1. INTRODUCTION 5

1 2 √ 1 √ 2
un + − 2 ≤ un − 2 ,
2 un 2
c’est-à-dire :
√ 1 √ 2
|un+1 − 2| ≤ un − 2 .
2
La convergence est donc quadratique.
1.1.3.2. Exemple 2. Considérons la fonction :
1
f (x) = a − ,
x
où a > 0, que nous étudierons pour x > 0. En appliquant la méthode de Newton
à l’équation f (x) = 0, nous obtiendrons une façon de calculer 1/a. Introduisons la
fonction :
F (x) = x(2 − ax).
Soit I l’intervalle ouvert ]0, 2/a[. Si u0 ∈ I alors 0 < F (u0 ) ≤ 1/a. Posons ensuite
u1 = F (u0 ) et par récurrence un = F (un−1 ). La suite (un )n≥1 est croissante,
majorée par 1/a, donc converge vers le point fixe 1/a de F . En renumérotant la
1
suite on peut supposer par exemple que |u1 − 1/a| ≤ 10a . Nous pouvons écrire :
|aun+1 − 1| = |2aun − a2 u2n − 1| = (aun − 1)2 .
Donc : 2n
n 1
|aun+1 − 1| = (au1 − 1)2 ≤ ,
10
ou encore : 2n
1 1 1
un+1 − ≤ .
a a 10
Là encore nous avons obtenu une convergence quadratique.
CHAPITRE 2

Les outils classiques d’approximation

2.1. Introduction
Lorsqu’on souhaite approcher une fonction, ses valeurs en des points particu-
liers, des sommes de séries, ou des valeurs d’intégrales etc. on est amené à utiliser
les outils de l’analyse et en particulier le calcul intégral. Nous verrons par la suite
diverses méthodes d’approximation d’objets divers de l’analyse, en particulier l’in-
terpolation en des points bien choisis par des fonctions élémentaires tout aussi bien
choisies. Ici nous commençons par un outil très simple mais particulièrement effi-
cace : l’intégration par partie. Si l’on veut bien laisser de côté le volet anecdotique
du calcul exact de certaines primitives, il faut comprendre l’intégration par partie
comme l’écriture d’une expression sous forme d’une partie principale (la partie tout
intégrée) et d’un reste (la deuxième intégrale)
Z x Z x
f ′ (t)g(t)dt = [f (t)g(t)]xa − f (t)g ′ (t)dt.
a a

Évidemment, si on espère que la partie

Z x
− f (t)g ′ (t)dt
a

puisse être considérée comme un reste, il faut qu’elle soit négligeable devant la
partie tout intégrée
[f (t)g(t)]xa .
Cette remarque nous indique comment choisir, quand on fait une intégration par
partie d’un produit de deux fonctions, celle qu’on intègre et celle qu’on dérive : celle
qu’on dérive est en général celle qui varie peu, de manière à obtenir une dérivée
petite. Attention ceci ne s’applique pas pour l’utilisation de l’intégration par partie
pour des calculs de primitives, ni pour l’établissement de formules théoriques où
l’on veut obtenir une forme particulière pour le terme tout intégré.
Dans la suite de ce chapitre, nous donnons deux applications très impor-
tantes de l’intégration par partie : la formule de Taylor et la formule d’Euler-
Maclaurin. Ces deux formules constituent des outils de base des méthodes d’ap-
proximation. Nous exprimerons toutes ces formules asymptotiques en utilisant des
restes intégraux qui donnent des formules exactes, et en évitant au maximum les
formules faisant intervenir des restes écrits avec des points dont on ne connait pas la
valeur mais juste une localisation. En général, dans un vrai problème on n’a jamais
vraiment besoin de ces formules et les formules avec reste intégral ainsi que l’outil
”intégration par partie” permettent d’obtenir les évaluations dont on a besoin.

2.1.1. Intégrons dans le bon sens.

7
8 2. LES OUTILS CLASSIQUES D’APPROXIMATION

Exemple 1. Il s’agit d’évaluer au voisinage de +∞ l’intégrale

Z +∞ −u
e
2
du.
ln(x) u

On choisit clairement de dériver la fonction 1/u2 et d’intégrer la fonction e−u , on

obtient
Z +∞ −u Z +∞ −u
e 1 e
2
du = 2
− 2 3
du.
ln(x) u x(ln(x)) ln(x) u

Cette manière de faire permet effectivement d’obtenir un reste négligeable devant

la partie intégrée. En effet :
−u
e−u e
3
=o ,
u u2
donc : !
Z +∞ Z +∞
e−u e−u
du = o du
ln(x) u3 ln(x) u2
et en conséquence :
Z +∞
e−u 1
2
du ∼ .
ln(x) u x(ln(x))2
Exemple 2 (communiqué par Raymond Raynaud). Il s’agit de calculer
Z e
Ip = x2 (ln(x))p dx.
1
3
On a I0 = (e − 1)/3.
2.1.1.1. Le piège. Le piège consiste ici à voir apparaı̂tre une formule très simple
et à se laisser aller à intégrer par partie ”à l’envers”. C’est-à-dire qu’on va dériver
(ln(x))p et intégrer x2 .
e3 p
Ip = − Ip−1 .
3 3
Si nous itérons le calcul nous sommes amenés à écrire la formule exacte

e3 p p(p − 1) p! p!
Ip = 1− + + · · · + (−1)p−1 p−1 + (−1)p p−1 (e3 − 1).
3 3 9 3 3
Hélas, comme à chaque étape on a pris un ”reste plus grand que la partie qui aurait
dû être principale”, cette formule est très instable numériquement. En effet si on
change un peu la valeur initiale I0 en J0 , alors le terme Jp calculé vérifie
p!
Jp − Ip = (−1)p (J0 − I0 ),
3p
ce qui fait que si J0 6= I0 , alors |Jp − Ip | → +∞. Donc une erreur d’arrondi va
complètement modifier le calcul.
2.1.1.2. L’intégration dans le bon sens. Bien sûr, si on intègre dans le bon sens
en exhibant une partie principale et un reste plus petit que cette partie principale,
p
alors ceci ne se produit plus. Intégrons donc (ln(x))
x et dérivons x3 . Nous obtenons
e3 3
Ip = − Ip+1 ,
p+1 p+1
2.3. LA FORMULE D’EULER-MACLAURIN 9

et cette fois-ci en itérant le calcul on tombe sur une formule stable qui nous permet
d’écrire tout de suite
e3
Ip ∼ .
p
Si on veut pousser le développement plus loin on obtient :

3 1 4 1
Ip = e − +o .
p p2 p2
Bien sûr la formule de récurrence trouvée est la même que dans le premier
calcul, écrite différemment. Mais justement ceci nous permet de voir que si on a
une vision ”à l’envers” de l’intégration par partie, alors la suite des calculs qu’on
est amené naturellement à faire conduit à de mauvaises situations.

2.2. La formule de Taylor

Théorème 2.1. Soit f une fonction à valeurs réelles définie sur le segment
[a − ǫ, a + ǫ] et n + 1 fois continument dérivable sur ce segment. Alors pour tout
point x du segment [a − ǫ, a + ǫ] on peut écrire
Z x
(x − a) (x − a)n (n) (x − t)n (n+1)
f (x) = f (a) + f (a) + · · · + f (a) + f (t)dt.
1! n! a n!
Preuve. On utilise une démonstration par récurrence. La formule
Z x
f ′ (t)dt = f (x) − f (a),
a
assure que le théorème est vrai pour n = 0. Si on suppose vraie la formule à l’ordre
n ≥ 0 alors
Z x x Z x
(x − t)n (n+1) −(x − t)n+1 (n+1) (x − t)n+1 (n+2)
f (t)dt = f (t) + f (t)dt,
a n! (n + 1)! a a (n + 1)!
Z x Z x
(x − t)n (n+1) (x − a)n+1 (n+1) (x − t)n+1 (n+2)
f (t)dt = f (a) + f (t)dt,
a n! (n + 1)! a (n + 1)!
ce qui nous donne la formule à l’ordre n + 1.

2.3. La formule d’Euler-Maclaurin

2.3.1. Un exemple à la main. Soit f une fonction de classe C 3 sur [0, 1].
Cherchons à exprimer
Z 1
f (t)dt.
0
Pour cela on peut commencer par dire en remplaçant f par sa valeur en 0 qu’une
valeur approchée de l’intégrale est f (0). Étudions alors
Z 1
W = f (t)dt − f (0).
0
on voit que
Z 1
W = (f (t) − f (0)) dt
0
ce qui par intégration par partie (on dérive f(t) -f(0) et on intégre 1) donne
Z 1
1
W = [P1 (t) (f (t) − f (0))]0 − f ′ (t)P1 (t)dt,
0
10 2. LES OUTILS CLASSIQUES D’APPROXIMATION

où P1 (t) est une primitive de 1 (donc de la forme t − a) à bien choisir.

Z 1
W = (1 − a) (f (1) − f (0)) − f ′ (t)P1 (t)dt.
0
L’intégrale qui reste dans le second membre peut à son tour être intégrée par partie
en dérivant f ′ et en intégrant P1 . Soit P2 (t) une primitive de P1 (t). Choisissons
P1 (t) tel que P2 (1) = P2 (0) ou encore
Z 1
P1 (t)dt = 0.
0
Ceci impose de prendre a = 1/2 (P1 (t) = t − 1/2). On a alors
Z 1
W = 1/2 (f (1) − f (0)) − P2 (0) (f1′ (t) − f ′ (0)) + P2 (t)f (2) (t)dt.
0
Intégrons de nouveau par partie l’intégrale qui subsiste au second membre. Notons
P3 (t) une primitive de P2 (t) et choisissons P2 (t) de telle sorte que P3 (1) = P3 (0).
Comme P1 (t) = t − 1/2 on a P2 (t) = t2 /2 − t/2 + C et la condition imposée
Z 1
P2 (t)dt = 0
0
donne C = 1/12. Alors P3 (t) = t /6 − t2 /4 + t/12 + C, et si on impose aussi la
3

condition Z 1
P3 (t)dt = 0
0
alors P3 (t) = t3 /6 − t2 /4 + t/12. On obtient après une nouvelle intégration par
partie :
Z 1
W = 1/2 (f (1) − f (0)) − 1/12 (f1′ (t) − f ′ (0)) − P3 (t)f (3) (t)dt.
0
Arrêtons là le développement et écrivons en conclusion
Z 1 Z 1
f (t)dt = 1/2 (f (0) + f (1)) − 1/12 (f ′ (1) − f ′ (0)) − P3 (t)f (3) (t)dt.
0 0

2.3.2. Polynômes de Bernoulli. Les calculs précédents mettent en évidence

la suite des polynômes de Bernoulli.
Théorème 2.2. Il existe une suite (Qn )n≥0 et une seule de polynômes telle
que
(1) Q0 = 1
(2) Q′n = Qn−1 , pour n ≥ 1
R1
(3) 0 Qn (u)du = 0, pour n ≥ 1.
Ces polynômes sont appelés polynômes de Bernoulli.
Preuve. Par récurrence, Q0 est bien défini, Qn−1 étant construit, la condition (2)
fixe Qn à une constante près. La condition (3) fixe cette constante.

Compte tenu du mode de construction de ces polynômes, il est facile de voir

que leurs coefficients sont rationnels.
Voici les premiers polynômes de Bernoulli :
2.3. LA FORMULE D’EULER-MACLAURIN 11

1 x2 x 1
Q0 = 1 Q1 = x − − +
Q2 =
2 2 2 12
x3 x2 x
Q3 = − + .
6 4 12
Proposition 2.3. Pour n ≥ 2, Qn (1) = Qn (0).
R1
Preuve. En effet on doit avoir 0 Qn−1 (u)du = 0. Mais comme Qn est une primitive
R1
de Qn−1 on a 0 Qn−1 (u)du = Qn (1) − Qn (0), d’où le résultat.
xn−1
Proposition 2.4. Si n ≥ 1, Qn (x + 1) − Qn (x) = (n−1)! .

Preuve. On constate que l’égalité est vraie pour n = 1. Supposons la vraie pour
n, par primitivation on obtient alors l’égalité à l’ordre n + 1 à une constante près.
Mais la proposition précédente permet de conclure à la nullité de cette constante
d’intégration.
Pp
Cette égalité permet de calculer des sommes du type k=1 k n . Par exemple si
n = 2 on obtient
p
1X 2
Q3 (p + 1) − Q3 (1) = k ,
2
k=1
ce qui donne
p
X p(p + 1)(2p + 1)
k2 = .
6
k=1

Proposition 2.5. Pour tout n ≥ 0, Qn (1 − x) = (−1)n Qn (x).

Preuve. Le résultat est vrai pour n = 0 et n = 1. Supposons le résultat vrai pour
n ≥ 2. Alors par primitivation on obtient au rang n + 1 la formule voulue à une
constante d’intégration près
Qn+1 (1 − x) = (−1)n+1 Qn+1 (x) + C.
Si n + 1 est pair en donnant à x la valeur 0 on calcule C = 0.
Si n + 1 est impair alors par primitivation
Qn+2 (1 − x) = Qn+2 (x) + Cx + D,
en prenant x = 0 on obtient d’abord D = 0, puis en prenant x = 1 on obtient
C = 0.
Proposition 2.6. Pour n ≥ 1,
Q2n+1 (0) = Q2n+1 (1/2) = Q2n+1 (1) = 0.
Preuve. Il suffit d’appliquer la proposition précédente avec x = 0, ce qui donne
Q2n+1 (0) = Q2n+1 (1) = 0, puis avec x = 1/2, ce qui donne Q2n+1 (1/2) = 0.

En étudiant par récurrence les variations des fonctions Qn , on pourra montrer

que
Proposition 2.7. Pour n ≥ 1, Q2n (0) 6= 0 et Signe(Q2n (0)) = (−1)n+1 .
12 2. LES OUTILS CLASSIQUES D’APPROXIMATION

On notera
Bk = (−1)k+1 (2k)!Q2k (0).
Ainsi les Bk (nombres de Bernoulli) sont des rationnels positifs. Les premiers
nombres de Bernoulli sont
B1 = 1/6 B2 = 1/30 B3 = 1/42 B4 = 1/30.
2.3.3. Formule d’Euler-Maclaurin. Soit f une fonction réelle de classe C ∞
sur [0, 1].
Z 1
f (1) − f (0) = f ′ (t)dt,
0
Z 1
′
f (1) − f (0) = [Q1 (t)f (t)]10 − Q1 (t)f (2) (t)dt,
0
Z 1
′ ′
f (1) − f (0) = 1/2(f (1) + f (0)) − Q1 (t)f (2) (t)dt,
0
et par récurrence on montre que :
Théorème 2.8. Soit f une fonction réelle de classe C ∞ sur [0, 1]. Alors pour
tout n ≥ 1 :
n
1 X (−1)k Bk (2k)
f (1) − f (0) = f ′ (1) + f ′ (0) + f (1) − f (2k) (0)
2 (2k)!
k=1
Z 1
− Q2n+1 (x)f (2n+2) (x)dx.
0

En particulier si on applique ce résultat à une primitive de f on obtient

Z 1 n
1 X (−1)k Bk (2k−1)
f (t)dt = f (1) + f (0) + f (1) − f (2k−1) (0)
0 2 (2k)!
k=1
Z 1
− Q2n+1 (x)f (2n+1) (x)dx.
0
On peut appliquer le théorème 2.8 sur un intervalle [a, b] au lieu de [0, 1]. Il
suffit comme toujours de faire le changement de variable affine qui envoie a sur 0
et b sur 1 : t = a + u(b − a). On obtient alors :
Théorème 2.9. Soit f une fonction réelle de classe C ∞ sur [a, b]. Alors pour
tout n ≥ 1 :
n
b−a ′ X (−1)k Bk (b − a)2k (2k)
f (b) − f (a) = f (b) + f ′ (a) + f (b) − f (2k) (a)
2 (2k)!
k=1
Z b
u−a
− Q2n+1 ( )(b − a)2n−1 f (2n+2) (u)du.
a b−a
Découpons maintenant l’intervalle [a, b] en q morceaux de même longueur h =
(b − a)/q sous la forme :
a = a0 < a1 < a2 < · · · < aq−1 < aq = b,
appliquons le théorème 2.9 sur chaque morceau et sommons les résultats. Nous
obtenons :
2.3. LA FORMULE D’EULER-MACLAURIN 13

Théorème 2.10. Soit f une fonction réelle de classe C ∞ sur [a, b]. Soit a =
a0 < a1 < · · · < aq1 < aq = b le partage de [a, b] en q intervalles de même longueur
h = (b − a)/q. Alors pour tout n ≥ 1 :
q−1
X
h ′
f (b) − f (a) = f (aq ) − f (a0 ) = (f (b) + f ′ (a)) + h f ′ (as )
2 s=1
n
X (−1)k Bk h2k
+ f (2k) (b) − f (2k) (a)
(2k)!
k=1
Z q−1
!
h v X
−h2n−1 Q2n+1 f (2n+2) (as + v) dv.
0 h s=0

Cette dernière formule est appelée la formule sommatoire d’Euler-Maclaurin

car elle permet de calculer la somme :
q−1
X
f ′ (as ).
s=1

2.3.4. Application à l’évaluation de restes de séries. Considérons la

série de terme général 1/n2 . On sait que cette série converge et que :
∞
X 1 π2
S= = .
k2 6
k=1

Écrivons la somme S sous la forme :

S = Sp + Rp
où :
p
X ∞
X
1 1
Sp = et Rp = .
k2 k2
k=1 k=p+1

Nous cherchons à évaluer Rp . Pour cela introduisons la fonction f (x) = −1/x

dont la dérivée est f ′ (x) = 1/x2 . Appliquons le théorème 2.10 à la fonction f sur
l’intervalle [p, r + 1] (où p et r sont des entiers) avec h = 1 et n = 1. Nous obtenons
successivement :

1 1 1
Rp − Rr = + + ··· + 2,
(p + 1)2 (p + 2)2 r

1 1 1 1 1 1 1 1
Rp − Rr = − − 2
+ + − + Tp,r ,
p r+1 2 p (r + 1)2 6 p 3 (r + 1)3
et en faisant tendre r vers +∞ :
1 1 1
Rp = − + 3 + Tp ,
p 2p2 6p
où Tp se calcule facilement en utilisant le théorème 2.10. Un calcul simple (compa-
raison d’une somme de série avec une intégrale) permet de voir que :

1
Tp = O ,
p4
14 2. LES OUTILS CLASSIQUES D’APPROXIMATION

ce qui donne pour Rp le développement asymptotique :

1 1 1 1
Rp = − 2 + 3 + O .
p 2p 6p p4
CHAPITRE 3

Interpolation des fonctions

3.1. Introduction
L’interpolation est un sujet très vaste lié aux questions d’approximation
des fonctions. Très grossièrement il s’agit de trouver dans une classe fixée de fonc-
tions (par exemple les fonctions polynomiales) un élément réalisant un certain
nombre de contraintes. Souvent ces contraintes sont liées à la donnée d’une fonc-
tion f qu’on cherche à approcher par un procédé d’interpolation (par exemple la
fonction cherchée doit prendre la même valeur que f en des points donnés). On se
trouve alors confronté à plusieurs problèmes de natures différentes. Tout d’abord
un problème algébrique, celui de trouver le ou les éléments de la classe choisie
qui réalise les contraintes. Ensuite un problème d’approximation qui consiste lors-
qu’on est parti d’une fonction f à mesurer la qualité de l’approximation théorique
obtenue. Enfin un problème algorithmique, celui de déterminer un algorithme per-
formant qui permette de calculer facilement et de manière aussi exacte que possible
la ou les solutions.
Dans un premier temps nous partirons d’un exemple important : l’interpolation
de Lagrange.

3.2. Interpolation de Lagrange

Soient x0 , x1 , ..., xn des nombres complexes distincts et y0 , y1 , ..., yn des nombres
complexes. Il s’agit de trouver un polynôme P (X) vérifiant P (xk ) = yk pour toutes
les valeurs de k comprises entre 0 et n.

3.2.1. L’aspect algébrique.

Existence de solutions. Notons C[X] l’espace des polynômes à coefficients com-
plexes et Cn [X] le sous espace des polynômes à coefficients complexes de degré
inférieur ou égal à n. Considérons alors l’application linéaire T de C[X] dans Cn+1
qui à un polynôme P (X) fait correspondre (P (x0 ), P (x1 ), ..., P (xn )). On voit que le
noyau Ker(T ) de l’application T est l’espace constitué des multiples du polynôme
N (X) = (X − x0 )(X − x1 )...(X − xn ) et qu’on peut écrire
M
C[X] = Cn [X] Ker(T ).
Ceci nous montre que la restriction de T à Cn [X] est une bijection de Cn [X]
sur Cn+1 .
Le résultat obtenu est donc le suivant

Théorème 3.1. Pour tout élément (y0 , y1 , ..., yn ) de Cn+1 il existe un polynôme
P (X) de degré ≤ n et un seul (polynôme d’interpolation de Lagrange) tel que
15
16 3. INTERPOLATION DES FONCTIONS

P (xk ) = yk (0 ≤ k ≤ n). Tout polynôme Q(X) (de degré quelconque) vérifiant

aussi Q(xk ) = yk (0 ≤ k ≤ n) s’écrit sous la forme

Q(X) = P (X) + K(X)N (X).

Ecriture sous la forme naturelle (base des monômes). Si on cherche à trouver
explicitement le polynôme de Lagrange écrit sous la forme habituelle P (X) = a0 +
a1 X + ... + an X n on est amené à résoudre le système de n + 1 équations en les n + 1
inconnues a0 , ..., an


 a0 + a1 x0 + ... + an xn0 = y0

a0 + a1 x1 + ... + an xn1 = y1

 ············ ··· ···

a0 + a1 xn + ... + an xnn = yn .
Ce système est un système de Vandermonde dont on sait bien sûr qu’il admet
une solution unique. Nous fournirons ultérieurement un algorithme pour résoudre ce
système plus rapidement que ne le ferait une méthode classique comme la méthode
du pivot par exemple.
Ecriture sous forme de Lagrange. Ainsi le calcul des coefficients ai du polynôme
cherché se ramène à la résolution d’un système de Vandermonde. Sans être très
compliqué ceci n’est cependant pas immédiat. Or il peut se faire qu’on n’ait pas
absolument besoin des coefficients ai et qu’on puisse se satisfaire d’exprimer le
polynôme d’interpolation dans une base mieux adaptée que la classique base des
monômes. Dans cet ordre d’idée, introduisons les n + 1 polynômes Lk (X) (où 0 ≤
k ≤ n) qui vérifient

Lk (xk ) = 1
Lk (xj ) = 0 (j 6= k).
D’après l’étude qui précède le polynôme Lk (X) existe et est unique. On vérifie
aisément que Lk (X) s’écrit explicitement sous la forme

(X − x0 )...(X − xk−1 )(X − xk+1 )...(X − xn )

Lk (X) = .
(xk − x0 )...(xk − xk−1 )(xk − xk+1 )...(xk − xn )
Ces polynômes forment une base de Cn [X] et le polynôme d’interpolation s’ex-
prime dans cette base
n
X
P (X) = yk Lk (X).
k=0
Il est intéressant de remarquer que le polynôme Lk (X) s’écrit aussi

N (X)
Lk (X) =
(X − xk )N ′ (xk )
où N (X) = (X − x0 )(X − x1 )...(X − xn ).
Ecriture sous la forme de Newton. L’inconvénient des polynômes Lk (X) est que
chacun d’eux fait intervenir tous les points d’interpolation et donc si on rajoute un
nouveau point tout calcul fait à partir des polynômes Lk (X) doit être entièrement
refait.
3.2. Interpolation DE LAGRANGE 17

Prenons alors les polynômes Nk (X) = (X −x0 )(X −x1 )...(X −xk−1 ) où 0 ≤ k ≤
n (avec N0 = 1). Ces polynômes forment aussi une base de Cn [X] et le polynôme
d’interpolation se décompose sur cette base sous la forme de Newton
n
X
P (X) = bk Nk (X).
k=0
Le problème est alors de calculer les coefficients bk . Pour ce faire définissons les
différences divisées successives des valeurs yi par rapport aux points xi

[y0 ] = y0
[y1 , ..., yk ] − [y0 , ..., yk−1 ]
[y0 , y1 , ..., yk ] =
xk − x0
Théorème 3.2. Les coefficients de la décomposition du polynôme d’ interpo-
lation de Lagrange de degré n dans la base de Newton sont donnés par

bk = [y0 , y1 , ..., yk ]
où
0 ≤ k ≤ n.
Preuve. La formule à démontrer est clairement vraie si on a n = 0. Supposons
la formule vraie pour les polynômes de degré n−1 interpolant en n points. Soit alors
Pn−1 (X) le polynôme d’interpolation de Lagrange associé aux points x0 , x1 , ..., xn−1
et aux valeurs y0 , y1 , ..., yn−1 , Qn−1 (X) le polynôme d’interpolation de Lagrange
associé aux points x1 , x2 , ..., xn et aux valeurs y1 , y2 , ..., yn et
Xn
P (X) = bk Nk (X)
k=0
le polynôme d’interpolation de Lagrange associé aux points x0 , x1 , ..., xn et aux
valeurs y0 , y1 , ..., yn . Il est facile de voir que
n−1
X
Pn−1 (X) = bk Nk (X)
k=0
si bien qu’il reste simplement en vertu de l’hypothèse de récurrence à établir la
formule pour le coefficient bn .
Pour cela définissons
(X − x0 )Qn−1 (X) − (X − xn )Pn−1 (X)
P̃n (X) = .
xn − x0
On vérifie que
P̃n (xi ) = yi
pour 0 ≤ i ≤ n. Donc
P̃n (X) = P (X).
En égalant les coefficients du terme de degré n dans l’expression des polynômes Pn
et P̃n on obtient la relation cherchée.
En comparant l’expression de Newton du polynôme Pk d’interpolation en k + 1
points avec celle de Lagrange on trouve en regardant les coefficients des termes de
degré k
18 3. INTERPOLATION DES FONCTIONS

k
X yj
[y0 , y1 , ..., yk ] = ′ .
j=0
Nk+1 (xj )
Il est clair dans cette représentation que le rajout d’un nouveau point ne fait
que rajouter un nouveau terme au polynôme, les autres termes restant identiques.

3.2.2. L’aspect approximation.

Le théorème de division des fonctions différentiables. Soit f une fonction réelle
définie sur R de classe C p+1 , où p est un entier naturel. On suppose que f s’annule
en un point a de R. Posons
Z 1
g(x) = f ′ (a + (x − a)u)du
0
alors g(x) est l’unique fonction continue telle que

f (x) = (x − a)g(x).
De plus d’après le théorème de dérivation sous le signe intégrale on voit que la
fonction g(x) est de classe C p et que pour tout 0 ≤ q ≤ p
Z 1
(q)
g (x) = uq f (q+1) (a + (x − a)u)du
0
et par suite

1
|g (q) (x)| ≤ supt∈[a,x]|f (q+1) (t)|.
q+1
Majoration de l’erreur. Soit f une fonction de classe C n+1 , P le polynôme d’in-
terpolation de Lagrange qui prend les mêmes valeurs que f aux point x0 , x1 , ..., xn
et I un intervalle compact contenant x, x0 , x1 , ..., xn . Appliquons alors le théorème
précédent à f (x) − P (x). On obtient

f (x) − P (x) = (x − x0 )g0 (x)

avec
(n) 1
|g0 (x)| ≤ supt∈I |f (n+1) (t)|
n+1
(ne pas oublier que P (n+1) (x) = 0)
puis

g0 (x) = (x − x1 )g1 (x)

avec
(n−1) 1 (n)
|g1 (x)| ≤ supt∈I |g0 (t)|
n
et ainsi de suite. Si bien que

1
|f (x) − P (x)| ≤ |(x − x0 )(x − x1 )...(x − xn )|supt∈I |f (n+1) (t)|
(n + 1)!
3.2.3. L’aspect algorithmique.
3.2. Interpolation DE LAGRANGE 19

Le calcul des différences divisées. L’algorithme des différences divisées est très
simple. Il utilise l’écriture du polynôme d’interpolation sous la forme de Newton. Les
coefficients sont alors calculés par la formule de récurrence établie précédemment

[y0 ] = y0
[y1 , ..., yk ] − [y0 , ..., yk−1 ]
[y0 , y1 , ..., yk ] = .
xk − x0
si bien que le calcul se fait conformément à la figure 1.
Le nombre d’opérations à effectuer est en O(n2 ).

[y0 ]
❍
❍❍
❍❍
❥
❍
[y0 , y1 ]
❍
❍❍
✯
✟ ❍❍
✟✟ ❥
❍
✟✟
[y1 ] ✟ [y0 , y1 , y2 ]
❍❍
❍❍ ❍
❍❍
❍❍ ✟✟
✯ ❍❍
❥ ✟ ❥
✟✟
[y1 , y2 ] ✟ [y0 , y1 , y2 , y3 ]
❍❍
❍❍ ✯
✟
✟
✯ ❍❍
✟✟ ❥ ✟✟
✟✟ ✟✟
[y2 ] ✟ [y1 , y2 , y3 ]
❍❍
❍❍
❍❍
❥ ✟✟
✯
✟✟
✟
[y2 , y3 ] ✟

✟✟
✯
✟✟
[y3 ] ✟✟

Fig. 1. Le calcul des différences divisées

Résolution d’un système de Vandermonde. Nous avons vu que le calcul effec-

tif des coefficients du polynôme d’interpolation de Lagrange dans la base naturelle
des monômes passe par la résolution d’un système de Vandermonde. Voici un algo-
rithme qui permet de résoudre un tel système. Cet algorithme est basé en fait sur
l’algorithme de Hörner pour l’évaluation de polynômes. Il est plus rapide que les al-
gorithmes directs de résolution des systèmes linéaires généraux comme la méthode
du pivot de Gauss ou la méthode de Householder qui sont en O(n3 ) alors que nous
obtenons ici un algorithme en O(n2 ).
20 3. INTERPOLATION DES FONCTIONS

Soit N un entier ≥ 2. Etant donné x = (x1 , x2 , ..., xN ) un N-uplet de réels deux

à deux distincts on note B la matrice
 
1 1 ··· 1
 x1 x2 ··· xN 
 
 x21 x 2
· · · x 2 
B= 2 N 
 .. .. .. 
 . . . 
N −1 N −1 N −1
x1 x2 · · · xN
On cherche à résoudre le système
BW = Q
où Q est la matrice colonne constituée des seconds membres q1 , q2 , · · · , qN du
système et où W est la matrice colonne constituée des inconnues w1 , w2 , · · · , wN
du système.

Pour tout entier j vérifiant 1 ≤ j ≤ N on pose

YN
x − xn
Pj (x) =
n=1
xj − xn
n6=j

Pj est donc un polynôme en x de degré N − 1 qui peut s’écrire

N
X
Pj (x) = Aj,k xk−1
k=1

En effectuant le produit de la matrice A = (Aj,k )j,k par la matrice B on constate

que
AB = (Pj (xk ))j,k
ce qui prouve que A est l’inverse de B.
On peut alors écrire que W=AQ. On obtient ainsi les formules
N
X
wj = Aj,k qk .
k=1

Nous allons dans la suite mettre en place une méthode de résolution qui calcule
les coefficients des polynômes Pj , donc qui calcule l’inverse de la matrice B.Pour
calculer les coefficients de Pj on sera amené à calculer les coefficients de
N
Y
Nj (x) = (x − xn )
n=1
n6=j

et aussi le dénominateur intervenant dans la formule qui définit Pj , c’est à dire le

nombre Nj (xj ).
Posons
P (x) = (x − x1 )(x − x2 )...(x − xN )
P (x) est donc un polynôme de degré N qui s’écrit sous la forme :
P (x) = xN + cN xN −1 + ... + c2 x + c1
3.2. Interpolation DE LAGRANGE 21

Montrons tout d’abord comment si on connaı̂t les coefficients cj on peut calculer

les coefficients du polynôme
N
Y
Nj (x) = (x − xn )
n=1
n6=j

Pour cela posons

Nj (x) = bN xN −1 + ... + b2 x + b1
On vérifie immédiatement sur l’expression de Nj (x) que le coefficient du terme
de plus haut degré est 1. En remarquant que P (x) = Nj (x)(x − xj ) on établit la
formule
bk−1 = ck + xj bk .
Si bien que

bN = 1
bk−1 = ck + xj bk
Connaissant les coefficients de Nj il est alors facile de calculer le dénominateur
Nj (xj ) intervenant dans la définition de Pj .
En effet posons tN = bN = 1 et définissons pour tout k ≤ N
tk−1 = xj tk + bk−1
On constate alors que t1 = Nj (xj ), le calcul proposé pour Nj (xj ) n’étant rien
d’ autre que l’algorithme de Horner.
Il reste maintenant à calculer les coefficients cj de P .
Pour tout entier k vérifiant 1 ≤ k ≤ N on définit
Qk (x) = (x − x1 )(x − x2 )...(x − xk )
et on écrit Qk sous la forme
Qk (x) = xk + αk,k xk−1 + αk,k−1 xk−2 + ... + αk,1
Il est facile de voir sur l’expression de Q1 (x) = x − x1 que α1,1 = −x1 .
De la formule
Qk (x) = Qk−1 (x)(x − xk ).
découlent pour k = 2, 3, ..., N les formules
αk,k = αk−1,k−1 − xk

αk,j = αk−1,j−1 − xk αk−1,j j = k − 1, ..., 2

ce qui achève l’algorithme.
On peut voir que le nombre d’opérations à faire dans cet algorithme est en
O(N 2 ), la partie la plus coûteuse étant le calcul des coefficients ck .

3.2.4. Un cas particulier : les points d’interpolation sont les racines

ne de l’unité.
22 3. INTERPOLATION DES FONCTIONS

Interpolation de Lagrange et transformée de Fourier discrète. Rappelons que

si a = (a0 , a1 , ..., an−1 ) est une suite finie de n nombres complexes on définit la
transformée de Fourier b a de la suite a comme étant la suite b
a = (b
a0 , b
a1 , ..., b
an−1 ) où
n−1
1X 2iπuv
av =
b au e − n .
n u=0
La transformation inverse se calcule facilement par
n−1
X 2iπuv
au = av e
b n .
v=0

On remarque que
b 1
au = an−u .
b
n
A un coefficient près le même algorithme permettra de calculer la transformée de
Fourier et son inverse.
Notons
Pba (X) = ba0 + b an−1 X n−1
a1 X + ... + b
On voit alors que
2iπu
au = Pba (e n ).
Cette dernière remarque met en évidence un aspect très important de la transformée
de Fourier discrète : l’aspect interpolation . En effet il est facile de trouver grâce à ce
que nous avons vu le polynôme d’interpolation de Lagrange qui prend les valeurs au
2iπu
aux points xu = e n ; C’est le polynôme Pba (X) dont les coefficients sont donnés
par la transformée de Fourier discrète de la suite a = (a0 , a1 , ..., an−1 ) des valeurs
prises aux points d’interpolation .
Le calcul explicite : transformée de Fourier rapide. Il existe divers façons
proches les une des autres de calculer une transformée de Fourier discrète. Toutes
ces variantes sont des algorithmes de transformée de Fourier rapides (FFT). Nous
nous placerons ici dans le cas où le nombre d’éléments de la suite à transformer est
n = 2m . Pour tout r > 0 et tout 0 ≤ k ≤ 2r − 1 posons
2ikπ
W2kr = e− 2r .
Remarquons que
r
W2kr = (W2kr+1 )2 = (W2k+2
r+1 )
2

r
W2kr+1 = −W2k+2
r+1

par exemple
(W83 )2 = (W87 )2 = W43
W83 = −W87 .
On rappelle que si
a = (a0 , a1 , ..., a2m −1 )
et si
1 2m −1

Pa (X) = a 0 + a 1 X + ... + a 2m −1 X
2m
alors
au = Pa (W2um ).
b
3.2. Interpolation DE LAGRANGE 23

Pour tout polynôme

r
P (X) = p0 + p1 X + ... + p2r −1 X 2 −1

notons
r−1
P0 (X) = p0 + p2 X + ... + p2r −2 X 2 −1

et
r−1
P1 (X) = p1 + p3 X + ... + p2r −1 X 2 −1

alors
P (X) = P0 (X 2 ) + XP1 (X 2 )
ce qui donne si 0 ≤ k ≤ 2r−1 − 1
P (W2kr ) = P0 (W2kr−1 ) + W2kr P1 (W2kr−1 )
et
r−1
P (W2k+2
r ) = P0 (W2kr−1 ) − W2kr P1 (W2kr−1 ).
Ces dernières formules vont nous donner un algorithme pour calculer les valeurs
de la transformée de Fourier.
Remarquons tout d’abord que si on a tabulé les valeurs de W2km alors on dispose
aussi des valeurs de W2kr pour tout r ≤ m.

W84 W85 W86 W87

W80 W81 W82 W83
−W80 −W81 −W82 −W83
W42 W43
W40 W41
−W40 −W41
W21
W20
−W20

Pratique du calcul. Le coefficient n1 n’interviendra qu’à la fin. Pour cela au

lieu de calculer avec le polynôme Pa nous calculerons avec P = nPa = a0 + ... +
m
a2m −1 X 2 −1 .
L’exemple m = 3 est suffisamment instructif pour décrire l’algorithme. Remar-
quons que

P000 (X) = a0 , P001 (X) = a4 , P010 (X) = a2 , P011 (X) = a6

P100 (X) = a1 , P101 (X) = a5 , P110 (X) = a3 , P111 (X) = a7 .

On commence donc à faire une permutation σ des éléments

a 0 , a1 , a2 , a3 , a4 , a5 , a6 , a7
pour les mettre dans l’ordre
a 0 , a4 , a2 , a6 , a5 , a3 , a7 .
Ceci se fait facilement en remarquant qu’à chaque indice supposé écrit en binaire
on fait correspondre l’indice obtenu en écrivant les bits dans l’ordre inverse. Ainsi
l’indice 4 = 100 est transformé en 1 = 001. la suite du calcul de la transformée
de Fourier se fait en trois étapes indiquées par la figure 2 et à la fin on divise les
coefficient obtenus par 8.
24 3. INTERPOLATION DES FONCTIONS
a0 a4 a2 a6 a1 a5 a3 a7
❝ ❝ ❝ ❝ ❝ ❝ ❝ ❝
❆ ✁ ❆ ✁ ❆ ✁ ❆ ✁
❆ ✁ ❆ ✁ ❆ ✁ ❆ ✁
❆ ✁ ❆ ✁ ❆ ✁ ❆ ✁
❆✁ ❆✁ ❆✁ ❆✁
✁❆ ✁❆ ✁❆ ✁❆
✁
W0
❆ −W20
✁
W0
❆ −W20
✁ ❆
W0 −W20
✁ ❆
W0 −W20
✁ 2 ❆ ✁ 2 ❆ ✁ 2 ❆ ✁ 2 ❆
☛
✁ ❯
❆ ✁☛ ❆❯ ✁☛ ❆❯ ✁☛ ❆❯
❝❄ ❄
❝ ❄
❝ ❝❄ ❄
❝ ❝❄ ❄
❝ ❄
❝
❅ ❅ ❅ ❅
❅ ❅ ❅ ❅
❅ ❅ ❅ ❅
❅ ❅ ❅ ❅
❅ ❅ ❅ ❅
❅ ❅ ❅ ❅
W40 ❅ W41 −W40❅ −W41 W40 ❅ W41 −W40❅ −W41
❅ ❅ ❅ ❅
❄✠ ❄✠ ❘ ❄
❅ ❘ ❄
❅ ❄ ✠ ❄ ✠ ❘ ❄
❅ ❘
❅ ❄
❝ ❝ ❝ ❝ ❝ ❝ ❝ ❝
❍❍ ❍❍ ❍❍ ❍❍✟✟ ✟✟ ✟✟ ✟✟
❍❍ ❍❍ ❍❍ ✟✟❍❍ ✟✟ ✟ ✟ ✟ ✟
❍❍ ❍❍ ✟❍ ✟❍ ✟✟ ❍❍ ✟✟ 1 ✟✟ 2
❍❍ ❍✟
✟❍ ❍
✟❍✟ −W80❍✟
✟❍ −W8 ✟
✟ −W8 −W83
❍❍ ✟✟ ❍❍ ✟✟ ❍✟ ❍✟ ❍✟ ❍✟
✟
✟W80 ❍ ✟
❍✟W81 ❍ ✟
❍✟W82 ❍ ✟
❍✟W83 ❍ ❍
✟✟ ✟✟❍❍ ✟✟❍❍ ✟✟❍❍ ❍❍
✟✟ ✟✟ ✟❍ ✟❍ ✟✟ ❍❍ ❍❍ ❍❍
❄ ✟ ✟ ❄ ✟ ✟ ❄ ✟ ✟ ❄ ✟❍❍ ✟ ❄ ❍ ❍ ❄ ❍ ❍ ❄ ❍❍ ❄
❝✟✙ ❝✟✙ ❝✟✙ ❝✟✙ ❥ ❝
❍ ❥ ❝
❍ ❥ ❝
❍ ❥
❍ ❝
ab0 ab1 ab2 ab3 ab4 ab5 ab6 ab7

Fig. 2. La FFT sur 8 points

Appelons M81 , M82 , M83 les matrices

 
1 1 0 0 0 0 0 0
 1 −1 0 0 0 0 0 0 
 
 0 0 1 1 0 0 0 0 
 
1
 0 0 1 −1 0 0 0 0 
M8 =   0


 0 0 0 1 1 0 0 
 0 0 0 0 1 −1 0 0 
 
 0 0 0 0 0 0 1 1 
0 0 0 0 0 0 1 −1
 
1 0 1 0 0 0 0 0
 0 1 0 1 0 0 0 0 
 
 1 0 −1 0 0 0 0 0 
 
 0 1 0 −1 0 0 0 0 
M82 =   0


 0 0 0 1 0 1 0 
 0 0 0 0 0 1 0 1 
 
 0 0 0 0 1 0 −1 0 
0 0 0 0 0 1 0 −1
3.4. QUELQUES EXEMPLES IMPORTANTS 25

 
1 0 0 0 1 0 0 0
 0 1 0 0 0 1 0 0 
 
 0 0 1 0 0 0 1 0 
 
 0 0 0 1 0 0 0 1 
M8 = 
3
 1


 0 0 0 −1 0 0 0 
 0 1 0 0 0 −1 0 0 
 
 0 0 1 0 0 0 −1 0 
0 0 0 1 0 0 0 −1
S1 , S2 , S3 les matrices diagonales définies par

S1 = Diag(1, W20 , 1, W20 , 1, W20 , 1, W20 )

S2 = Diag(1, 1, W40, W41 , 1, 1, W40 , W41 )
S3 = Diag(1, 1, 1, 1, W80, W81 , W82 , W83 )
et enfin Σ la matrice de la permutation ”reverse bit” σ.
Dans ces conditions la matrice F8 de la transformation de Fourier sur 8 points
s’écrit
1
F8 = M83 S3 M82 S2 M81 S1 Σ.
8
Ceci se généralise facilement pour n = 2m . Le nombre d’opérations à effectuer
pour calculer cette transformation est de l’ordre de nlog(n).

3.3. Le problème général de l’interpolation

Interprétons de manière un peu plus algébrique le problème d’interpolation de
Lagrange. En particulier notons gi la forme linéaire sur Cn [X] qui à tout polynôme
Q(X) fait correspondre Q(xi ). On cherche alors un élément P (X) de Cn [X] qui
réalise gi (P ) = yi pour tout i. Remarquons que les gi forment une base du dual de
Cn [X] et que cette base est la base duale de la base constituée par les Li .
Nous poserons le problème général de l’interpolation en ces termes :
Soit E un espace vectoriel de dimension n, E ∗ le dual de E, (gi )i une base de
E ∗ et (yi )i des nombres. Trouver un élément P de E tel que gi (P ) = yi pour tout
1 ≤ i ≤ n.
Il est clair que si (ei )i est la base de E dont (gi )i est la base duale, alors
n
X
P = yi ei .
i=1
Nous allons voir par la suite quelques exemples qui entrent dans ce cadre
général.

3.4. Quelques exemples importants

3.4.1. Interpolation d’Hermite. Soient x0 < x1 et y0 , y1 , y0′ , y1′ des
nombres réels. On cherche un polynôme P de degré ≤ 3 tel que

 P (x0 ) = y0


P (x1 ) = y1
 P ′ (x0 ) = y0′

 ′
P (x1 ) = y1′
26 3. INTERPOLATION DES FONCTIONS

Notons E l’espace vectoriel de dimension 4 des polynômes de degré ≤ 3. Définissons

les 4 formes linéaires sur E

 δ0,x0 (P ) = P (x0 )


δ0,x1 (P ) = P (x1 )

 δ1,x0 (P ) = P ′ (x0 )

δ1,x1 (P ) = P ′ (x1 )
La question posée est donc de résoudre le problème linéaire suivant : trouver P tel
que


 δ0,x0 (P ) = y0

δ0,x1 (P ) = y1

 δ1,x0 (P ) = y0′

δ1,x1 (P ) = y1′
Pour cela cherchons si on peut trouver 4 polynômes H0,x0 , H0,x1 , H1,x0 , H1,x1 tels
que

1 si i = k et j = l
δi,xj (Hk,xl ) = .
0 sinon
Si on arrive à trouver ces polynômes cela prouvera à la fois que ces polynômes
forment une base de E, que les formes linéaires introduites forment une base de E ∗ ,
qui est la base duale de la base polynômiale trouvée.
Un simple calcul nous permet de trouver effectivement ces polynômes (ce sont
en fait des solutions dans des cas particuliers bien choisis du problème posé) ;

(x − x1 )2 (3x0 − x1 − 2x)
H0,x0 (x) =
(x0 − x1 )3

(x − x0 )2 (3x1 − x0 − 2x)
H0,x1 (x) =
(x1 − x0 )3

(x − x1 )2 (x − x0 )
H1,x0 (x) =
(x0 − x1 )2

(x − x0 )2 (x − x1 )
H1,x1 (x) = .
(x0 − x1 )2
Il ressort de toutes ces considérations que le problème a une solution unique
(polynôme d’interpolation de Hermite) donnée par

P (x) = y0 H0,x0 (x) + y1 H0,x1 (x) + y0′ H1,x0 (x) + y1′ H1,x1 (x).
Si nous sommes partis d’une fonction f de classe C 4 sur un intervalle com-
pact I contenant les points x0 , x1 , et si nous appelons P le polynôme d’interpola-
tion de Hermite associé aux points x0 , x1 et aux valeurs f (x0 ), f (x1 ), f ′ (x0 ), f ′ (x1 ),
alors comme dans l’exemple de l’interpolation de Lagrange, l’application répétée
du théorème de division des fonctions différentiables nous donne l’approximation

1
|f (x) − P (x)| ≤ (x − x0 )2 (x − x1 )2 sup |f (4) (t)|.
4! t∈I
3.4. QUELQUES EXEMPLES IMPORTANTS 27

3.4.2. Interpolation par les splines cubiques. Soit ∆ un partage d’un

segment [a, b],
a = x1 < x2 < · · · < xn = b.
Définition 3.3. Une fonction f définie sur [a, b] est une fonction spline cubique
relativement au partage ∆ si les conditions suivantes sont satisfaites :
1) f est de classe C 2 [a, b] ;
2) f coı̈ncide avec un polynôme de degré 3 sur chaque intervalle [xj , xj+1 ] ;
Nous noterons S∆ l’ensemble de ces fonctions.
Ainsi, S∆ est un sous-espace vectoriel de l’espace des fonctions définies sur [a, b].
Soit Y = (y1 , y2 , · · · , yn ) une suite de n nombres réels. Nous noterons S∆,Y
l’ensemble des fonctions splines cubiques qui vérifient la condition :
3) f (xj ) = yj pour 1 ≤ j ≤ n.

Remarque 3.4. Soit Y0 le n-uple particulier :

Y0 = (0, 0, · · · , 0).
On voit tout de suite que S∆,Y0 est un sous-espace vectoriel de S∆ .
Nous nous proposons de déterminer des fonctions splines cubiques répondant à
certaines conditions supplémentaires.
Soit donc f ∈ S∆,Y et posons Mj = f ′′ (xj ). Sur l’intervalle [xj , xj+1 ] la dérivée
seconde de f est de degré 1 et on a sur cet intervalle
xj+1 − x x − xj
f ′′ (x) = Mj + Mj+1 .
xj+1 − xj xj+1 − xj
En intégrant deux fois,
Mj (xj+1 − x)3 Mj+1 (x − xj )3
f (x) = + + Q(x)
6 (xj+1 − xj ) 6 (xj+1 − xj )
où Q(x) est un polynôme de degré 1 que l’on peut présenter sous la forme
Q(x) = α(xj+1 − x) + β(x − xj ).
Evaluons α et β en utilisant f (xj ) = yj et f (xj+1 ) = yj+1 . On trouve
!
(xj+1 − xj )2 1
α = yj − M j
6 (xj+1 − xj )
!
(xj+1 − xj )2 1
β = yj+1 − Mj+1 ,
6 (xj+1 − xj )
si bien que si on pose
hj+1 = xj+1 − xj
on obtient pour x ∈ [xj , xj+1 ],
Mj (xj+1 − x)3 Mj+1 (x − xj )3
f (x) = + +
6 hj+1 6 hj+1
! !
h2j+1 (xj+1 − x) h2j+1 (x − xj)
yj − M j + yj+1 − Mj+1 .
6 hj+1 6 hj+1
28 3. INTERPOLATION DES FONCTIONS

Nous allons exploiter maintenant les conditions sur la dérivée première. Pour 2 ≤
j ≤ n − 1 on a
hj+1 hj+1 yj+1 − yj
f ′ (x+
j ) = −Mj − Mj+1 +
3 6 hj+1
hj hj yj − yj−1
f ′ (x−
j ) = Mj−1 + Mj +
6 3 hj
si bien que
hj hj + hj+1 hj+1 yj+1 − yj yj − yj−1
Mj−1 + Mj + Mj+1 = − .
6 3 6 hj+1 hj
On dispose donc de n − 2 équations linéaires pour déterminer les n inconnues
M1 , M2 , · · · , Mn . Il convient donc si on espère avoir une solution unique de donner
deux conditions supplémentaires. Pour la suite du calcul nous imposerons donc les
valeurs de la dérivée de f en x1 et en xn , c’est-à-dire
f ′ (x1 ) = y1′
f ′ (xn ) = yn′ .
On a alors !
6 y2 − y1
2M1 + M2 = − y1′
h2 h2
et !
6 ′ yn − yn−1
Mn−1 + 2Mn = y − .
hn n hn
Si on pose

 λ1 = 1

 λ
n = 0
hj+1

 λj = hj +hj+1 2≤j ≤n−1

µj = 1 − λj 1≤j≤n
et  !

 y2 −y1

 b1 = 6
− y1′

 h2 h2

 !


6 yn −yn−1
bn = hn yn′ − hn



 !



 6 yj+1 −yj yj −yj−1
 bj
 = hj +hj+1 hj+1 − hj 2≤j ≤n−1

alors le système à résoudre est

    
2 λ1 0 0 ··· ··· 0 M1 b1
 µ2 2 λ2 0 ··· ··· 0  M2   b2 
    
 ..  ..   .. 
 0 µ3 2 λ3 .  .   . 
    
 .. ..  ..   .. 
 0 0 µ4 2 . .  . = . .
    
 . . .. ..  ..   .. 
 .. .. . . λn−2 0  .   . 
    
 . .  ..   .. 
 .. .. µn−1 2 λn−1  .   . 
0 0 ··· ··· 0 µn 2 Mn bn
3.4. QUELQUES EXEMPLES IMPORTANTS 29

Ce système est tridiagonal, à diagonale strictement dominante ; il possède une so-

lution et une seule. On peut employer pour un tel système une simplification de la
méthode du pivot, qui dans ce cas particulier s’exécute en temps linéaire. Ainsi on
peut énoncer le théorème suivant :
Théorème 3.5. Soient a′ et b′ deux nombres réels. Il existe une unique spline
cubique f ∈ S∆,Y telle que f ′ (a) = a′ et f ′ (b) = b′ .
Remarque 3.6. En particulier la seule fonction spline cubique qui s’annule sur
tous les nœuds du partage ∆ et dont la dérivée est nulle en a et en b est la fonction
nulle.

Les fonctions splines cubiques minimisent approximativement l’énergie de

flexion d’une tige. En effet si la fonction g(x) représente l’équation d’une tige par-
faitement élastique son énergie de flexion est
2
Z b g ′′ (t) dt
2 5/2
a
1 + g (t)′

et donc si la dérivée de g(x) est petite on peut approcher cette énergie de flexion
par
Z b 2
g ′′ (t) dt.
a
Nous allons voir que les fonctions splines cubiques minimisent cette dernière
expression. Sur C 2 [a, b] on définit un semi-produit scalaire par
Z b
< f, g >= f ′′ (t)g ′′ (t)dt
a
qui donne la semi-norme
Z !1/2
b
|f | = |f ′′ (t)|2 dt .
a

Soit N∆ le sous-espace de C 2 [a, b] constitué des fonctions nulles aux points xj

et telles que f ′ (a) = f ′ (b) = 0.
Lemme 3.7. Les espaces N∆ et S∆ sont orthogonaux et d’intersection réduite
à {0}.
Preuve. Soit f ∈ N∆ et g ∈ S∆ . Calculons :
Z b
< f, g >= f ′′ (t)g ′′ (t)dt.
a
Pour cela calculons pour chaque intervalle [xj , xj+1 ] du partage :
Z xj+1
Ij = f ′′ (t)g ′′ (t)dt.
xj

Par intégration par partie on obtient :

Z xj+1
′ ′′ x
Ij = [f (t)g (t)]xj+1
j
− f ′ (t)g (3) (t)dt,
xj
30 3. INTERPOLATION DES FONCTIONS

ce qui donne en intégrant de nouveau par partie le reste et en tenant compte du

fait que la dérivée quatrième de g est nulle et que f vaut 0 aux points xj et xj+1 :
x
Ij = [f ′ (t)g ′′ (t)]xj+1
j
.
On en conclut que :
Z b n−1
X b
f ′′ (t)g ′′ (t)dt = Ij = [f ′ (t)g ′′ (t)]a ,
a j=1
′
et puisque f (t) s’annule en a et b on en conclut que < f, g >= 0. Le fait que
l’intersection est réduite à la fonction nulle est une conséquence de la remarque 3.6.

Lemme 3.8. Toute fonction f ∈ C 2 [a, b] s’écrit d’une façon unique sous la
forme f = s + h où s ∈ S∆ et h ∈ N∆ .
Preuve. Posons Y = (f (x1 ), f (x2 ), · · · , f (xn )). Comme h doit appartenir à N∆ ,
s’il existe une spline cubique répondant à la question alors nécessairement s ∈
S∆,Y , s′ (a) = f ′ (a) et s′ (b) = f ′ (b). On sait qu’il existe une unique spline cubique
répondant à la question. La fonction h est alors uniquement définie par h = f − s,
et elle appartient bien à N∆ .
Théorème 3.9. Parmi toutes les fonctions f ∈ C 2 [a, b] qui vérifient f (xj ) = yj
et les conditions f ′ (a) = y1′ et f ′ (b) = yn′ la fonction qui minimise l’énergie de
flexion est la spline cubique (l’élément de S∆,Y ) qui vérifie en outre f ′ (a) = y1′ et
f ′ (b) = yn′ .
Preuve. Soit f ∈ C 2 [a, b]. Écrivons d’après le lemme précédent f = s+h, où s ∈ S∆
et h ∈ N∆ . On a en vertu de l’orthogonalité de s et h :
< f, f >=< s, s > + < h, h > .
Donc :
< s, s >=< f, f > − < h, h >≤< f, f > .

On pourrait aussi montrer le résultat suivant
Théorème 3.10. Parmi toutes les fonctions f ∈ C 2 [a, b] qui vérifient f (xj ) =
yj la fonction qui minimise l’énergie de flexion est la spline cubique (l’élément de
S∆,Y ) qui vérifie en outre f ′′ (a) = f ′′ (b) = 0.
3.4.3. Interpolation par des polynômes trigonométriques. Etant
donnés 2n + 1 nombres réels distincts
−π ≤ x1 < x2 < · · · < x2n+1 < π
et 2n+1 nombres complexes y1 , · · · , y2n+1 on cherche un polynôme trigonométrique
X n
P (x) = ap cos(px) + bp sin(px)
p=0

tel que P (xi ) = yi .

On pourrait étudier directement ce problème, mais ici nous allons plutôt le
relier au problème de l’interpolation de Lagrange que nous avons déjà étudié.
3.4. QUELQUES EXEMPLES IMPORTANTS 31

Pour cela posons pour tout −n ≤ k ≤ n

1
ck = a|k| − iSigne(k)b|k|
2
alors n
X
P (x) = ck eikx
k=−n
et aussi
2n
X
einx P (x) = ck−n eikx .
k=0
En posant X = eix on obtient
2n
X
einx P (x) = ck−n X k ,
k=0
ce qui ramène notre problème à un problème d’interpolation de Lagrange. Remar-
quons que dès que les ck sont connus, on calcule facilement les ap et les bp grâce
aux formules
ap = cp + c−p
bp = i(cp − c−p ).
CHAPITRE 4

Calcul numérique des intégrales définies

4.1. Introduction
Nous voulons calculer numériquement les intégrales définies
Z b
f (x)dx
a
où f est une fonction suffisamment régulière pour assurer l’existence des dérivées
et des bornes supérieures dont nous aurons besoin.
Les méthodes que nous décrivons dans un premier temps s’appuient sur la
démarche suivante :
– On découpe l’intervalle d’intégration en N morceaux de même longueur
a = x0 < x1 < · · · < xN = b
b−a
où xi+1 − xi = N .
– Sur chaque morceau on calcule la valeur approchée de l’intégrale
Z xi+1
f (x)dx
xi

en remplaçant sur l’intervalle [xi , xi+1 ] la fonction f par une fonction po-
lynômiale qui interpole f .
Les diverses méthodes qui entrent dans ce cadre diffèrent par le degré des po-
lynômes d’interpolation utilisés et par le choix des points des segments [xi , xi+1 ] en
lesquels on interpole f . En particulier pour ce dernier point, nous verrons comment
choisir au mieux les points d’interpolation. Il faut éviter de confondre les points
xi du partage en N morceaux du segment initial, avec les points des partages des
intervalles [xi , xi+1 ], que nous serons amenés à introduire pour appliquer sur ces
intervalles des méthodes interpolatoires.
Comme nous le verrons ces méthodes peuvent être éventuellement complétées
par une méthode d’accélération de convergence. Enfin nous donnerons un abord
purement analytique d’une classe de formules de calcul approchée d’intégrales.

4.2. Mise en œuvre de méthodes interpolatoires

Dans toute cette section pour la clarté des calculs et des méthodes, nous tra-
vaillerons sur l’intervalle [−1, 1] au lieu de [xi , xi+1 ], puis donnerons les formules
qu’on obtient de la même façon sur les intervalles [xi , xi+1 ], et enfin nous som-
merons ces diverses formules pour obtenir l’approximation de l’intégrale sur [a, b].
Dans chaque cas nous noterons e l’erreur de la méthode employée sur [−1, 1], ei
l’erreur de la méthode employée sur [xi , xi+1 ] et enfin E l’erreur globale introduite
sur le segment [a, b].
33
34 4. CALCUL NUMÉRIQUE DES INTÉGRALES DÉFINIES

Signalons enfin que nous supposerons la fonction f à intégrer sur le seg-

ment [a, b] continument dérivable autant de fois qu’il le faut pour pouvoir
appliquer les théorèmes permettant l’établissement des majorations des erreurs. En
particulier on aura à appliquer souvant les théorèmes de majoration de la différence
d’une fonction et d’un de ses polynômes interpolateurs.
Nous noterons
m(p) (f ) = sup f (p) (t) ,
t∈[−1,1]
(p)
mi (f ) = sup f (p) (t) ,
t∈[xi ,xi+1 ]

M (p) (f ) = sup f (p) (t) .

t∈[a,b]

4.2.1. Méthode des rectangles. La méthode que nous présentons mérite à

peine le nom de méthode numérique, car elle ne donne pas un calcul efficace et n’est
pas employée à cet usage. Cependant elle est importante pour la clarté de l’exposé,
car elle initialise en quelque sorte un processus qui va nous conduire à des méthodes
plus efficaces. Grâce à elle nous pourrons présenter et illustrer un certain nombre
de problèmes de fond.

1111111111
0000000000
0000000000
1111111111
0000000000
1111111111
00000000001
1111111111
−1

Fig. 1. Méthode des rectangles à gauche

La méthode des rectangles, inspirée de la définition même de l’intégrale de

Riemann consiste à supposer la fonction constante sur tout l’intervalle [−1, 1], la
valeur de la constante étant celle prise en un point déterminé par la fonction.
Nous prendrons ici la valeur de la fonction au point −1 (méthode des rectangles à
gauche). Remarquons qu’il s’agit bien d’une méthode d’interpolation de Lagrange
en un point par un polynôme de degré zéro. Ainsi
Z 1 Z 1
f (t)dt = f (−1)dt + e = 2f (−1) + e.
−1 −1
En écrivant Z u
f (u) = f (−1) + f ′ (t)dt,
−1
4.2. MISE EN ŒUVRE DE MÉTHODES INTERPOLATOIRES 35

on obtient Z Z
1 1
(1)
f (u)du − 2f (−1) ≤ m (f ) (u + 1)du,
−1 −1
ou encore
|e| ≤ 2m(1) (f ).
Un calcul analogue sur l’intervalle [xi , xi+1 ] donne
Z xi+1
f (t)dt = (xi+1 − xi )f (xi ) + ei ,
xi
avec
(xi+1 − xi )2 (1)
|ei | ≤ mi (f ).
2
Enfin sur l’intervalle [a, b] nous obtenons la formule globale
Z b
(b − a)
f (u)du = [f (x0 ) + · · · + f (xN −1 )] + E,
a N
avec
(b − a)2 (1)
|E| ≤ M (f ).
N

4.2.2. Méthode du milieu. La méthode des rectangles utilise la valeur de

la fonction f en un point déterminé de l’intervalle. Il faudrait voir si un bon choix
de ce point ne conduirait pas à une optimisation de la ”qualité de la méthode”.
Encore faudrait il donner un sens à cette expression, ce que nous ferons par la
suite. En attendant exposons la méthode du milieu qui consiste à utiliser comme
valeur approchée de la fonction sur l’intervalle [−1, 1], sa valeur au point milieu
(c’est-à-dire en 0).

1111111111
0000000000
0000000000
1111111111
0000000000
1111111111
0000000000
1111111111
0000000000
1111111111
111111111111111
000000000000000
00000000001
1111111111
−1

Fig. 2. Méthode du milieu

Ainsi Z Z
1 1
f (t)dt = f (0)dt + e = 2f (0) + e
−1 −1
En écrivant Z u
f (u) = f (0) + f ′ (t)dt,
0
36 4. CALCUL NUMÉRIQUE DES INTÉGRALES DÉFINIES

on obtient
Z 1 Z 1 Z u
′
f (u)du = 2f (0) + f (t)dt du.
−1 −1 0

Or la méthode d’intégration par partie appliquée aux fonctions 1 et f ′ (t) en dérivant

f ′ (t) et en intégrant 1 (on choisit t−u comme primitive de la fonction 1) nous donne
Z u Z u
f ′ (t)dt = uf ′ (0) − (t − u)f (2) (t)dt.
0 0

Donc
Z 1 Z u Z 1 Z u
′ (2)
f (t)dt du = − (t − u)f (t)dt du.
−1 0 −1 0

On en conclut que
Z 1
|e| = f (u)du − 2f (0) ≤ 1/3m(2) (f ).
−1

Un calcul analogue sur l’intervalle [xi , xi+1 ] donne

Z xi+1
xi + xi+1
f (t)dt = (xi+1 − xi )f ( ) + ei
xi 2
avec
(xi+1 − xi )3 (2)
|ei | ≤ mi (f ).
24
Sur l’intervalle [a, b] on obtient par sommation
Z b
(b − a)
f (t)dt = [f (x′0 ) + · · · + f (x′N −1 )] + E
a N
où
xi + xi+1
x′i =
2
et où
(b − a)3 (2)
|E| ≤ M (f ).
24N 2

4.2.3. Remarques concernant les deux méthodes précédentes. Ces

deux méthodes sont toutes les deux basées sur une interpolation de degré 0 (c’est à
dire par une constante). Cependant la deuxième a été en quelque sorte optimisée en
choisissant un bon point d’interpolation. Pourquoi dire que cette deuxième méthode
est meilleure que la première ? Ceci pour au moins deux raisons : d’une part l’erreur
E est en 1/N dans la première méthode et en 1/N 2 dans la deuxième, d’autre part
la première méthode donne un calcul exact pour la classe des polynômes constants
alors que la deuxième donne un calcul exact pour une classe plus vaste, celle des
polynômes de degré 1 (ce qui ne veut pas dire qu’au hasard des choix, on ne puisse
pas tomber sur des fonctions particulières autres, où l’une ou l’autre des formules
donne une valeur exacte). Nous serons amenés plus tard à préciser ces critères de
comparaison.
4.2. MISE EN ŒUVRE DE MÉTHODES INTERPOLATOIRES 37

1111111111
0000000000
0000000000
1111111111
0000000000
1111111111
0000000000
1111111111
0000000000
1111111111
00000000001
1111111111
−1

Fig. 3. Méthode des trapèzes

4.2.4. Méthode des trapèzes. Passons maintenant à une interpolation de

Lagrange de degré 1 (c’est-à-dire par une fonction affine) aux bornes de l’intervalle.

Ainsi sur l’intervalle [−1, 1] on approchera f (u) par

(u + 1)f (1) + (u − 1)f (−1)
.
2
Alors on aura Z 1
f (u)du = f (1) + f (−1) + e.
−1
Le calcul de l’erreur e s’effectue en introduisant l’erreur dûe à l’interpolation de
Lagrange, erreur calculée dans le chapitre portant sur l’interpolation :
(u + 1)f (1) + (u − 1)f (−1) 1
f (u) − ≤ |(u − 1)(u + 1)|m(2) (f ).
2 2
En intégrant cette inégalité sur [−1, 1] on obtient la majoration suivante de l’erreur
2 (2)
|e| ≤ m (f ).
3
Des calculs analogues sur l’intervalle [xi , xi+1 ] nous donnent
Z xi+1
xi+1 − xi
f (u)du = [f (xi ) + f (xi+1 )] + ei
xi 2
avec
(xi+1 − xi )3 (2)
|ei | ≤ mi (f ).
12
Enfin sur l’intervalle [a, b] on peut écrire
Z b
b − a f (x0 ) + f (xN )
f (u)du = + f (x1 ) + · · · + f (xN −1 ) + E
a 2 2
avec
(b − a)3 (2)
|E| ≤ M (f ).
12N 2
38 4. CALCUL NUMÉRIQUE DES INTÉGRALES DÉFINIES

4.2.5. Méthode des tangentes. Cette méthode consiste à approcher la fonc-

tion par sa tangente au point milieu de l’intervalle. C’est donc une interpolation
d’Hermite par un polynôme de degré 1.

1111111111
0000000000
0000000000
1111111111
0000000000
1111111111
0000000000
1111111111
0000000000
1111111111
0000000000
1111111111
0000000000
1111111111
0000000000
1111111111
00000000001
1111111111
−1

Fig. 4. Méthode des tangentes

Des considérations géométriques très simples montrent que cette méthode se

ramène exactement à la méthode du milieu. Ainsi on dispose d’une piste pour
expliquer que la méthode du milieu soit meilleure que celle des rectangles : en bien
choisissant le point d’interpolation de degré 0, tout se passe comme si on disposait
d’une interpolation de degré 1.
4.2.6. Méthode de Gauss. Comme dans le cas de l’interpolation de La-
grange de degré 0, nous allons essayer dans le cas de l’interpolation de Lagrange de
degré 1, de choisir au mieux les deux points d’interpolation de manière à améliorer
la méthode. Ainsi pour la méthode des trapèzes les point d’interpolation étaient aux
bornes de l’intervalle, pour la méthode de Gauss les points d’interpolation seront
des points bien choisis de l’intervalle.

On se place dans l’espace R1 [X] des polynômes de degré ≤ 1 sur R. Soient α et

β tels que −1 ≤ α < β ≤ 1. On note δα et δβ les formes linéaires sur R1 [X] définies
par δu (P ) = P (u). Si on note
X −β
Pα (X) =
α−β
et
X −α
Pβ (X) =
β−α
alors δu (Pv ) = δu,v (0 si u 6= v, 1 si u = v). On en conclut que (δα , δβ ) est une base
du dual de R1 [X]. Donc la forme linéaire I su R1 [X] définie par
Z 1
I(P ) = f (u)du
−1
se décompose sous la forme
I = λα δα + λβ δβ .
4.2. MISE EN ŒUVRE DE MÉTHODES INTERPOLATOIRES 39

Fig. 5. Méthode de Gauss

Cette formule appliquée aux polynômes Pα et Pβ donne les valeurs explicites des
constantes λα et λβ , ce qui permet d’écrire
2β 2α
I= δα + δβ .
β−α α−β
Ceci veut dire que pour tout polynôme P ∈ R1 [X] on a
Z 1
2β 2α
P (u)du = P (α) + P (β).
−1 β − α α −β
Pour une fonction f on va donc avec cette méthode écrire
Z 1
2β 2α
f (u)du = f (α) + f (β) + e,
−1 β−α α−β
et l’erreur e est nulle pour les polynômes de degré ≤ 1. Peut-on choisir les points α
et β pour que la formule reste exacte (e = 0) pour un degré supérieur ? La réponse
est donnée par :
Proposition 4.1. Il existe un couple et un seul de points α et β tels que pour
tout polynôme P de degré ≤ 3 on ait
Z 1
2β 2α
P (u)du = P (α) + P (β);
−1 β−α α−β
ces points sont
−1 1
α= √ β= √ .
3 3
De plus, quels que soient les points α et β il existe un polynôme de degré 4 qui ne
vérifie pas la formule précédente.
Preuve. le calcul se fait à partir de la formule exacte pour les polynômes de degré
≤ 1 en imposant que cette formule reste exacte pour le polynôme X 2 et le polynôme
X 3 (ce qui est nécessaire et suffisant pour que la formule reste exacte pour les
polynômes de degré ≤ 3). Un calcul simple nous donne les points α et β attendus
ainsi que les coefficients explicites de la formule.
40 4. CALCUL NUMÉRIQUE DES INTÉGRALES DÉFINIES

Proposition 4.2. Pour tout polynôme de degré ≤ 3 on a

Z 1
1 1
P (u)du = P (− √ ) + P ( √ ).
−1 3 3
La méthode de Gauss consiste donc à utiliser cette formule pour approcher
l’intégrale :
Z 1
1 1
f (u)du = f (− √ ) + f ( √ ) + e.
−1 3 3
La majoration de l’erreur e se fait maintenant en utilisant un polynôme de
degré ≤ 3 qui interpole f . Plus précisément on considère le polynôme d’Hermite H
de degré ≤ 3 tel que
1 1 1 1
H(− √ ) = f (− √ ) H( √ ) = f ( √ )
3 3 3 3
1 1 1 1
H ′ (− √ ) = f ′ (− √ ) H ′ ( √ ) = f ′ ( √ ).
3 3 3 3
Dans ces conditions Z 1
1 1
H(u)du = f (− √ ) + f ( √ ),
−1 3 3
donc Z Z Z
1 1 1
|e| = f (u)du − H(u)du ≤ |f (u) − H(u)|du.
−1 −1 −1
On sait alors (cf. chapitre sur l’interpolation) que sur [−1, 1]
2
1 2 1
|f (u) − H(u)| ≤ u − m(4) (f ),
24 3
d’où on tire par intégration
m(4) (f )
|e| ≤ .
135
Sur chaque intervalle [xi , xi+1 ] on obtient de manière analogue
Z xi+1
xi+1 − xi xi+1 + xi xi+1 − xi xi+1 + xi xi+1 − xi
f (t)dt = f + √ +f − √ +ei
xi 2 2 2 3 2 2 3
avec
5 (4)
xi+1 − xi mi (f )
|ei | ≤ .
2 135
Par sommation on obtient sur [a, b]
Z b N
X −1 Z xi+1
f (u)du = f (u)du + E
a i=0 xi

avec
(b − a)5 M (4) (f )
|E| ≤ .
N4 4320
Remarque 4.3. En fait bien qu’au départ on ait travaillé sur une interpolation
de Lagrange de degré 1, tout se passe ensuite de la même façon que pour une
interpolation d’Hermite de degré 3. C’est ce qui fait la performance de la méthode.
4.2. MISE EN ŒUVRE DE MÉTHODES INTERPOLATOIRES 41

−1 1

Fig. 6. Métode de Simpson

4.2.7. Méthode de Simpson. Nous partons maintenant de l’interpolation

de Lagrange de degré 2 (par une fonction parabolique) aux bornes de l’intervalle et
en un point γ tel que −1 < γ < 1. Nous employons une démarche analogue à celle
du paragraphe précédent. Notons δ−1 , δγ , δ1 les formes linéaires sur l’espace R2 [X]
des polynômes de degré ≤ 2 définies par
δu (P ) = P (u).
Les polynômes
(X − 1)(X − γ)
P−1 (X) = ,
2(1 + γ)
(X − 1)(X + 1)
Pγ (X) = ,
(γ − 1)(γ + 1)
(X + 1)(X − γ)
P1 (X) = ,
2(1 − γ)
vérifient δu (Pv ) = δu,v . On en déduit que (δ−1 , δγ , δ1 ) est une base du dual de R2 [X]
et donc que la forme linéaire I définie par
Z 1
I(P ) = P (t)dt
−1
se décompose sur cette base
I = λ−1 δ−1 + λγ δγ + λ1 δ1 .
Ceci veut dire que pour tout polynôme P de degré ≤ 2 on a
Z 1
P (t)dt = λ−1 P (−1) + λγ P (γ) + λ1 P (1).
−1
Par un bon choix du point γ on va voir que cette formule persiste pour les polynômes
de degré ≤ 3. On laisse au lecteur d’établir la proposition suivante
Proposition 4.4. Il existe un point γ et un seul tel que pour tout polynôme
de degré ≤ 3 on ait
Z 1
P (t)dt = λ−1 P (−1) + λγ P (γ) + λ1 P (1).
−1
42 4. CALCUL NUMÉRIQUE DES INTÉGRALES DÉFINIES

Cette valeur de γ est 0. De plus il existe un polynôme de degré 4 pour lequel cette
formule n’a pas lieu.
Ce point γ = 0 étant choisi on obtient la formule dite des 3 niveaux.
Proposition 4.5. Pour tout polynôme de degré ≤ 3
Z 1
1
P (t)dt = [P (−1) + P (1) + 4P (0)].
−1 3
La méthode de Simpson consiste donc à utiliser cette formule pour approcher
l’intégrale :
Z 1
1
[f (−1) + f (1) + 4f (0)] + e.
f (u)du =
−1 3
Soit H le polynôme de degré ≤ 3 tel que
H(−1) = f (−1) H(0) = f (0) H(1) = f (1) H ′ (0) = f ′ (0).
On a alors
Z 1 Z 1 Z 1
|e| = f (u)du − H(u)du ≤ |f (u) − H(u)| du.
−1 −1 −1

On sait (cf. chapitre sur l’interpolation) que sur [−1, 1]

1 2
|f (u) − H(u)| ≤ u |(u − 1)(u + 1)|m(4) (f )
4!
ce qui donne par intégration
m(4) (f )
|e| ≤ .
90
Sur les intervalles [xi , xi+1 ] on obtient
Z xi+1
xi+1 − xi xi + xi+1
f (u)du = f (xi ) + f (xi+1 ) + 4f + ei
xi 6 2
avec
5 (4)
xi+1 − xi mi (f )
|ei | ≤ .
2 90
Sur l’intervalle [a, b] tout entier on peut écrire
Z b " N −1 N −1 #
b−a X X xi + xi+1
f (u)du = f (x0 ) + f (xN ) + 2 f (xi ) + 4 f +E
a 6N i=1 i=0
2
avec
(b − a)5 M (4) (f )
|E| ≤ .
N4 2880
L’erreur est là aussi comme pour la méthode de Gauss en 1/N 4 , mais avec des
coefficients un peu moins bons. Cependant la formule
√ ne fait pas intervenir comme
dans celle de Gauss des nombres ”compliqués” ( 3) ce qui explique qu’elle ait été
plus prisée avant l’apparition des ordinateurs.
Peut on là aussi optimiser mieux la méthode en choisissant de manière astu-
cieuse les trois points d’interpolation ? Ceci est en partie l’objet de l’étude plus
générale du paragraphe suivant.
4.2. MISE EN ŒUVRE DE MÉTHODES INTERPOLATOIRES 43

4.2.8. Méthodes interpolatoires et polynômes orthogonaux. Soient

P0 , P1 , · · · , Pn les n + 1 premiers polynômes de Legendre. Rappelons que les po-
lynômes de Legendre forment l’unique suite de polynômes telle que
a) deg(Pi ) = i
R1
b) −1 Pi (u)Pj (u)du = 0 (i 6= j)
c) Pi (1) = 1.
Appelons a1 , · · · , an les racines de Pn ; elles sont distinctes et appartiennent à
l’intervalle [−1, 1]. Les formes linéaires (δai , · · · , δan ) définies par δai (P ) = P (ai )
forment une base du dual de l’espace Rn−1 [X] des polynômes de degré ≤ n − 1
(espace de dimension n). On en conclut que la forme linéaire I définie sur Rn−1 [X]
par
Z 1
I(P ) = P (u)du
−1
se décompose sur cette base sous la forme
I = λa1 δa1 + · · · + λan δan ,
ce qui veut dire que pour tout polynôme P de degré ≤ n − 1 on a
Z 1
P (u)du = λa1 P (a1 ) + · · · + λan P (an ).
−1

Le choix des racines du polynôme de Legendre Pn pour les points ai est un très
bon choix dans le mesure où
Proposition 4.6. Pour tout polynôme de degré ≤ 2n − 1 on a
Z 1
P (u)du = λa1 P (a1 ) + · · · + λan P (an ).
−1

Preuve. Il suffit de montrer que la formule est vraie pour les polynômes
d’une base de l’espace R2n−1 [X] des polynômes de degré ≤ 2n − 1. La famille
(P0 , · · · , Pn−1 , P0 Pn , · · · , Pn−1 Pn ) est une base de R2n−1 [X] (tous les degrés sont
représentés une fois et une seule). Il est clair que P0 , · · · , Pn−1 vérifient la formule
puisque ces polynômes sont de degré ≤ n − 1. Quand aux Pi Pn (0 ≤ i ≤ n − 1)
il estR clair en utilisant l’orthogonalité de la suite des polynômes de Legendre
1
que −1 Pi (u)Pn (u)du = 0 et que sachant que les ai sont les racines de Pn ,
δi (Pj Pn ) = Pj (ai )Pn (ai ) = 0. Ils vérifient donc eux ausi la formule.
Toute autre famille de points est moins bonne que la famille (ai ).
Proposition 4.7. Soient α1 , · · · , αn des points distincts de l’intervalle [−1, 1]
dont l’un au moins n’est pas racine du polynôme de Legendre de degré n. Alors il
existe un polynôme Q de degré ≤ 2n − 1 tel que
Z 1
Q(u)du 6= λα1 Q(α1 ) + · · · + λαn Q(αn ).
−1

Preuve. Posons Qn (X) = (X − α1 ) · · · (X − αn ).Par hypothèse, Qn n’est pas

proportionnel au polynôme de Legendre Pn , il n’est donc pas orthogonal au sous
R 1 polynômes de degré ≤ n − 1. Il existe un polynôme R de ce sous espace
espace des
tel que −1 R(u)Qn (u)du 6= 0. Il suffit alors de prendre Q = RQn .
44 4. CALCUL NUMÉRIQUE DES INTÉGRALES DÉFINIES

Cette propriété optimale des racines des polynômes de Legendre peut être uti-
lisée pour avoir une meilleure majoration de l’erreur quand on remplace l’intégrale
à calculer par l’intégrale du polynôme d’interpolation de Lagrange de degré ≤ n − 1
en n points de l’intervalle. C’est ainsi √ que nous
√ avons procédé pour la méthode
de Gauss, avec n = 2. Les point −1/ 3 et 1/ 3 sont les racines du polynôme de
2
Legendre de degré 2 (P2 (X) = 3X2 −1 ).

Dans le cas de l’interpolation de degré 2 (avec les notations précédentes n = 3),

nous pouvons obtenir mieux que par la méthode de Simpson et par la méthode de
Gauss. La méthode de Simpson en effet est une méthode utilisant l’interpolation de
degré 2, partiellement optimisée par le choix du point milieu de l’intervalle, mais
qui n’est pas complètement optimisée puisque les deux autres points d’interpolation
sont les bords de l’intervalle et non pas des racines du polynôme de Legendre de
degré 3. Le polynôme de Legendre de degré 3 est

5X 2 − 3X
P3 (X) = ,
2
ce qui donne pour racines
√ √
3 3
a1 = − √ a2 = 0 a3 = √ .
5 5
Un calcul simple donne alors pour tout polynôme P de degré ≤ 5
Z 1 " √ ! √ !#
1 3 3
P (u)du = 5P − √ + 8P (0) + 5P √ .
−1 9 5 5

Ainsi pour toute fonction f de classe C 5 on a

Z 1 " √ ! √ !#
1 3 3
f (u)du = 5P − √ + 8P (0) + 5P √ +e
−1 9 5 5

où |e| se majore en utilisant l’inégalité

1 2 2
|f (u) − P (u)| ≤ u (u − 3/5)2 m(6) (f ),
6!
ce qui par intégration fournit
1
|e| ≤ m(6) (f ).
15750
Sur chaque intervalle [xi , xi+1 ] on a une erreur
7 (6)
xi+1 − xi mi (f )
ei ≤ ,
2 15750

et sur [a, b] une erreur

(b − a)7 M (6) (f )
|E| ≤ .
N 6 15750 × 128
4.4. ACCÉLÉRATION DE CONVERGENCE - MÉTHODE DE ROMBERG 45

4.3. Présentation générale des quadratures élémentaires - Ordre d’une

méthode
Toutes les méthodes vues jusqu’à présent entrent dans le cadre général suivant
(méthodes de quadrature élémentaires). On a sur l’intervalle [xi , xi+1 ] une
formule du type
Z xi+1 Xli
f (u)du = (xi+1 − xi ) ωi,j f (ηi,j ) + ei ,
xi j=0

où
li
X
ηi,j ∈ [xi , xi+1 ] et ωi,j = 1.
j=0
On essaie d’adapter les paramètres pour que la méthode soit la meilleure possible.
Définition 4.8. On dit qu’une méthode de quadrature est d’ordre N si la
formule approchée est exacte pour les polynômes de degré ≤ N et inexacte pour
au moins un polynôme de degré N + 1.
P li
Remarquons que puisque j=0 ωi,j = 1, une méthode de quadrature
élémentaire est toujours au moins d’ordre 0.
Le lecteur est invité à préciser les paramètres qui correspondent aux diverses
méthodes décrites précédemment.

4.4. Accélération de convergence - Méthode de Romberg

La méthode que nous allons décrire maintenant part d’une méthode
d’accélération de convergence sur les suites.
4.4.1. Formule de Richardson. Soit u une suite de nombres réels conver-
geant vers un nombre réel a. On suppose que

un − a = λk n + O(k ′n )
où k et k sont des réels tels que 0 < |k ′ | < |k| < 1 et où λ est un réel non nul.
′

Alors
un+1 − a = λk n+1 + O(k ′n )
kun − ka = λk n+1 + O(k ′n )
et par suite
un+1 − kun
− a = O(k ′n ).
1−k
On construit de cette manière une suite
un+1 − kun
vn =
1−k
qui converge vers a plus vite que la suite initiale u.
Plus généralement on peut supposer que
un − a = λ1 k1n + λ2 k2n + · · · + λp kpn + O(kp+1
n
)
où
0 < |kp+1 | < |kp | < · · · < |k1 | < 1
46 4. CALCUL NUMÉRIQUE DES INTÉGRALES DÉFINIES

et où λ1 , · · · , λn sont des réels non nuls. Pour toute suite w, notons Rk w la suite
définie par
wn+1 − kwn
Rk wn = .
1−k
Alors
un+1 − k1 un
Rk1 un = .
1 − k1
Puisque
un+1 − a = λ1 k1n+1 + λ2 k2n+1 + · · · + λp kpn+1 + O(kp+1
n
)
et que
k1 un − k1 a = λ1 k1n+1 + λ2 k1 k2n + · · · + λp k1 kpn + O(kp+1
n
)
on obtient
un+1 − k1 un − (a − k1 a) = (λ′2 k2n + · · · + λ′p kpn )(1 − k1 ) + O(kp+1
n
)
où
(k2 − k1 ) (kp − k1 )
λ′2 = λ2 , · · · , λ′p = λp .
1 − k1 1 − k1
On en déduit que
Rk1 un − a = λ′2 k2n + · · · + λ′p kpn ) + O(kp+1
n
).
Par conséquent on peut itérer le procédé et on obtient alors
n
Rkp Rkp−1 · · · Rk1 un − a = O(kp+1 ).
Remarque 4.9. Reprenons la formule de Richardson à l’ordre 1. Il se peut que
k ne soit pas connu. On pose alors
un+1 − µn un
vn =
1 − µn
où
un+1 − un
µn = .
un − un−1
Dans ces conditions
µn − k = O ((k ′ /k)n ) ,
donc
vn − a = O(k ′n ).
Remarque 4.10. On peut donner une version continue de la formule de Ri-
chardson.
Soit f (y) telle que limy→0 f (y) = α0 . On suppose que
f (y) = α0 + α1 y + · · · + αp y p + O(y p+1 ).
Soit alors 0 < r < 1 et y0 > 0.
Formons la suite
Am,0 = f (rm y0 ),
et remarquons que
lim Am,0 = α0 ,
m→∞

Am,0 = α0 + α1 y0 rm + · · · + αp y0p (rp )m + O (rp+1 )m .
4.4. ACCÉLÉRATION DE CONVERGENCE - MÉTHODE DE ROMBERG 47

Donc en posant k1 = r, · · · , kp = rp on se ramène à la formule de Richardson

précédente. Plus précisément, pour 0 ≤ n ≤ p − 1 on définit
Am,n − rn+1 Am−1,n
Am,n+1 = .
1 − rn+1
Dans ces conditions
Am,p − α0 = O (rp+1 )m .
4.4.2. Méthode de Romberg. Nous allons montrer comment marche cette
méthode .
Découpons l’intervalle [a, b] par un partage équidistant en N = 2n morceaux et
posons h = (b − a)/2n . Appelons
n
2X −1
I2n = 1/2h[f (a) + f (b) + 2 f (xi )]
i=1
la valeur approchée de l’intégrale donnée par la méthode des trapèzes.
La formule d’Euler-Maclaurin nous donne
Z b
f (u)du − I2n = α(1/4)n + O ((1/8)n ) .
a
On peut donc appliquer la méthode de Richardson à l’ordre 1 avec k = 1/4 et
k ′ = 1/8. On pourra remarquer que dans ce cas on obtient la formule donnée
par la méthode de Simpson sur 2n−1 intervalles. Il faut appliquer la méthode de
Richardson à un ordre ≥ 3 pour trouver une formule qui n’est pas donnée par les
considérations des paragraphes précédents.
CHAPITRE 5

Analyse numérique des équations différentielles

5.1. Introduction
5.1.1. Position du problème. Nous nous intéressons ici aux méthodes à
un pas d’intégration des équations différentielles du premier ordre. Précisons le
problème que nous nous posons. Nous cherchons à résoudre numériquement le
problème de Cauchy

y ′ = f (x, y)
(5.1)
y(x0 ) = y0 .
Dans toute la suite on supposera qu’il existe un voisinage compact V du point
(x0 , y0 ) tel que la fonction f soit continue sur V et lipschitzienne par rapport à y
uniformément par rapport à x. C’est-à-dire qu’il existe une constante L telle que
pour tout x, tout y1 , tout y2 tels que (x, y1 ) ∈ V et (x, y2 ) ∈ V on ait
|f (x, y2 ) − f (x, y1 )| ≤ L|y2 − y1 |.
Nous noterons aussi
M = sup |f (x, y)|.
(x,y)∈V
Soit T > 0 tel que
[x0 − T, x0 + T ] × {y | |y − y0 | ≤ M T } ⊂ V.
On sait alors qu’il existe une solution unique qu’on notera φ du problème de Cauchy
(5.1) définie sur l’intervalle [x0 − T, x0 + T ].
Nous allons chercher à approcher la fonction φ sur l’intervalle [x0 , x0 + T ] (ou
sur l’intervalle [x0 − T, x0 ]).
Dans certains cas on sera amené à faire des hypothèses plus fortes sur la
régularité de la fonction f , ce qui entrainera une plus forte régularité de la fonction
solution φ.
5.1.2. Notations. On choisit un partage
x0 < x1 < · · · < xn = x0 + T
de l’intervalle d’étude [x0 , x0 + T ], et on pose
hk = xk+1 − xk .
On définit aussi
H= max hk .
0≤k≤n−1
La fonction φ étant la solution du problème (5.1) on pose pour 0 ≤ k ≤ n
yk = φ(xk ).
49
50 5. ANALYSE NUMÉRIQUE DES ÉQUATIONS DIFFÉRENTIELLES

Une méthode numérique à un pas permettra de calculer une approximation de

yk à partir d’une approximation de yk−1 .

5.2. Généralités sur les méthodes

Partons de la formule
Z h
φ(x + h) = φ(x) + φ′ (x + t)dt.
0

Une première idée pour calculer y1 = φ(x1 ) = φ(x0 +h0 ) est d’utiliser cette formule,
donc d’écrire
Z h0
y1 = y0 + φ′ (x0 + t)dt,
0
puis d’évaluer l’intégrale
Z h0
φ′ (x0 + t)dt
0
par une méthode de calcul numérique approché d’intégrales.

5.2.1. Méthode de la tangente d’Euler.

5.2.1.1. Description. Utilisons comme méthode de calcul la méthode de Rie-
mann. Nous obtenons alors
Z h0
φ′ (x0 + t)dt = h0 φ′ (x0 ) + O(h20 ),
0

ce qui donne
y1 = y0 + h0 φ′ (x0 ) + O(h20 ).
Ceci nous conduit à prendre pour approximation du point y1 le point u1 défini par
u1 = u0 + h0 f (x0 , u0 ),
où u0 est une approximation de y0 . Cette méthode nous donne au rang k
uk = uk−1 + hk−1 f (xk−1 , uk−1 ).
Cette méthode est appelée méthode de la tangente d’Euler en raison de l’in-
terprétation géométrique suivante.
5.2.1.2. Interprétation géométrique. On peut voir une équation différentielle
comme la donnée d’un champ de vecteurs : à tout point (x, y) on fait correspondre le
vecteur (1, f (x, y)) de coefficient directeur f (x, y). Le problème de Cauchy consiste
à trouver la ligne de champ qui passe par le point (x0 , y0 ). On la construit de
manière approchée en traçant une ligne polygonale partant de (x0 , y0 ) et suivant
tout d’abord la direction de la tangente de coefficient directeur f (x0 , y0 ) à la trajec-
toire en ce point. On arrive à un point (x1 , y1 ) et en ce point on prend la direction
indiquée par le champ de vecteurs, c’est-à-dire la droite de coefficient directeur
f (x1 , y1 ). Évidemment en général dès le point (x1 , y1 ) on a quitté la bonne trajec-
toire. Mais on peut espérer que si on fait de tous petits pas, on ne décolle pas trop
(cf. figure 1).

5.2.2. Métode d’Euler modifiée.

5.2. GÉNÉRALITÉS SUR LES MÉTHODES 51

x0 x0+h x0+2h

Fig. 1. Méthode de la tangente d’Euler

5.2.2.1. Description. En améliorant la méthode d’intégration de

Z h0
φ′ (x0 + t)dt
0

par la méthode des trapèzes par exemple

Z h0
h0 ′
φ′ (x0 + t)dt = (φ (x0 ) + φ′ (x0 + h0 )) + O(h30 ),
0 2
alors
h0
φ(x0 + h) = φ(x0 ) + (f (x0 , y0 ) + f (x0 + h0 , φ(x0 + h0 )) + O(h30 ),
2
soit encore
h0
y1 = y0 + (f (x0 , y0 ) + f (x1 , y1 )) + O(h30 ).
2
Malheureusement au second membre intervient la valeur y1 que l’on veut justement
calculer. On peut penser à utiliser une méthode d’approximations successives pour
calculer y1 : on injecte une valeur approchée de y1 dans le second membre de la
formule précédente et on obtient une nouvelle valeur approchée de y1 qu’on espère
meilleure dans le premier membre de la formule. Pour cela on part de la valeur
approchée v1 de y1 donnée par la méthode précédente de la tangente d’Euler

v1 = u0 + h0 f (x0 , y0 )
52 5. ANALYSE NUMÉRIQUE DES ÉQUATIONS DIFFÉRENTIELLES

x0 x1= x0+h x0+2h

Fig. 2. Méthode d’Euler modifiée

(dans ce calcul on considèrera que v0 = y0 ). Est-il la peine d’aller plus loin dans la
méthode des approximations successives ? Une évaluation de l’erreur nous donne
∂f
f (x1 , φ(x0 + h0 )) − f (x1 , v1 ) = (y1 − v1 ) (x1 , y1 ),
∂y
ce qui donne en utilisant l’évaluation de l’erreur de la méthode d’Euler
|f (x1 , φ(x0 + h0 )) − f (x1 , v1 )| = O(h20 ).
En conséquence, il n’est pas utile d’aller plus loin pour rester dans les limites d’une
erreur en 0(h30 ) et on obtient en définitive
h0
y1 = y0 + (f (x0 , y0 ) + f (x1 , v1 )) + 0(h30 ),
2
où
v1 = y0 + h0 f (x0 , y0 ).
En conclusion nous sommes amenés à prendre le schéma suivant

 u0 est une valeur initiale proche de y0

(5.2) vk = uk−1 + hk−1 f (xk−1 , uk−1 )
 uk = uk−1 + hk−1 (f (xk−1 , uk−1 ) + f (xk , vk ))

2
Cette méthode est appelée méthode d’Euler modifiée et admet l’interprétation
géométrique suivante.
5.2. GÉNÉRALITÉS SUR LES MÉTHODES 53

5.2.2.2. Interprétation géométrique. Si partant du point (x0 , y0 ) la méthode de

la tangente d’Euler faisant un bout de chemin sur la droite de coefficient directeur
f (x0 , y0 ) on arrive au point (x1 , v1 ) alors on fait la moyenne entre le coefficient
directeur en (x0 , y0 ) et le coefficient directeur en (x1 , v1 ) et finalement on emprunte
à partir de (x0 , y0 ) la droite ayant pour coefficient directeur cette moyenne (cf.
figure 2).

5.2.3. Généralisation. Le nombre u0 étant une valeur approchée de y0 , on

construit la suite (uk )k≥0 en partant de u0 et en calculant uk pour k ≥ 1 par une
formule récurrente du type :
(5.3) uk = uk−1 + hk−1 Φ(xk−1 , uk−1 , hk−1 ).
C’est ce qu’on appelle une méthode à un pas. Par exemple dans la méthode
d’Euler on prend :
Φ(x, y, h) = f (x, y)
et dans la méthode d’Euler modifiée :
1
Φ(x, y, h) = (f (x, y) + f (x + h, y + hf (x, y))) .
2
La question est alors la suivante : comment obtenir une fonction Φ intéressante.
5.2.3.1. Première idée : la formule de Taylor. On peut penser à la formule de
Taylor, ce qui nous donne :
hp0 (p)
y1 = y0 + hφ′ (x0 ) + · · · + φ (x0 ) + O(hp+1
0 ).
p!
Posons :
f (0) (x, y) = f (x, y),

∂(f (0) ) ∂(f (0) )

f (1) (x, y) = (x, y) + (x, y)f (x, y),
∂x ∂y
et plus généralement :
∂(f (k−1) ) ∂(f (k−1) )
f (k) (x, y) = (x, y) + (x, y)f (x, y).
∂x ∂y
Il est clair que :
f (k) (x, φ(x)) = φ(k+1) (x),
si bien que :

h h(p−1) (p−1)
y1 = y0 + h f (0) (x0 , y0 ) + f (1) (x0 , y0 ) + · · · + f (x0 , y0 ) + O(hp+1 ).
2 p!
On est donc amené à poser :
h (1) h(p−1) (p−1)
Φ(x, y, h) = f (0) (x, y) + f (x, y) + · · · + f (x, y).
2 p!
Il est visible que dans cette méthode, le calcul des dérivées successives peut s’avérer
compliqué. De plus l’erreur est difficile à contrôler et peut être catastrophique si les
dérivées de f ne sont pas ”petites”.
54 5. ANALYSE NUMÉRIQUE DES ÉQUATIONS DIFFÉRENTIELLES

5.2.3.2. Deuxième idée : une formule d’intégration. Cette idée, que nous avons
évoquée au début de cette section, consiste à partir de q > 0 nombres réels
c1 , c2 , · · · , cq distincts ou non, d’écrire des formules de quadrature numérique du
type :
Z ci Xq
(5.4) ψ(t)dt ≈ ai,j ψ(cj )
0 j=1

et :
Z 1 q
X
(5.5) ψ(t)dt ≈ bj ψ(cj ).
0 j=1

Rappelons que l’intervalle [x0 , x0 + T ] a été découpé sous la forme :

x0 < x1 · · · < xn = x0 + T.
Chaque sous-intervalle [xk , xk+1 ] où xk+1 = xk + hk contient lui-même un sous-
découpage basé sur q points. Ces q points sont les points :
xk,i = xk + ci hk ,
où 1 ≤ i ≤ q. Alors les formules de quadrature numérique choisies nous permettent
d’écrire : Z xk,i
φ(xk,i ) = yk + f (t, φ(t))dt,
xk
ce qui donne l’approximation :
q
X
φ(xk,i ) = yk + hk ai,j f (xk,j , φ(xk,j )).
j=1

On peut aussi écrire :

Z xk +hk
φ(xk + hk ) = yk + f (t, φ(t))dt,
xk

ce qui donne l’approximation :

q
X
yk+1 ≈ yk + hk bj f (xk,j , φ(xk,j )).
j=1

Ces considérations permettent de penser au schéma suivant :

 q
 X

 uk,i = u k + hk ai,j f (xk,j , uk,j )



 j=1
Xq
(5.6)

 uk+1 = u k + hk bj f (xk,j , uk,j )



 j=1

xk,j = xk + cj hk ,
où 
 0 ≤ k ≤ n−1
1 ≤ j ≤ q

1 ≤ i ≤ q.
En posant :
Kk,i = f (xk,i , uk,i )
5.2. GÉNÉRALITÉS SUR LES MÉTHODES 55

on obtient :
  
 q
X



 Kk,i
 = f xk,i , uk + hk ai,j Kk,j 
j=1
(5.7) q

 X


 uk+1
 = u k + hk bj Kk,j .
j=1

Le schéma ainsi décrit peut être explicite : le calcul de uk,i (ou de Kk,i ) ne fait
intervenir que des uk,j (des Kk,j ) avec j < i, c’est-à-dire déjà calculés, ou implicite
dans le cas contraire.
Il reste à savoir comment choisir les coefficients ai,j , bj , cj , c’est-à-dire quelles
formules de quadrature choisir. Toutes ces méthodes rentrent sous la dénomination
commune de méthode de Runge-Kutta. Cependant une de ces méthodes porte le
nom de méthode de Runge-Kutta classique.
5.2.4. Les méthodes de Runge-Kutta. Nous retrouvons bien entendu par
ce principe général les cas déjà traitées.
• Métode d’Euler. Dans ce cas on prend q = 1, c1 = 0, a1,1 = 0, b1 = 1.
Avec ces valeurs on retrouve effectivement la méthode d’Euler.
• Méthode d’Euler modifiée. Dans cette méthode aussi appelée méthode
de Heun on prend q = 2, c1 = 0, c2 = 1, a1,1 = 0, a1,2 = 0, a2,1 = 1, a2,2 = 0,
b1 = 1/2, b2 = 1/2. Les méthodes d’intégration utilisées sont d’une part la
méthode des rectangles pour la méthode d’intégration numérique (5.4) où
i = 2 (pour i = 1 il n’y a rien à faire) et celle des trapèzes pour la méthode
d’intégration numérique (5.5).
Toujours pour q = 2 on peut généraliser cette méthode en prenant un
nombre 0 < α ≤ 1 puis c1 = 0, c2 = α, a1,1 = 0, a1,2 = 0, a2,1 = 1, a2,2 = 0,
b1 = 1 − 1/(2α)), b2 = 1/(2α). Pour α = 1 on obtient la méthode de Heun,
pour α = 1/2 on obtient une méthode basée sur la méthode du milieu pour
la méthode d’intégration numérique (5.5) et la métode des rectangles pour
le calcul (5.4).
• Méthode de Runge-Kutta classique. On prend q = 4, c1 = 0, c2 =
1/2, c3 = 1/2, c4 = 1, a2,1 = 1/2, a3,2 = 1/2, a4,3 = 1, les autres ai,j sont
nuls, b1 = 1/6, b2 = 1/3, b3 = 1/3, b4 = 1/6.
On peut alors écrire le shéma de calcul de la façon suivante :


 Kk,1 = f (x
k , uk )

 hk hk



 Kk,2 = f xk + , uk + Kk,1

 2 2
 hk hk
(5.8) Kk,3 = f xk + , uk + Kk,2

 2 2



 Kk,4 = f (xk+1 ,uk + hk Kk,3 )

 Kk,1 Kk,2 Kk,3 Kk,4


 uk+1 = u k + hk + + + .
6 3 3 6
ANNEXE A

Représentation des nombres

A.1. Introduction
Nous nous intéressons à la façon de définir les objets que nous allons utiliser
dans les calculs sur machine pour représenter les nombres. Ceci passe par au moins
deux stades :
• La définition de la nature mathématique de ces objets
• Leur représentation interne
Nous regarderons les cas des réels et des entiers.

A.2. Les nombres réels

Les nombres réels sont approchés par les nombres à virgule flottante (en
simple précision, double précision, précisions étendues).
Les normes IEEE-754 et IEEE-854 (Institute of Electrical and Electronics
Engineers) √ définissent ces nombres ainsi que certaines façons d’opérer dessus
(+, −, ÷, ×, ).

A.2.1. Première approche. On fixe une base de numération b, un nombre

fixe p de digits, un signe S = (−1)s (avec s = 0 ou s = 1) un exposant e compris
entre deux valeurs fixées emin ≤ e ≤ emax . On considère alors les nombres de la
forme :

x = (−1)s x0 .x1 x2 · · · xp−1 be ,

où 0 ≤ xi < b.
la partie S est le signe, x0 x1 · · · xp−1 la mantisse et e l’exposant.

A.2.2. Les flottants normalisés. Un flottant normalisé est un réel non

nul qui peut s’écrire sous la forme précédente avec x0 6= 0. Les nombres réels qui
seront exactement représentables en machine seront 0 et les flottants normalisés.
Un réel pourra ne pas être exactement représenté pour les raisons suivantes : man-
tisse trop longue (infinie pour les nombres non b-adiques), exposant trop grand
(overflow) ou trop petit (underflow).

A.2.3. Cas concrets. Les bases utilisées sont 2 et 10. Chacune a ses avan-
tages.
La base 2 sera choisie quand surtout des pas de calcul et peu d’affichages
(langages de programmation généralistes par exemple).
La base 10 sera choisie si on a des affichages à chaque pas (calculette par
exemple) et donc beaucoup de traductions en décimal à faire.
La norme IEEE-754 régit le cas de la base 2.
57
58 A. REPRÉSENTATION DES NOMBRES

A.2.4. La base 2, simple précision. • La définition : b = 2, p = 24,

emax = 127, emin = −126 = −emax + 1. On introduit aussi l’exposant biaisé
E = e + emax .
• Le stockage en machine : 32 bits,

sE7 · · · E0 x1 · · · x23 .

Remarquons que x0 qui vaut 1 n’est pas écrit.

exemple : (−1.00110001000111111110011)2 × 2(11)2 se représente par les 32
bits suivants :
1 | 10000010 | 00110001000111111110011.

A.2.5. Bits non utilisés. Remarquons que 1 ≤ E ≤ 2emax = 254. Il reste

donc la valeur E = 0 et E = 255 qui sont non encore utilisées.
On utilise E = 0 avec une mantisse nulle pour 0.
On utilise E = 255 avec une mantisse nulle pour ±∞ (overflow).
On utilise E = 255 avec une mantisse non nulle pour ”Not a Number”.
On utilise E = 0 avec une mantisse non nulle pour l’underflow.

A.2.6. La base 2, double précision. b = 2, p = 53, emax = 1023, emin =

−1022 = −emax + 1, E = e + emax . Le stockage se fait sur 64 bits suivant le même
principe que pour la simple précision, avec 1 bits de signe, puis 11 bits d’exposant,
puis 52 bits de mantisse (x0 n’est pas écrit).

Il existe aussi une notion de double précision étendue avec un stockage sur
≥ 80 bits.

A.2.7. La base 10. Les calculatrices utilisent en général la base 10 de manière

à avoir un affichage facile. Les digits sont alors stockés en Décimal Codé Binaire.
Chaque digit décimal est codé sur 4 bits par sa valeur en binaire. Ainsi 7 par
exemple est codé 0111. On a de ce fait une petite perte (car sur 4 bits on pourrait
stocker 16 symboles).

A.2.8. Taille du stockage. L’exposant e vérifie en général −99 ≤ e ≤ 99 et

est stocké en DCB sur 1 octet. Les signes du nombre et de l’exposant sont stockés
sur 1/2 octet.
Si on dispose de n octets la mantisse aura une longueur p = 2n−3 (par exemple
si n = 8 on aurra 13 chiffres pour la mantisse). Attention il ne faut pas confondre
p avec le nombre de chiffres affichés qui est souvent plus petit.
Exercice : Monter un calcul qui permette de déduire la véritable taille de la
mantisse.

A.2.9. Arrondi. La norme IEEE définit quatre façons d’arrondir un résultat :

• Arrondi vers +∞ : x est arrondi au plus petit nombre représentable ≥ x.
• Arrondi vers −∞ : x est arrondi au plus grand nombre représentable ≤ x.
• Arrondi vers 0 : la valeur absolue de x est arrondie vers −∞ (le signe est
conservé).
• Arrondi au plus près : x est arrondi au nombre représentable le plus proche.
A.3. LES ENTIERS 59

A.2.10. Arrondi correct. On choisit un mode d’arrondi. Soit f une fonction

d’une ou de plusieurs variables réelles à valeurs réelles. On veut réaliser le calcul
de f de telle manière que si les valeurs de u = (u1 , · · · , uk ) sont des nombres
représentables alors la valeur approchée obtenue pour f (u) soit l’arrondi de la vraie
valeur de f (u). (Table Maker’s dilemma) √
La norme IEEE impose l’arrondi correct pour les fonctions de base +, −, ÷, ×,
et les conversions. Ceci n’est pas imposé pour sin, exp · · · .

A.2.11. Quelques difficultés. Supposons u, v, w des flottants. L’addition

n’est pas associative : notons F (x) le nombre flottant qui représente x.
F ((u + v) + w) = F (F (u + v) + w) ,

F (u + (v + w)) = F (u + F (v + w)) .
Or il se peut que F (u + v) = u et F (u + w) = u tandis que F (u + F (v + w)) 6= u.
Par exemple si b = 10, p = 11 on prend u = 1, v = w = 5 × 10−11 .

A.3. Les entiers

On veut représenter des nombres entiers en machine. Supposons par exemple
que nous ayons un octet pour le faire. Avec un octet on dispose de 28 =
256 écritures différentes, ces écritures pouvant être considérées comme les
développements binaires des entiers de l’intervalle {0..255}.
Comme on veut répartir équitablement les entiers que l’on représente entre des
entiers positifs et des entiers négatifs, on décide de s’intéresser aux 256 entiers de
l’intervalle {−128..127}. Il convient donc d’établir une bijection qui permette des
calculs commodes, entre l’intervalle {−128..127} des entiers qu’on veut représenter,
et l’intervalle {0..255} des représentations.
En résumé, à tout entier x de l’intervalle {−128..127} on va faire correspondre sa
représentation R(x) qui sera un entier de l’intervalle {0..255}. En outre comme R(x)
doit être stocké dans la mémoire d’une machine, on regardera plus spécialement les
propriétés du développement binaire (sur un octet) de R(x).

A.3.1. Représentation en ”complément à 2”. Rappelons que si n est un

entier, n mod 256 est le reste de la division de n par 256 ou encore l’unique entier
m tel que 0 ≤ m ≤ 255 et m congru à n modulo 256.
Notons I l’intervalle {−128..127} et J l’intervalle {0..255}. Soit R l’application
de I dans J définie par
R(x) = x mod 256.

A.3.2. Premières propriétés. a) Montrer que R est une application bi-

jective.
b) Calculer R(0), R(100), R(127), R(−1), R(−100), R(−128). Donner les
développements binaires des résultats obtenus.
c) Calculer R(x) en fonction de x.
d) Déterminer l’image par R de l’ensemble des x ≥ 0 de I ainsi que l’image
par R de l’ensemble des x < 0 de I.
60 A. REPRÉSENTATION DES NOMBRES

A.3.3. L’opposé. Comment reconnaı̂tre sur le développement binaire de

R(x) le signe de x ?
e) On suppose que x ∈ I \ {−128, 0} . Calculer R(−x) en fonction de R(x).
On constatera que R(−x) = (255 − R(x)) + 1. En déduire un algorithme simple
permettant de calculer le développement binaire de R(−x) à partir de celui de R(x)
(algorithme dit de complémentation à 2).
f) Pour écrire en binaire la représentation R(x) d’un entier x de l’intervalle
I on applique la stratégie suivante :
• Si x ≥ 0, on développe x en binaire.
• Si x < 0, on développe −x en binaire et on applique l’algorithme de
complémentation à 2 (cf. e) ).
Appliquer cette méthode pour calculer l’écriture binaire de R(18), R(−20).
A.3.4. Addition des entiers et représentation. Soit T l’application de
P∞ P7
N dans {0..255} qui à tout n = j=0 aj 2j fait correspondre T (n) = j=0 aj 2j
(troncature limitée aux 8 premiers bits).
a) Quel est le lien entre n mod 256 et T (n) ?
b) Montrer que si x1 , x2 , x1 + x2 sont des éléments de I alors

R(x1 + x2 ) = R(x1 ) + R(x2 ) mod 256,

R(x1 + x2 ) = T R(x1 ) + R(x2 ) .
c) Il est bien entendu que l’addition de deux éléments de I ne sera valide que
si le résultat est aussi dans I. Comme la machine ne connaı̂t que les représentants
des nombres qu’on additionne, nous mettons en place ici une méthode (bien adaptée
aux circuits électroniques) qui opère sur les représentations et permette à la fois de
détecter si l’opération d’addition est valide et de calculer dans ce cas le résultat.
Soient x1 et x2 deux éléments de I. Soit C le carry, retenue du 8ieme bit vers
l’exterieur, et α la retenue du 7ieme bit vers le 8ieme, obtenus en faisant
l’addition R(x1 ) + R(x2 ). On pose V = C ⊕ α.
• Si V = 0 l’addition est valide et R(x1 + x2 ) s’obtient en faisant en binaire
l’addition de R(x1 ) avec R(x2 ) et en négligeant tout débordement au delà
du huitième bit (cf. b) ).
• Si V = 1 l’addition n’est pas valide.
En effet : c1) Supposons 0 ≤ x1 ≤ 127 et 0 ≤ x2 ≤ 127. Examiner les deux
cas x1 + x2 ≤ 127 (opération valide) et x1 + x2 > 127 (opération invalide), et dans
chaque cas calculer C, α, V .
c2) Supposons 0 ≤ x1 ≤ 127 et −128 ≤ x2 < 0 (opération toujours valide).
On examinera suivant les valeurs possibles de R(x1 ) + R(x2 ) quelles sont les valeurs
possibles de C, α, V .
c3) Supposons −128 ≤ x1 < 0 et −128 ≤ x2 < 0. Examiner les deux
cas x1 + x2 < −128 (opération invalide) et x1 + x2 ≥ −128 (opération valide), et
dans chaque cas calculer C, α, V . (Indication : pour calculer α on pourra regarder si
l’addition des deux nombres de 7 bits (R(x1 ) − 128) et (R(x2 ) − 128) a une retenue
vers le huitième bit.)
c4) En conclure la validité de l’algorithme annoncé.

Vous aimerez peut-être aussi

X Racine
Pas encore d'évaluation
X Racine
7 pages
Maths 5
Pas encore d'évaluation
Maths 5
33 pages
Equations 07
Pas encore d'évaluation
Equations 07
6 pages
Epuchap1 PDF
Pas encore d'évaluation
Epuchap1 PDF
6 pages
Résolution d'équations non linéaires
100% (1)
Résolution d'équations non linéaires
7 pages
Cours Nonlin
Pas encore d'évaluation
Cours Nonlin
15 pages
Modulecalcul Num
Pas encore d'évaluation
Modulecalcul Num
41 pages
Chapitre 1 Résolution Numérique Des Équations Non Linéaires
Pas encore d'évaluation
Chapitre 1 Résolution Numérique Des Équations Non Linéaires
7 pages
Analyse Ch3 TD
Pas encore d'évaluation
Analyse Ch3 TD
12 pages
Résolution numérique d'équations
100% (1)
Résolution numérique d'équations
12 pages
COURS Calcul Numerique 2023-2024 UTA
Pas encore d'évaluation
COURS Calcul Numerique 2023-2024 UTA
53 pages
CoursMethodesNumérique Zougab
Pas encore d'évaluation
CoursMethodesNumérique Zougab
34 pages
Bernoussi PDF
Pas encore d'évaluation
Bernoussi PDF
18 pages
Cours Et TD Methode Numerique 2ST 2020
100% (2)
Cours Et TD Methode Numerique 2ST 2020
29 pages
Approximation Num Erique Des Solutions D' Equations: Ecole Nationale Des Sciences Appliqu Ees - F' Es
Pas encore d'évaluation
Approximation Num Erique Des Solutions D' Equations: Ecole Nationale Des Sciences Appliqu Ees - F' Es
6 pages
An Non Lineaire 2
Pas encore d'évaluation
An Non Lineaire 2
86 pages
Calculs Approches Point Fixe
Pas encore d'évaluation
Calculs Approches Point Fixe
10 pages
Fonctions Nonlineaires
Pas encore d'évaluation
Fonctions Nonlineaires
9 pages
Cours Analyse Numerique
Pas encore d'évaluation
Cours Analyse Numerique
20 pages
Chapitre1 20
Pas encore d'évaluation
Chapitre1 20
10 pages
Méthodes pour Trouver les Zéros d'Équations Non-Linéaires
Pas encore d'évaluation
Méthodes pour Trouver les Zéros d'Équations Non-Linéaires
6 pages
ANA NUM (1) - Copie
Pas encore d'évaluation
ANA NUM (1) - Copie
51 pages
Cours Methode Numerique 2ST 2024
Pas encore d'évaluation
Cours Methode Numerique 2ST 2024
28 pages
Résolution numérique des équations
Pas encore d'évaluation
Résolution numérique des équations
9 pages
Méthode de Newton et convergence
Pas encore d'évaluation
Méthode de Newton et convergence
13 pages
MP3 Preparation For The Day and I Love You Too Baby Girl I Love
Pas encore d'évaluation
MP3 Preparation For The Day and I Love You Too Baby Girl I Love
6 pages
Poly MNRO 2022-23
Pas encore d'évaluation
Poly MNRO 2022-23
138 pages
Amar D Amar - Djema@univ-Bejaia - DZ 13 Février 2022
Pas encore d'évaluation
Amar D Amar - Djema@univ-Bejaia - DZ 13 Février 2022
14 pages
Chap2 Test
100% (1)
Chap2 Test
16 pages
Résolution d'équations non linéaires
Pas encore d'évaluation
Résolution d'équations non linéaires
10 pages
Algorithmes pour Équations Non Linéaires
50% (2)
Algorithmes pour Équations Non Linéaires
10 pages
Numerique 2023
Pas encore d'évaluation
Numerique 2023
9 pages
Méthodes pour trouver les racines d'une fonction
Pas encore d'évaluation
Méthodes pour trouver les racines d'une fonction
7 pages
Résolution d'Équations Non Linéaires
Pas encore d'évaluation
Résolution d'Équations Non Linéaires
2 pages
1er Chapitre Analyz Numérik
Pas encore d'évaluation
1er Chapitre Analyz Numérik
25 pages
Chapitre 2
Pas encore d'évaluation
Chapitre 2
7 pages
Cours de Méthodes Numériques et Recherche Opérationnelle
Pas encore d'évaluation
Cours de Méthodes Numériques et Recherche Opérationnelle
58 pages
Méthode Numérique de Base
Pas encore d'évaluation
Méthode Numérique de Base
44 pages
CALPHA - Mathematiques Et Raisonnement Logique - 2012
Pas encore d'évaluation
CALPHA - Mathematiques Et Raisonnement Logique - 2012
19 pages
Chapitre 1 - Equations Non Linéaires
Pas encore d'évaluation
Chapitre 1 - Equations Non Linéaires
9 pages
Chapitre III PDF
Pas encore d'évaluation
Chapitre III PDF
45 pages
Résumé de l'analyse numérique
Pas encore d'évaluation
Résumé de l'analyse numérique
9 pages
CNAM-Calcul Scientifique A1 - Résolution Des Fonctions
Pas encore d'évaluation
CNAM-Calcul Scientifique A1 - Résolution Des Fonctions
10 pages
Approximations et Convergences
Pas encore d'évaluation
Approximations et Convergences
56 pages
Nonlineaire 1
Pas encore d'évaluation
Nonlineaire 1
84 pages
MN - Chapitre 1 - Equation - Non - Lineaire
Pas encore d'évaluation
MN - Chapitre 1 - Equation - Non - Lineaire
9 pages
Notes Dec Ours
Pas encore d'évaluation
Notes Dec Ours
37 pages
Chapitre4 Ana Num Eq Non Lineaires 22 23
Pas encore d'évaluation
Chapitre4 Ana Num Eq Non Lineaires 22 23
7 pages
Résolution numérique d'équations
Pas encore d'évaluation
Résolution numérique d'équations
9 pages
Cours Saisie AN1 Word
Pas encore d'évaluation
Cours Saisie AN1 Word
106 pages
Chapitre 01 Méthode Numérique
Pas encore d'évaluation
Chapitre 01 Méthode Numérique
10 pages
Chapitre 1
Pas encore d'évaluation
Chapitre 1
22 pages
Equations Non Lineaires
Pas encore d'évaluation
Equations Non Lineaires
7 pages
CalculNum2425 VF
Pas encore d'évaluation
CalculNum2425 VF
56 pages
Support de Cours Chapitre 4
Pas encore d'évaluation
Support de Cours Chapitre 4
14 pages
Cours d'Analyse Numérique 2
Pas encore d'évaluation
Cours d'Analyse Numérique 2
15 pages
Emploi du temps automne 2023-2024 ENCG
Pas encore d'évaluation
Emploi du temps automne 2023-2024 ENCG
9 pages
Chap 3 Les Resistances Correction
Pas encore d'évaluation
Chap 3 Les Resistances Correction
4 pages
Rôle de L'enseignant Montessori
Pas encore d'évaluation
Rôle de L'enseignant Montessori
4 pages
Cours Et Exercice Calcul Vectoriel Dans Le Plan 1
Pas encore d'évaluation
Cours Et Exercice Calcul Vectoriel Dans Le Plan 1
4 pages
Matrices de permutation et symétriques
Pas encore d'évaluation
Matrices de permutation et symétriques
11 pages
Typologie et gestion des conflits
Pas encore d'évaluation
Typologie et gestion des conflits
10 pages
TAUVEL CorpsGalois Exercicesv03
Pas encore d'évaluation
TAUVEL CorpsGalois Exercicesv03
21 pages
PNL - SM276 - 05 - Thiry
Pas encore d'évaluation
PNL - SM276 - 05 - Thiry
3 pages
Les Reconstitutions Corono-Radiculaires
100% (1)
Les Reconstitutions Corono-Radiculaires
14 pages
Communication Optique et Multiplexage
Pas encore d'évaluation
Communication Optique et Multiplexage
3 pages
2ème Devoir Du 2ème Semestre Histoire-Geographie 3ème 2020-2021 Ceg Agassa-Godomey
Pas encore d'évaluation
2ème Devoir Du 2ème Semestre Histoire-Geographie 3ème 2020-2021 Ceg Agassa-Godomey
2 pages
Série 1 - 2 Bac SM-2026
Pas encore d'évaluation
Série 1 - 2 Bac SM-2026
5 pages
Theories Du Developpement. Version 2
Pas encore d'évaluation
Theories Du Developpement. Version 2
16 pages
MVA006 ndc08
Pas encore d'évaluation
MVA006 ndc08
10 pages
Terminologie Économique et Méthodologie
Pas encore d'évaluation
Terminologie Économique et Méthodologie
6 pages
Cisaillement Simple: Concepts et Applications
Pas encore d'évaluation
Cisaillement Simple: Concepts et Applications
12 pages
Signo Cafop
Pas encore d'évaluation
Signo Cafop
5 pages
Lutter Contre La Désinformation - Cortecs
Pas encore d'évaluation
Lutter Contre La Désinformation - Cortecs
33 pages
Projet Fin D'etude
Pas encore d'évaluation
Projet Fin D'etude
12 pages
3e SQ4 Sujet R
Pas encore d'évaluation
3e SQ4 Sujet R
7 pages
AUTRES Directives élaboration-PA-version-finale
Pas encore d'évaluation
AUTRES Directives élaboration-PA-version-finale
58 pages
T.P. Cellule Photovoltaïque: Notice Notice
Pas encore d'évaluation
T.P. Cellule Photovoltaïque: Notice Notice
12 pages
Langue Voix Active Et Voix Passive Candide
Pas encore d'évaluation
Langue Voix Active Et Voix Passive Candide
2 pages
TD 1 - Eecs
100% (1)
TD 1 - Eecs
1 page
Dossier - KLP - Emmanuell Ferran
Pas encore d'évaluation
Dossier - KLP - Emmanuell Ferran
17 pages
Corrigé de La Série 1 de TD
Pas encore d'évaluation
Corrigé de La Série 1 de TD
3 pages
Lettre Mission - Inspections Modele Type
Pas encore d'évaluation
Lettre Mission - Inspections Modele Type
2 pages
Phobie des Maladies chez Infirmiers
Pas encore d'évaluation
Phobie des Maladies chez Infirmiers
11 pages
Programme D'assainissement Des Petites Communes
Pas encore d'évaluation
Programme D'assainissement Des Petites Communes
50 pages
Capteur Metrologie Cours 23 11
Pas encore d'évaluation
Capteur Metrologie Cours 23 11
189 pages