Méthode d’analyse numérique.

Pascal Viot

To cite this version:

Pascal Viot. Méthode d’analyse numérique.. 2006. cel-00092946

HAL Id: cel-00092946

[Link]
Submitted on 12 Sep 2006

HAL is a multi-disciplinary open access

L’archive ouverte pluridisciplinaire HAL, est
archive for the deposit and dissemination of sci-
destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub-
scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from
émanant des établissements d’enseignement et de
teaching and research institutions in France or
recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers.
publics ou privés.
 M´ethodes d’analyse num´erique
Cours du DEA “Mod´elisation Dynamique et
Statistique des Syst`emes Complexes”

Pascal Viot
Laboratoire de Physique Th´eorique des Liquides, Boˆıte
121, 4, Place Jussieu, 75252 Paris Cedex 05
Email : viot@[Link]

17 janvier 2003
 Ce cours est une introduction aux m´ethodes d’analyse num´erique
tr`es lar- gement utilis´ees en physique afin de r´esoudre les ´equations alg
´ebriques ou dif- f´erentielles que l’on rencontre dans la mod´elisation de
ph´enom`enes physiques, chimiques ou biologiques.
Ce domaine particuli`erement vaste n´ecessite simultan´ement des
connais- sances math´ematiques, informatiques et physiques. De larges
classes de pro- bl`emes num´eriques sont abord´ees dans ces notes et
montrent la n´ecessit´e de bien caract´eriser les propri´et´es math
´ematiques du probl`eme consid´er´e afin de choisir la m´ethode num´erique
la mieux adapt´ee pour le traitement num´erique.

2
Chapitre 1

Int´egration et
sommes discr`etes

Contenu
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Les m´ethodes de Cˆotes . . . . . . . . . . . . . . . . . 4
1.2.1 Trap`eze ...................................................................................5
1.2.2 Simpson.....................................................................................5
1.3 M´ethode de Romberg . . . . . . . . . . . . . . . . . . 6
1.4 M´ethodes de Gauss . . . . . . . . . . . . . . . . . . . 6
1.5 Int´egrales multiples . . . . . . . . . . . . . . . . . . . 9

1.1 Introduction
Les sciences physiques se fixent l’objectif de pr´edire les ph´enom`enes `a
partir de la connaissance d’un nombre fini de grandeurs microscopiques
intervenant dans la mod´elisation. Celle-ci doit, dans un premier temps, d
´ecrire les faits exp´erimentaux observ´es et, dans un second temps,
permettre d’en pr´edire de nouveaux. La description des ph´enom`enes peut
ˆetre soit microscopique, soit ph´enom´enologique.
La mod´elisation se divise naturellement en trois ´etapes : une premi`ere
´etape consiste `a d´eterminer les param`etres microscopiques essentiels qui
interviennent pour d´ecrire le ph´enom`ene puis `a choisir le mod`ele adapt´e
pour d´ecrire le ph´e- nom`ene. La second ´etape consiste `a ´etablir les
´equations (tr`es souvent diff´eren- tielles) qui d´ecrivent le ph´enom`ene.
La troisi`eme ´etape consiste `a r´esoudre les
´equations pr´ec´edemment ´etablies afin de donner des r´eponses
quantitatives. Ces r´eponses permettent de valider ou d’invalider le mod`ele,
soit en comparant les pr´edictions avec l’exp´erience, soit en analysant la
coh´erence interne du mod`ele
`a travers ses pr´edictions.
Dans la situation la plus optimiste ou` il a ´et´e possible de r´ealiser les trois
´etapes pr´ec´edentes, l’obtention d’un r´esultat num´erique n´ecessite au
moins le calcul d’une int´egrale simple. Pour r´ealiser ce ’simple’ travail, il
3
existe bien souvent plusieurs m´ethodes et la taˆche principale consiste `a s
´electionner celle

4
 Int´egration et sommes discr`etes

qui est la mieux adapt´ee pour satisfaire l’objectif. Ce principe est tr`es g
´en´eral et s’applique `a l’ensemble des probl`emes num´eriques que nous
allons aborder tout au long de ce cours. Ce cours a donc pour objectifs de
pr´esenter quelques algorithmes associ´es `a chacun des probl`emes rencontr
´es. Il existe des logiciels gratuits ou payants qui utilisent les m´ethodes
num´eriques que nous allons ren- contrer dans ce cours ; la connaissance de
ces m´ethodes est donc un pr´ealable pour choisir la mani`ere de r´esoudre
un probl`eme num´erique donn´e le plus effi- cacement possible.
Le cas typique d’une int´egrale d´efinie `a ´evaluer est
∫ b

a
I f (x)dx.
Comme l’´evaluation de la fonction f pour une infinit´e de points est
impossible, l’int´egration num´erique consiste `a remplacer l’int´egrale Eq.
(1.1) par une somme discr`ete sur un nombre fini de points.
N
Σ
IN = aif (xi) (1.2)
i=1

ou` ai et xi sont des variables que nous allons pr´eciser dans la suite. Pour
que l’´evaluation num´erique soit correcte, il est n´ecessaire d’imposer que
toute m´e- thode v´erifie que
lim IN = I (1.3)
N
→∞

Au del`a de la v´erification de ce crit`ere Eq. (1.3), la qualit´e d’une

m´ethode sera ´evalu´ee par la mani`ere dont la convergence vers le r´esultat
exact s’effectue. Dans la suite nous allons consid´erer des fonctions de classe
C1 1 (continuˆment d´erivable) sur le support [a, b]2 . La derni`ere restriction
impos´ee `a la fonction est que

|f J (x)| < K ∀x ∈ [a, b] (1.4)

ou` K est une constante finie. Cela revient `a supposer que la d´eriv´ee de la
fonction f n’est jamais singuli`ere sur le support [a, b].

1.2 Les m´ethodes de Coˆtes

Math´ematicien anglais, contemporain et collaborateur de Newton,
Roger Coˆtes s’est int´eress´e `a des m´ethodes de calculs num´eriques et
exacts pour l’in- t´egration et explique que les m´ethodes suivantes portent
son nom.
1
On peut calculer l’int´egrale d’une fonction plus g´en´erale a` condition que cette
fonction ait la mˆeme mesure ( mesure d´efinie au sens de Lebesgue) qu’une fonction de
classe C 1 , c’est- a`-dire que les fonctions ne diff`erent que sur un ensemble de mesure
nulle. Comme exemple simple d’ensemble de mesure nulle, citons les ensembles d
´enombrables.
2
Si la fonction est continuˆment d´erivable par morceaux sur un nombre fini d’intervalles
sur l’intervalle [a, b], on peut se ramener au cas pr´ec´edent sur chaque intervalle, et

4
donc ´evaluer l’int´egrale sur l’ensemble de l’intervalle [a, b].

5
 1.2 Les m´ethodes de Coˆtes

Les m´ethodes les plus simples que l’on peut utiliser pour calculer une
int´e- grale simple sont celles ou` les abscisses sont choisies de mani`ere r
´eguli`erement espac´ees. Si on a N + 1 abscisses, on rep`ere celles-ci
simplement par la relation
xi = x0 + ih (1.5)
avec x0 = a, xN = b et h est appel´e le pas de l’int´egration. Pour
simplifier les notations, on pose
fi = f (xi) (1.6)

1.2.1 Trap`eze
La m´ethode des trap`ezes consiste `a approximer la fonction entre
deux abs- cisses successives par une droite, ce qui donne.
∫ xi+1
h 3 JJ

f (x)dx = (fi + fi+1) + ).

x 2
Le terme d’erreur indique la qualit´e de l’´evaluation de l’int´egration et d
´epend de mani`ere cubique du pas d’int´egration ; f JJ se r´ef`ere `a un point
situ´e `a l’int´erieur de l’intervalle. Pour que cette m´ethode converge
rapidement il est n´ecessaire de choisir un pas h inf´erieur `a f JJ . A noter que
cette formule devient exacte quand la fonction est un polynˆome de degr´e 1
sur l’intervalle [x1 , x2 ].
Sur l’intervalle [a, b], on a
N
∫b (b − a)3f
h −1 h
f (x) = Σfi + 2 (f0 + fN ) + JJ
(1.8)
a
i=1 O
N2
ou` on a utilis´e que h = (b − a)/N

1.2.2 Simpson
La m´ethode de Simpson consiste `a remplacer la fonction par un
polynˆome de degr´e 2 sur un intervalle constitu´e de trois abscisses cons
´ecutives
∫ xi+2
f (x)dx 1 fi +4 fi+1 + 5 (4)
3 3 3 ).
1 +
x
Il se trouve que cette formule est exacte jusqu’`a des polynˆomes de degr´e
3 ce qui implique que l’erreur d´epende de h `a la puissance 5.
On peut aussi d´eterminer la formule `a 4 points qui est aussi exacte pour
les polynˆomes de degr´e 3. Cette formule s’appelle aussi Simpson 3/8 `a
cause des coefficients du d´eveloppement
∫ xi+3
3 9 9 3 5 (4)
f (x)dx fi + fi+1 + fi+2 + + ).
x 8 8 8 8
Sur un intervalle complet, en choisissant N pair, la m´ethode de
Simpson donne une estimation de l´int´egrale sur l´intervalle [a, b].

6
 
∫ h N/2−1
bf (x) = f + f + 2 Σ (2f2i−1 + f2i ) + O 1
(1.11)
0 N N
a 3 i=1 4

7
Int´egration et sommes discr`etes

La m´ethode de Simpson est donc de deux ordres de grandeur plus

efficace que la m´ethode des trap`ezes. Mais il est possible de mieux utiliser
la m´ethode des trap`ezes. Sachant que cette derni`ere m´ethode converge
en 1/N 2 , on peut
´evaluer l’int´egrale deux fois sur le mˆeme intervalle par la m´ethode des
trap`ezes ; la premi`ere fois avec N points et la seconde avec 2N points, pui
en combinant les deux r´esultats de la mani`ere suivante
4 1
S= S2N − SN (1.12)
3 3
Sachant que le d´eveloppement asymptotique de la m´ethode des trap`ezes est
une fonction paire de 1/N 2 , on en d´eduit que la formule (1.12) donne une
estimation de l’int´egrale en 1/N 4 , et ce r´esultat redonne une ´evaluation
analogue `a la m´ethode de Simpson.

1.3 M´ethode de Romberg

L’id´ee de la m´ethode de Romberg s’inspire directement de la remarque
faite au paragraphe pr´ec´edent. Si on calcule successivement par la m´ethode
−
des tra- p`ezes les int´egrales avec un nombres de points N/2k , N/2k 1 , . . . ,
N , on peut rapidement avoir une estimation de l’int´egrale avec une erreur en
(1/N 2k ) ou` k est le nombre de fois que l’on a calcul´e l’int´[Link] formule
it´erative utilis´ee est la suivante
Sk+1 (Sk(h) − Sk(2h)) (1.13)
(h) = (h) + (4k − 1)
Sk
Le tableau 1.1 r´esume la proc´edure r´ecursive employ´ee dans la m´ethode
de Romberg.

Pas Trap`eze Simpson Boole troisi`eme am

s ´elioration
h S1(h) S2(h) S3(h) S4(h)
2h S1(2h) S2(2h) S3(2h)
4h S1(4h) S2(4h)
8h S1(8h)

TAB. 1.1 – Table de Romberg

1.4 M´ethodes de Gauss

Dans les m´ethodes pr´ec´edentes, nous avons vu qu’en changeant les
coeffi- cients de pond´eration des valeurs de la fonction `a int´egrer aux
abscisses r´e- guli`erement espac´ees, on pouvait grandement am´eliorer la
convergence de la m´ethode. Les m´ethodes de Gauss ajoutent aux m
´ethodes pr´ec´edentes de pou- voir utiliser des abscisses non r´eguli`erement

8
espac´ees.
Soit W (x) une fonction strictement positive sur l´intervalle [a, b], appel
´ee fonction de poids, on choisit une suite de points xi telle que l’int´egrale
soit

9
 1.4 M´ethodes de Gauss

approch´ee par une somme discr`ete de la forme

b
∫ ΣN
W (x)f (x)dx = wifi (1.14)
a

i=1

Quand les abscisses sont r´eguli`erement espac´ees, les coefficients

inconnus sont les poids wi ; cela implique que pour N points d´int
´egration, on peut obtenir une ´evaluation exacte de l´int´egrale jusqu’`a un
polynˆome de degr´e N 1 si N est pair (trap`ezes) et N si N est impair −
(Simpson). Les m´ethodes de Gauss utilise le fait si les abscisses et les poids
sont des inconnues `a d´eterminer ; la formule d´int´egration de Gauss `a N
points devient exacte jusqu’`a des polynˆomes de degr´e 2N 1, ce qui
augmente la −pr´ecision de l´´evaluation sans qu’il soit n´ecessaire d
´augmenter le nombre de points `a calculer.
Pour d´eterminer ces 2N param`etres, on s´appuie sur la construction
de polynˆomes orthogonaux. Le principe de cette construction remonte `a
Gauss et Jacobi et a ´et´e largement d´evelopp´e par Christoffel. Soit un
intervalle (a, b), on introduit le produit scalaire de deux fonctions f et g avec
la fonction de poids W par :
∫b
< f |g >≡
W (x)f (x)g(x)dx (1.15)
a
Les fonctions sont dites orthogonales si leur produit scalaire est nul. La
fonc- tion est normalis´ee quand < f f >= 1. Un ensemble orthonorm´e de
fonctions est un ensemble de |fonctions toutes normalis´ees et orthogonales
deux `a deux.
On peut construire de mani`ere syst´ematique une suite de polynˆomes
telle que le coefficient du monoˆme de degr´e le plus ´elev´e soit ´egal `a un
et telle qu’ils soient tous orthogonaux.
La relation de r´ecurrence est la suivante.

p−1(x) = 0 (1.16)
p0(x) = 1 (1.17)
pi+1 (x) = (x − ai )pi (x) − bi pi−1 (x) (1.18)

avec les coefficients ai et bi d´etermin´es de la mani`ere suivante

< x p i |p i >
a i= < pi |pi > (1.19)
< x p i | p i− 1 >
b i =< pi−1 |pi−1 > (1.20)

(1.21)

Si l’on divise chaque polynˆome par < pi pi >1/2 , on obtient alors des
|
polynˆomes normalis´es.
Une propri´et´e fondamentale de ces polynˆomes ainsi construits est la

1
sui- vante : Le polynˆome pi a exactement j racines distinctes plac´ees sur
l’intervalle [a, b]. Chaque racine du polynˆome pi se trouve entre deux
racines cons´ecutives du polynˆome pi+1 .

1
Int´egration et sommes discr`etes

Les N racines du polynˆome pN sont choisies comme abscisses dans l’´evalua-

tion de l’int´egrale de Gauss et les poids wi sont calcul´es `a partir de la
formule
< pN −1| 1 |pN −1 >
wi = x−xi (1.22)
pN−1(xi)pJN (xi)
dN < pN −1 |pN −1 > (1.23)
= J
dN −1 pN −1 (xi )p N (xi)
ou` le symbole prime designe´ la derivee´´ du polynˆome et dn le coefficient du
monoˆme le plus ´elev´e du polynˆome pN 3 .
Avec ce choix, on peut montrer en utilisant la relation de r´ecurrence, Eq.
(1.18), que < pi p1 >= 0.
|
A titre d’exemples voici les fonctions de poids classiques que l’on utilise pour
calculer une int´egrale par la m´ethode de Gauss
– Gauss-Legendre
La fonction de poids est W = 1 sur l’intervalle [ 1, 1]. La relation de r
−
´ecurrence pour les polynˆomes de Legendre est

(i + 1)Pi+1 = (2i + 1)xPi − iPi−1 (1.24)

– Gauss-Hermite
La fonction de poids est W = exp( x2 ) sur la droite r´eelle. La
−
relation de r´ecurrence pour les polynˆomes d’Hermite est

Hi+1 = 2xTi − 2iTi−1 (1.25)

– Gauss-Laguerre
La fonction de poids est W = xαe−x sur l’intervalle [0, + ∞[. La relation
de r´ecurrence pour les polynˆomes de Laguerre est
(i + 1)Lα = (−x + 2i + α + 1)Lα − (i + (1.26)
α)Lα
i+1 i i−1

– Gauss-Jacobi
La fonction de poids est W = (1 —x)α(1 + x)β sur l’intervalle ] 1, 1[.
La relation de r´ecurrence pour les polynˆomes de Jacobi est

(α,β) (α,β)
c(α,β)
iP = (di + eix)P − f iP (1.27)
i+1 i i−1

ou` les coefficients ci , di , ei et fi sont donn´es par les relations

ci = 2(i + 1)(i + α + β + 1)(2i + α + β) (1.28)

di = (2i + α + β + 1)(α − β ) 2 2
(1.29)
ei = (2i + α + β)(2i + α + β + 1)(2i + α + β + 2) (1.30)
fi = 2(i + α)(i + β)(2i + α + β + 2) (1.31)

On peut utiliser des fonctions de poids qui sont int´egrables sur l´intervalle,
sans ˆetre n´ecessairement born´ees.

1
 3
Attention, la formule 4.5.9 de la r´ef´erence[1] est incorrecte et doit ˆetre remlac´ee par
l’´equa- tion (1.22)

1
 1.5 Int´egrales multiples

−
– Gauss-Chebyshev4 La fonction de poids est W = (1 − x2) 1/2 sur l’inter-
valle [−1, 1]. La relation de r´ecurrence pour les polynˆomes de
Chebyshev.

Ti+1 = 2xTi − Ti−1 (1.32)

1.5 Int´egrales multiples

Le probl`eme de l’´evaluation des int´egrales multiples est ´etroitement li
´e a` la difficult´e de l’´evaluation num´erique d’un tr`es grand nombre de
points pour la fonction consid´er´ee. Par exemple dans un espace `a trois
dimensions, si on utilise 30 points dans chacune des directions, il est n
´ecessaire de calculer la fonction pour 303 points. La situation s’aggrave tr`es
rapidement quand la dimension de l’espace augmente ! Le calcul des int
´egrales multiples est n´eanmoins possible dans un certain nombre de cas. Si
la fonction poss`ede une sym´etrie importante, par exemple la sym´etrie sph
´erique dans un espace de dimension d, on peut se ramener de mani`ere
analytique `a une int´egrale `a une dimension (Voir appendice A). De mani`ere
g´en´erale, en utilisant une sym´etrie de la fonction, on peut rame- ner le calcul
de l’int´egrale de dimension n `a celui d´une int´egrale de dimension nJ << n
ou` les m´ethodes pr´ec´edentes peuvent encore s’appliquer.
Dans le cas ou` il n’existe pas de sym´etrie, la m´ethode la plus efficace
est la m´ethode dite de Monte Carlo dont le principe est d´ecrit dans le cours
Simulation num´erique en Physique Statistique.(http ://[Link])

4
Pafnuty Lvovich Chebyshev (1821-1894) a un nom dont l´orthographe varie un peu selon
les langues, puisqu´il s’agit d´une traduction phon´etique. En Franc¸ais, son nom est g´en

1
´erale- ment orthographi´e Tchebytchev.

1
Int´egration et sommes discr`etes

1
Chapitre 2

Fonctions sp´eciales et
´evaluation de fonctions

Contenu
2.1 Introduction....................................................................... 11
2.2 Fonctions transcendantes simples..................................12
2.3 Fonction Gamma................................................................12
2.3.1 D´efinition et propri´et´es ...................................................12
2.3.2 Fonctions reli´ees : Ψ, B ......................................................14
2.4 Fonctions de Bessel...........................................................15
2.5 Fonctions Hyperg´eom´etriques ...........................................17
2.5.1 Fonction Hyperg´eom´etrique Gaussienne ........................17
2.5.2 Fonctions Hyperg´eom´etriques g´en´eralis´ees ...............18
2.6 Fonction erreur, exponentielle int´egrale .........................18
2.7 Conclusion.......................................................................... 20

2.1 Introduction
Les fonctions sp´eciales sont d´efinies de mani`ere assez impr´ecise,
puisqu’elles regroupent les fonctions que l’usage (ou la fr´equence d’utilisation)
a fini par as- socier `a un nom. Parmi ces fonctions, on trouve un grand
nombre de fonctions qui sont des solutions d’´equations diff´erentielles du
second ordre, sans que cette propri´et´e soit exclusive. Ces fonctions sont
toutefois tr`es utiles, car elles appa- raissent tr`es souvent, d`es que l’on
cherche `a r´esoudre des ´equations diff´erentielles du second ordre dont les
coefficients ne sont pas constants. Les fonctions sp´e- ciales sont disponibles
en programmation sous la forme de biblioth`eques. Elles sont aussi d´efinies
pour un grand nombre d’entre elles dans les logiciels de cal- cul symbolique
(Maple, Mathematica,...). Dans la suite de ce cours, nous allons d´efinir une
partie d’entre elles et d´ecrire les m´ethodes num´eriques utilis´ees dans les
biblioth`eques de programmes pour le calcul de ces fonctions.

1
Fonctions sp´eciales et ´evaluation de fonctions

2.2 Fonctions transcendantes simples

Les fonctions les plus simples que l’on rencontre lors de l’apprentissage des
math´ematiques sont tout d’abord les monoˆmes, puis les polynˆomes, et enfin
les fractions rationnelles. Le calcul de la valeur de la fonction pour un argument
r´eel ou complexe n´ecessite un nombre fini des quatre op´erations ´el
´ementaires que sont l’addition, la soustraction, la multiplication et la
division.
Les premi`eres fonctions transcendantes que l’on rencontre sont alors les
fonc- tions trigonom´etriques (sin, cos,tan, arccos, arcsin) ainsi que leurs
fonctions in- verses. Les fonctions exp et log repr´esentent g´en´eralement le
reste de l’arsenal des fonctions transcendantes d´efinies dans un cursus de
premier, voire de second cycle universitaire.

2.3 Fonction Gamma

2.3.1 D´efinition et propri´et´es
La fonction Gamma est g´en´eralement d´efinie par l’int´egrale suivante
∞
∫
Γ(z) = tz−1e−tdt (2.1)
0

quand la partie r´eelle de z est strictement positive, Re(z) > 0).

La formule d’Euler donne une expression de la fonction Γ pour toute valeur
de z complexe hormis les valeurs de z enti`eres n´egatives ou` la fonction
poss`ede des poˆles :
Γ(z) = lim n!nz (2.2)
n→∞ z(z + 1) . . . (z + n)
En int´egrant par parties l’´equation (2.1), on peut facilement montrer que

Γ(z + 1) = zΓ(z) (2.3)

En v´erifiant que Γ(1) = 1, on obtient par r´ecurrence que

Γ(n + 1) = n! (2.4)

Avec cette d´efinition, la fonction Γ apparaˆıt comme un prolongement

analytique de la fonction factorielle d´efinie sur N. D’apr`es l’´equation (2.2),
la fonction Γ a un poˆle en 0 et pour toutes les valeurs enti`eres n´egatives.
La formule suivante permet de relier la fonction entre les valeurs situ
´ees dans le demi-plan complexe ou` Re(z) > 1 et celui ou` Re(z) < 1 :
π
Γ(1 )=
(2.5)
—
Γ(z) sin(πz)
Pour calculer num´eriquement la fonction Γ pour une valeur de z en dehors
des poˆles, il est n´ecessaire de d´evelopper cette fonction sur la base des

1
polynˆomes et des exponentielles. La formule la plus pr´ecise est celle de
Lanczo`s. Ce d´evelop- pement est sp´ecifique `a la fonction Γ. La formule
qui s’inspire de la formule de

1
 2.3 Fonction Gamma

20

15
ln(

10

0
0 2 4 6 8 10
x

FIG. 2.1 – La fonction ln(Γ(x)) en fonction de x.

Stirling bien connue pour la fonction factorielle n’est valable que pour Re(z) > 0
et est donn´ee par
1
z+
z+γ+
Γ(z + 1) 2
e−(z+γ+ 2 )
1
1
=
2
√ c1 c2 cN
× 2π c 0 + + +...+ +ϵ (2.6)
z+ z+ z+
ou` ϵ est le param`etre estimant l’erreur. Pour le choix particulier γ = 5,
−10
N = 6 et c0 tr`es voisin de 1, | on a ϵ < 2.10 .
Il est difficile de calculer la fonction Γ pour des valeurs de z un peu im-
portantes. Cela r´esulte de la croissance tr`es rapide de la fonction Γ. On
peut montrer que la fonction Γ croˆıt plus vite que toute exponentielle, comme de
ma- ni`ere analogue on montre que l’exponentielle croˆıt plus vite que tout
polynˆome. On dit parfois que la fonction Γ a une croissance super-
exponentielle.
Dans de nombreuses formules, la fonction Γ apparaˆıt `a la fois au num
´erateur et au d´enominateur d’une expression. Chacun des termes peut ˆetre
tr`es impor- tant, mais le rapport est souvent un nombre relativement modeste.
Pour calculer num´eriquement ce type d’expression, il est pr´ef´erable de
calculer ln(Γ(z)) (voir Fig. 2.1). Ainsi la fraction est alors l’exponentielle de
la diff´erence de deux logarithmes. Tous les nombres qui interviennent dans

2
ce calcul sont exponen- tiellement plus petits que ceux qui apparaissent dans
un calcul direct, on ´evite ainsi le d´epassement de capacit´e de l’ordinateur.

2
Fonctions sp´eciales et ´evaluation de fonctions

0.8 a=1
a=3

a=10
0.6
P(a

0.4

0.2

0
0 2 4 8 10 12 14
6 x

FIG. 2.2 – La fonction P (a, x) en fonction de x pour trois valeurs de a (a =

1, 3, 10).

On d´efinit la fonction γ incompl`ete1 comme

∫x
γ(a, x) = −t
ta−1 e dt (2.7)
0

La fonction normalis´ee suivante P (a,

x)

γ(a, x)
P (a, x) = Γ(a) (2.8)

est parfois appel´ee aussi fonction Gamma incompl`ete. On peut montrer

que P (a, x) est monotone croissante avec x. La fonction est tr`es proche de 0
quand x est inf´erieur `a a 1 et proche de 1 quand x est tr`es sup
−
´erieur. La variation entre ces deux valeurs apparaˆıt autour de l’abscisse 1 et
a
sur une largeur de
√
l’ordre de a (voir figure 2.2).

2.3.2 Fonctions reli´ees : Ψ, B

A partir de la fonction Γ, on d´efinit des fonctions d´eriv´ees. En raison
de leur grande fr´equence d’utilisation, elles ont “re¸cu” un nom. Ainsi, la
fonction
2
 1
Noter que dans le logiciel Maple, la fonction Gamma et les fonctions Gamma incompl`etes
sont les seules fonctions d´efinies avec des lettres capitales.

2
 2.4 Fonctions de Bessel

Ψ, appel´ee aussi fonction Digamma est d´efinie comme la d´eriv´ee

logarithmique de la fonction Gamma2 :
d ln(Γ(x)) (2.10)
Ψ(x) = dx
Parmi les propri´et´es remarquables de la fonction Ψ, notons que, pour
des valeurs enti`eres, on a
n−1
Ψ(n) = Σ
γ+ (2.11)
1 i
i=1
−
ou` γ = 0.577 . . . est la constante d´Euler.
Les fonctions Beta qui sont not´ees paradoxalement avec un B sont d
´efinies par la relation :
B(z, w) = Γ(z)Γ(w)
Γ(z + w) (2.12)

2.4 Fonctions de Bessel

Les fonctions de Bessel sont d´efinies de la mani`ere suivante : consid
´erons l’´equation diff´erentielle du second ordre

x2yJJ + xyJ + (x2 − ν2)y = 0 (2.13)

Les solutions de cette ´equation sont appel´ees fonctions de Bessel de

premi`ere et de deuxi`eme esp`ece : La solution finie `a l’origine et not´ee
Jν (x) est appel´ee fonction de Bessel de premi`ere esp`ece et la seconde
solution not´ee Yν (x) est appel´ee fonction de Bessel de deuxi`eme esp`ece.
Si ν n’est pas un entier, ces fonctions sont reli´ees par la relation suivante :

Jν(x) cos(νπ) − J−ν(x)

Y (x) = sin(νπ) (2.14)
ν
La figure 2.3 repr´esente graphiquement les fonctions de Bessel de
premi`ere et de seconde esp`eces pour les quatre premi`eres valeurs enti`eres
de ν
Le comportement asymptotique des fonctions de Bessel de premi`ere et
de seconde esp`ece est le suivant
r
2
J (x) ' (cos(x − νπ/2 − π/4)) (2.15)
ν r π
2
Y (x) ' (sin(x − νπ/2 − π/4)) (2.16)
ν
π
Soit l’´equation diff´erentielle du second ordre

x2yJJ + xyJ − (x2 − ν2)y = 0 (2.17)

2
On d´efinit aussi les fonctions polygamma comme une g´en´eralisation de la fonction
Di- gamma
n
Ψ(n, x) = d Ψ(x)
(2.9)
dxn

2
Fonctions sp´eciales et ´evaluation de fonctions

0.5

Y(x),J 0

-0.5
J0 J1 J2 J3 Y0 Y1 Y2 Y3

-1

-1.5

-2
0 2 4
x
6 8 10

FIG. 2.3 – Les quatre premi`eres fonctions de Bessel enti`eres de premi`ere

es- p`ece et de deuxi`eme esp`ece. Ces fonctions sont pr´esentes dans toutes
les biblio- th`eques math´ematiques de programmation, dans les logiciels de
calcul symbo- lique comme Maple et dans un logiciel graphique comme
xmgrace.

Les solutions de cette ´equation sont appel´ees fonctions de Bessel modifi´ees .

La solution finie `a l’origine et not´ee Iν (x) est appel´ee fonction de Bessel
modifi´ee de premi`ere esp`ece et la seconde solution not´ee Kν (x) est
appel´ee fonction de Bessel modifi´ee de seconde esp`ece. Ces fonctions
sont reli´ees par la relation suivante :
π (Iν (−x) − Iν (x))
K ν(x) = (2.18)
2 sin(νπ)

La figure 2.4 repr´esente graphiquement les fonctions de Bessel modifi´ees

de premi`ere et de deuxi`eme esp`ece pour les quatre premi`eres valeurs
enti`eres de ν. Le comportement asymptotique des fonctions de Bessel modifi
´ees est le sui-
vant :

ez
Iν(x) ' √ (2.19)
r 2πx
π −z
ν '
K (x) e (2.20)
2

2
Les fonctions de Hankel H1,ν and H2,ν sont appel´ees fonctions de Bessel de

2
 2.5 Fonctions Hyperg´eom´etriques

4
I0 I1 I2 I3 K0 K1 K2 K3

3.5

2.5
K(x),I

1.5

0.5

0
0 2 4
x

FIG. 2.4 – Les quatre premi`eres fonctions de Bessel enti`eres modifi´ees

de pre- mi`ere et de deuxi`eme esp`ece. Ces fonctions sont pr´esentes dans
toutes les bi- blioth`eques math´ematiques de programmation, dans les
logiciels de calcul sym- bolique comme Maple et dans un logiciel graphique
comme xmgrace.

troisi`eme esp`ece et sont d´efinies par la relation

H1,ν(x) =Jν(x) + iYν(x) (2.21)

H2,ν(x) =Jν(x) − iYν(x) (2.22)

2.5 Fonctions Hyperg´eom´etriques

2.5.1 Fonction Hyperg´eom´etrique Gaussienne
Les fonctions hyperg´eom´etriques gaussiennes sont d´efinies comme
´etant les solutions de l’´equation diff´erentielle suivante.

x(1 − x)yJJ + [c − (a + b + 1)x]yJ − aby = 0 (2.23)

ou` a, b et c sont des constantes.

Si c, a − b et c − a − b sont non entiers, la solution g´en´erale de cette ´equation
est
y = F (a, b; c; x) + Bx1−cF (a − c + 1, b − c + 1; 2 − c; x) (2.24)
2
Fonctions sp´eciales et ´evaluation de fonctions

La fonction F peut ˆetre exprim´ee sous la forme d’une s´erie

F (a, b; c; z) ≡ 2F1(a, b, c;
z)
(2.25)
Γ(c) Σ∞ Γ(a + n)Γ(b + n)
=
zn Γ(a)Γ( Γ(c + n
n
Cette s´erie converge uniform´ement `a l’int´erieur du disque unit´e. D`es
que a, b ou c sont entiers, la fonction hyperg´eom´etrique peut se r´eduire
`a une fonction transcendante plus simple. Par exemple, 2F1(1, 1, 2; z) = −z−1
ln(1 − z)

2.5.2 Fonctions Hyperg´eom´etriques g´en´eralis´ees

On d´efinit des fonctions hyperg´eom´etriques g´en´eralis´ees de la
mani`ere sui- vante : soit le rapport

∞
Σ (a1)k(a2)k
p
a1, a2, . . . , ap ; =k=0 1 k 2 k . . . (ap)k
q k xk (2.2
( ( ... 6)
ou` on a utilis´e la notation de
Pochhammer

(a) k = Γ(a + k) (2.27)

Γ(a)

2.6 Fonction erreur, exponentielle int´egrale

La fonction erreur et la fonction erreur compl´ementaire sont d´efinies comme
2
erf (x) = √ ∫ −t2 dt
x

π e (2.28)
0
erfc(x) = 1 − erf (x)
2 ∫∞
=√ e−t dt
2
(2.29)
π x
La fonction erf est aussi pr´esente dans toutes les biblioth`eques standard de
programmation.
La fonction exponentielle int´egrale Ei est d´efinie comme la valeur
principale de l´int´egrale suivante pour x > 0.

Ei(x) = − ∫ ∞ t
− dt (2.30)
−x t
∫x
et (2.31)
= −∞ t

2
Le d´eveloppement en s´erie de cette fonction donne
∞ n
Σx
Ei(x) = γ + ln(x) + (2.32)
n=1 nn!

2
 2.6 Fonction erreur, exponentielle int´egrale

E1
E2
E3
3 E4

2
En

0
0 1 2 3 4
x

FIG. 2.5 – Les quatre premi`eres fonctions exponentielles int´egrales

(fonctions En . Ces fonctions sont pr´esentes dans toutes les biblioth`eques
math´ematiques de programmation, dans les logiciels de calcul symbolique
comme Maple et dans un logiciel graphique comme xmgrace.

Pour des grandes valeurs de x, on a le d´eveloppement asymptotique suivant

ex 1 Ei(x 1 +. (2.3
..
x
De mani`ere g´en´erale, on d´efinit les exponentielles int´egrales En (x) comme
∫ ∞ −zt
e
En(z) dt
La Figure 2.5 repr´esente les quatre premi`eres exponentielles int´egrales.
1
Le d´e- veloppement en s´erie de cette fonction donne
∞ n
Σ x
E 1(x) = −(γ + ln(x)) + (−1)n (2.35)
n
n
La fonction Ei (1, x) n´est d´efinie que pour des arguments r´eels : Pour x <
0, on a
Ei (x) = −Ei (1, −x) (2.36)
On peut noter que les exponentielles int´egrales En (x) sont reli´ees `a la
fonc- tion γ par la relation
En(x) = xn−1γ(1 − n, x) (2.37)

3
Fonctions sp´eciales et ´evaluation de fonctions

2.7 Conclusion
Cette introduction aux fonctions sp´eciales est tr`es loin d’ˆetre exhaustive
; il existe de nombreuses autres fonctions dites sp´eciales : les fonctions
elliptiques, les fonctions de Fresnel, les fonctions de Meier,. Le d
´eveloppement de biblio-
th`eques qui permettent de calculer les valeurs de ces fonctions est un
secteur tr`es actif et nous disposerons dans les ann´ees futures de
biblioth`eques encore plus performantes.

3
Chapitre 3

Racines d’´equations

Contenu
3.1 Introduction....................................................................... 21
3.2 Dichotomie.........................................................................22
3.3 M´ethode de Ridder .................................................................23
3.3.1 M´ethode de la position fausse ...........................................23
3.3.2 M´ethode de Ridder .............................................................24
3.4 M´ethode de Brent ...................................................................25
3.5 Newton-Raphson...............................................................25
3.6 Racines de Polynˆomes ............................................................26
3.6.1 R´eduction polynomiale .......................................................26
3.6.2 M´ethode de Laguerre ..........................................................27

3.1 Introduction
L’une des taˆches rencontr´ees fr´equemment lors d’un calcul est la
recherche de la racine d’une ´equation. Sans perte de g´en´eralit´e, on peut
toujours ´ecrire une ´equation ou` le membre de droite est ´egal `a z´ero,

f (x) = 0 (3.1)

Si x est une variable scalaire, le probl`eme est unidimensionnel. Si x est

une variable vectorielle (`a N dimensions) et que l’on a N ´equations `a
satisfaire, on peut formellement ´ecrire sous une notation vectorielle

f(x) = 0 (3.2)

Malgr´e la similarit´e des ´equations (3.1) et (3.2), un syst`eme

d’´equations `a N variables est consid´erablement plus compliqu´e `a r
´esoudre qu’un syst`eme unidi- mensionnel. La raison vient du fait que la m
´ethode g´en´erale pour la recherche de racines est li´ee `a la capacit´e
d’encadrer num´eriquement la r´egion ou` le syst`eme d’´equations poss`ede
une racine particuli`ere.
On exclut de ce chapitre le cas des syst`emes lin´eaires qui sera trait´e
dans le chapitre de l’alg`ebre lin´eaire. Le principe dominant la recherche

3
de racines

3
Racines d’´equations

d’´equations est celui de m´ethodes it´eratives, ou` en partant d’une valeur

d’essai (ou un couple de valeurs d’essai), on s’approche de plus en plus
pr`es de la solution exacte. Il est ´evident qu’une estimation de d´epart
raisonnable associ´ee
`a une fonction f qui varie suffisamment lentement est n´ecessaire pour
obtenir une convergence vers la solution recherch´ee.
Nous allons consid´erer le probl`eme unidimensionnel pour lequel
plusieurs m´ethodes sont disponibles afin de choisir la m´ethode la mieux
adapt´ee compte tenu des informations que l’on dispose sur la fonction f . Une
derni`ere partie de ce chapitre sera consacr´ee aux m´ethodes plus sp
´ecifiques pour la recherche de racines de polynˆomes

3.2 Dichotomie
Comme nous l’avons mentionn´e ci-dessus, la cl´e de la recherche de
racines d’´equations repose sur l’existence d’un encadrement pr´ealable de
cette racine. S’il existe un couple (a, b) tel que le produit f (a)f (b) < 0 et si la
fonction est continue, le th´eor`eme de la valeur interm´ediaire nous dit que
fonction s’annule au moins une fois `a l’int´erieur de cet intervalle.
La m´ethode de dichotomie est une m´ethode qui ne peut pas ´echouer,
mais sa rapidit´e de convergence n’est pas la meilleure en comparaison avec les
autres m´ethodes. L’id´ee de cette m´ethode est la suivante : soit une fonction f
monotone sur un intervalle [a0, b0] telle que f (a0)f (b0) < 0 , on sait alors qu’il
existe une et une seule racine comprise dans cet intervalle.
L´algorithme de la m´ethode de dichotomie est le suivante : tout
d’abord, on calcule f ( a0+b0 ).
2
– Si f ( a0 +b0 )f (a0 ) < 0, on d´efinit un nouvel encadrement de la racine
2
par le couple (a1, b1) tel que
a1 = a0 (3.3)
a0 + b0
b1 = . (3.4)
2
– Si f ( a0 +b0 )f (a0 ) > 0, alors on d´efinit un nouvel encadrement de la racine
2
par le couple (a1, b1) tel que
a0 + b0
a1 = (3.5)
2
b1 = b0. (3.6)
En it´erant cette m´ethode, on obtient une suite de couple (an , bn ) telle que
ϵn = bn − an v´erifie la relation
ϵn+1 ϵn
= 2 (3.7)

ou` ϵ0 = (b0 a0 )/2 Cela signifie que si l’on se fixe la tol´erance ϵ qui repr
−
´esente la pr´ecision `a laquelle on souhaite obtenir la racine, on a un nombre
d’it´erations
`a effectuer ´egal `a
n = ln2
3
ϵ
|b0 − a 0 |
(3.8)

3
 3.3 M´ethode de Ridder

ou` la notation ln2 signifie le logarithme en base 2.

Le choix de la valeur de la tol´erance n´ecessite quelques pr´ecautions.
Si la racine recherch´ee est de l’ordre de l’unit´e, on peut tr`es raisonnablement
− −
choisir ϵ 0 de l’ordre de 10 6 `a 10 13 selon que l’on travaille en simple ou
double pr´ecision.
Par contre pour une racine dont la valeur est de l’ordre de 1010 , une pr
−
´ecision de 10 4 sera la valeur maximale que l’on peut atteindre en double
pr´ecision. Inversement, pour une racine proche de z´ero, la pr´ecision peut
−
ˆetre meilleure que 10 14.

3.3 M´ethode de Ridder

3.3.1 M´ethode de la position fausse
La m´ethode de dichotomie sous-exploite le fait que l’on peut mieux utiliser
la fonction `a l’int´erieur de l’encadrement effectu´e `a chaque it´eration. La
m´ethode de la position fausse approche la fonction de mani`ere lin´eaire
dans l’intervalle consid´er´e.
Soit la droite y = cx+d passant par f (a0) et f (b0) en a0 et b0 respectivement,
on obtient facilement que
f (b0) − f (a0)
c= b0 − a0 (3.9)
b0f (a0) − a0f (b0)
d= b0 − a0 (3.10)
La nouvelle abscisse estim´ee pour la racine de l’´equation est donn´ee par y =
cx + d = 0, ce qui donne

d
x=− (3.11)
c
b0f (a0) − a0f (b0)
= f (b ) − f (a ) (3.12)
0 0

soit encore
f (a0)
x = a0 − (b0 − a0) (3.13)
f (b )0 − f (a 0)
= b + (b f (b0)
—a) (3.14)
0 0 0
f (b0) − f (a0)
On reprend `a ce stade le principe de l’algorithme pr´ec´edent si f (x)f (a0 ) >
0 alors
La figure

sinon

3
a1 = x

b1 = b0

a1 = a0

b1 = x

3
Racines d’´equations

2
4

4
0

1
-2

0 0.5 1 1.5 2 2.5 3

FIG. 3.1 – Sch´ema illustrant le principe de la position fausse. Les lignes en

pointill´e correspondent aux interpolations lin´eaires. A partir de
l’encadrement rep´er´e par les points 1 et 2, le processus it´eratif conduit aux
points 3, puis 4...

3.3.2 M´ethode de Ridder

Une variante de la m´ethode pr´ec´edente qui est tr`es efficace est bas´ee sur
l’algorithme suivant. On ´evalue la fonction au point x = a0 +b0 et on resout´
2
l’´equation en
z f (a0) − 2f (x)z + f (b0)z2 = 0 (3.19)
La solution positive est donn´ee par
( )+ ( ( )) ( )2 ( ) ( )
z= √
f x sgn f xf (b0) f x −f (3.20)
a f b
En appliquant la m´ethode de0 0 la position fausse non pas `a f (a0 ), f (x) et f
(b0 ), mais `a f (a0 ), f (x)z et f (b0 )z 2 , on obtient une approximation de la
racine, not´ee x4 et donn´ee par

x4 = x + (x −
sgn(f (a0) − f (b0))f (3.21)
a ) (x)
0 √
f (x)2 − f (a )f
0
(b )
0
Parmi les propri´et´es remarquables, notons que x4 est toujours situ´ee `a
l’in- t´erieur de l’intervalle [a0 , b0 ]. Si le produit f (x3 )f (x4 ) est n´egatif,
on prend l´intervalle [x3 , x4 ] comme nouvelle encadrement, sinon si on

3
consid`ere le pro- duit f (x1 )f (x4 ) ; si celui est n´egatif, le nouvel
encadrement est [x1 , x4 ], sinon on prend [x4 , x2 ]. On it`ere ensuite le proc´ed
´e.

3
 3.4 M´ethode de Brent

3.4 M´ethode de Brent

Le principe de cette m´ethode est de combiner les avantages des m
´ethodes pr´ec´edemment expos´ees en utilisant, le principe de l’encadrement
de la racine, la dichotomie, et l’interpolation quadratique inverse. Cela n
´ecessite de connaˆıtre trois valeurs de la fonction f dont la racine est `a d
´eterminer. Soit (a, f (a)), (b, f (b)), et (c, f (c)) la formule d’interpolation est
donn´ee par

(y − f (a))(y − f (b))c
x = (f (c) − f (a))(f (c) − f (b))
(y − f (b))(y − f (c))a (y − f (c))(y − f (a))b (3.22)
+ +
(f (a) − f (b))(f (a) − f (c)) (f (b) − f (c))(f (b) − f (a))
En choisissant y = 0, on peut ´ecrire l’´equation (3.22)
comme

x=b+ P (3.23)
Q
ou` P et Q sont donn´es par

P = S[T (R − T )(c − b) − (1 − R)(b − a)] (3.24)

Q = (T − 1)(R − 1)(S − 1) (3.25)

ou` R, S et T s’expriment comme

f (b)
R=
f (c) (3.26)
f (b)
S=
f (a) (3.27)
f (a)
T =
f (c) (3.28)

En pratique, b est une premi`ere estimation de la racine et P une petite correc-

Q
tion. Quand Q → 0 la valeur de P peut devenir tr`es grande et l’iteration´ par la
Q
m´ethode de Brent est remplac´ee par une it´eration de dichotomie.

3.5 Newton-Raphson
Toutes les m´ethodes pr´ec´edentes ne n´ecessitaient que la
connaissance de la fonction en diff´erents points de l’intervalle encadrant la
racine. Sous r´eserve que la variation de la fonction ne soit pas trop rapide,
seule une hypoth`ese de continuit´e est n´ecessaire.
La m´ethode de Newton-Raphson n´ecessite de plus que la fonction f
dont on cherche `a d´eterminer une racine, soit d´erivable au voisinage de
celle-ci.
Les it´erations successives de la m´ethode de Newton-Raphson sont bas´ees
sur le d´eveloppement limit´e de la fonction autour d’un point

4
 JJ
f (x)
f (x + δ) = f (x) + f (x)δ +
J
δ2 + . . . (3.29)
2

4
Racines d’´equations

Si δ est suffisamment petit, on peut n´egliger les termes non lin´eaires et

une estimation de la racine est donn´ee par f (x + δ) = 0.
f (x)
δ=− (3.30)
f J(x)
On voit imm´ediatement qu’il est n´ecessaire que la d´eriv´ee de la
fonction ne s’annule pas dans le voisinage de x, sous peine que l’estimation de δ
devienne tr`es grande et ne permette pas `a la m´ethode de converger.
Si les conditions pr´ec´edemment ´enonc´ees sont v´erifi´ees, on a une
m´ethode qui converge de mani`ere quadratique.
En effet, la relation de r´ecurrence entre estimations successives est donn´ee
par
f (xi)
xi+1 = xi − (3.31)
f J(x )i
En posant ϵi+1 = xi+1 − x, ou` x est la racine exacte, on a
f (xi)
ϵi+1 = ϵi − (3.32)
f J(x )i
Si on utilise un d´eveloppement limit´e de f au deuxi`eme ordre au point
xi (ce qui suppose que la fonction est deux fois d´erivable au voisinage de
la racine), on obtient
ϵi+1 f JJ(xi)
= − 2i 2f (3.33)
J
(xi)
La m´ethode converge donc tr`es rapidement par comparaison avec les m
´ethodes pr´ec´edentes.
A noter que si la d´eriv´ee de la fonction n’est pas connue
analytiquement, l’´evaluation num´erique de sa d´eriv´ee est possible par une
formule d’accroisse-
ment
f (x + ∆x) − f (3.34)
f(xJ(x)
) ' ∆

Dans ce cas, la m´ethode de Newton-Raphson se r´eduit `a une m´ethode

d’in- tersection et la convergence de celle-ci est moins rapide que la
convergence quadratique.

3.6 Racines de Polynoˆmes

3.6.1 R´eduction polynomiale
La recherche de racines d’un polynˆome se construit de la mani`ere suivante :
soit Pn (x) un polynˆome de degr´e n. Si on obtient une premi`ere racine, on
peut
´ecrire
Pn(x) = (x − x1)Pn−1(x) (3.35)
ou` Pn−1 (x) est un polynˆome de degr´e n 1. Ainsi, th´eoriquement, une fois
−
4
obte- nue une premi`ere racine, on peut recommencer la recherche d’une
autre racine pour une polynˆome de degr´e strictement inf´erieur.
Successivement, on poursuit cette proc´edure jusqu’`a l’obtention de
l’ensemble des n racines du polynˆome

4
 3.6 Racines de Polynˆomes

Pn (x). Rappelons que les polynˆomes `a coefficients complexes se

factorisent en un produit de monoˆmes de degr´e 1. Cette propri´et´e
exprime le fait que les po- lynoˆmes `a coefficients complexes ont l’ensemble
de leurs racines dans le plan complexe,
n
Y
Pn(x) = (x − xi) (3.36)
i=1
(C est un corps alg´ebriquement clos).

3.6.2 M´ethode de Laguerre

Les m´ethodes de recherche de z´eros de polynˆomes sont nombreuses
et une pr´esentation d´etaill´ee d´epasse largement le cadre de ce cours. Nous
avons choisi de pr´esenter une m´ethode dont le principe est assez simple.
La m´ethode de Laguerre utilise le fait que les d´eriv´ees logarithmiques
successives d’un polynˆome divergent au voisinage d’une racine. En prenant
le logarithme de l’´equation (3.36), on obtient
n
Σ
ln(|Pn (x)|) = ln(|x − xi |) (3.37)
i=1
En d´erivant l’´equation (3.37), on obtient
d ln( |Pn(x)|) Σ 1
n
dx = x − xi
i=1
J
Pn (x)
= Pn(x) (3.38)
=G (3.39)
En d´erivant l’´equation (3.38), il
vient n
2

d ln(|Pn(x)|) Σ 1
− d = (x −
i 2 JJ
J
nP Pn (x
= −)
Pn(x)
(x)
Pn(x)
=H (3.40)
Soit la racine x1 `a d´eterminer, on suppose que la valeur de d´epart x est
situ´ee
`a une distance a de x1 et que l’ensemble des autres racines sont situ´ees
`a une distance suppos´ee identique et qui vaut b
x − x1 = a (3.41)
x − xi = b i ∈ [2, n] (3.42)
En ins´erant les ´equations (3.41) et (3.42) dans les ´equations (3.38),
(3.40), on en d´eduit respectivement les relations suivantes
1 n− 1
+ =G (3.43)
a b
1 n− 1
+ =H (3.44)
4
a2 b2

4
Racines d’´equations

Apr`es ´elimination de b, la valeur de a

est (3.45)
n
=
a √
G ± (n − 1)(nH − G2)
Le signe plac´e devant la racine du d´enominateur est choisi tel que le d
´enomina- teur soit le plus grand possible. x a devient alors la nouvelle valeur
−
de d´epart et on it`ere le processus. En combinant cette m´ethode avec celle de
la r´eduction po- lynomiale, on peut calculer l’ensemble des racines. En
pratique, comme chaque racine n’est d´etermin´ee qu’avec une pr´ecision
finie, il est n´ecessaire d’ajouter une proc´edure dite de lissage pour ´eviter
les probl`emes d’instabilit´e num´erique.

4
Chapitre 4

Equations diff´erentielles

Contenu
4.1 Introduction....................................................................... 29
4.2 D´efinitions ...................................................................................30
4.3 M´ethodes d’int´egration `a pas s´epar´e ......................31
4.3.1 Introduction............................................................................31
4.3.2 M´ethode d’Euler ..................................................................32
4.3.3 M´ethode RK explicites `a un point ..................................32
4.3.4 M´ethodes RK implicites `a un point ................................33
4.3.5 M´ethodes RK explicites `a 2 points interm´ediaires .....33
4.3.6 M´ethodes RK explicites `a 3 points interm´ediaires .....33
4.3.7 Formule g´en´erale des m´ethodes RK explicites ...........34
4.4 M´ethode d’int´egration `a pas variables .......................34
4.4.1 Introduction............................................................................34
4.5 M´ethodes de Runge-Kutta “embarqu´ees” ....................35
4.6 M´ethode de Bulirsh-Stoer ....................................................36
4.7 Conclusion.......................................................................... 37

4.1 Introduction
La r´esolution num´erique d’´equations diff´erentielles est tr`es souvent
n´eces- saire, faute de l’existence de solutions analytiques. Le but de ce chapitre
est de montrer que la meilleure m´ethode, ou la plus efficace `a utiliser pour
obtenir une solution, n´ecessite de connaˆıtre la nature de l’´equation diff
´erentielle `a r´esoudre. Les m´ethodes les plus standards que nous allons pr
´esenter sont largement pr´e- sentes dans les logiciels de calcul comme Maple,
Matlab, Scilab, Octave, ou Mathematica, et surtout dans les biblioth`eques pour
la programmation (IMSL, NAG, GSL). Il est donc pr´ef´erable d’utiliser ces
biblioth`eques plutˆot que de r´e´ecrire un code peu performant et
probablement faux dans un premier temps.

4
Equations diff´erentielles

4.2 D´efinitions
Soit une fonction num´erique not´ee y(x) d´efinie sur un intervalle ferm´e
[a, b] de R et de classe Cp (continuˆment d´erivable d’ordre p). On appelle
´equation diff´erentielle d’ordre p une ´equation de la forme

F (x, y, yJ, yJJ, . . . y(p)) = 0 (4.1)

ou` y J repr´esente la derivee´´ premi`ere par rapport `a x, y JJ

la derivee´´ seconde,
etc... Plus g´en´eralement, on appelle syst`eme diff´erentiel un ensemble
d’´equations diff´erentielles reliant une variable x et un certain nombre de
fonction yi (x) ainsi que leurs d´eriv´ees. L’ordre du syst`eme diff´erentiel
correspond `a l’ordre de d´erivation le plus ´elev´e parmi l’ensemble des
fonsctions.
On appelle solution de l’´equation diff´erentielle (4.1) toute fonction y(x) d
´e- finie sur l’intervalle [a, b] de classe Cp qui v´erifie l’´equation (4.1).
Une classe restreinte d’´equations diff´erentielles ont des solutions math
´ema- tiques “simples”. Les ´equations diff´erentielles `a coefficients constants,
mais aussi (voir chapitre 2) l’´equation hyperg´eom´etrique, les ´equations de
Bessel,...
On appelle forme canonique d’une ´equation diff´erentielle une expression
du type
y(p) = f (x, y, yJ, yJJ, . . . y(p−1)) (4.2)
Seul ce type d’´equations sera consid´er´e dans ce chapitre.
Il est facile de v´erifier que toute ´equation diff´erentielle canonique peut ˆetre
´ecrite comme un syst`eme d’´equations diff´erentielles du premier ordre.
Si on introduit p − 1 fonctions d´efinies comme
y1 = y
y2 = yJ
...
yp = y(p−1) (4.3)

on peut exprimer l’´equation (4.2) sous la forme

y 1J = y 2
y 2J = y
3

...
ypJ = f (x, y, y1, y2, . . . yp) (4.4)

Par souci de commodit´e, on note dans la suite de ce chapitre y une

fonction dont telle que y(x) Rp ou` p est entier positif.
Soit une ´equation ∈ ou syst`eme diff´erentiel (canonique) d´efini sur un
intervalle [a, b] de R

yJ = f (x, y) x ∈ [a, b] y(x) et yJ(x) ∈ Rp (4.5)

4
On appelle condition de Cauchy la condition initiale

y(a) = y0 (4.6)

4
 4.3 M´ethodes d’int´egration `a pas s´epar´e

ou` y0 est un vecteur de Rp .

Le probl`eme de Cauchy est la recherche d’une solution v´erifiant la
condition initiale (4.6).
Pour tout couple de vecteurs y et y∗, et ∀x, la fonction f (x, y) satisfait
la condition de Lipschitz si
ǁf (x, y) − f (x, y ∗ )ǁ ≤ ǁy − y ∗ ǁL (4.7)
ou` L est appel´ee constante de Lipschitz et ou` la notation ǁǁ d´esigne la norme.
Si dans le domaine [a, b] Rp, la fonction f (x, y) est continue par rapport
×
au couple (x, y) et si pour tout x appartenant `a l’intervalle [a, b] elle
satisfait une condition de Lipschitz, alors (x0, y0), il existe un voisinage de ce
∀
point tel que l’´equation (4.5) satisfaisant la condition initiale (4.6) a une
solution unique dans ce voisinage.
La condition de Lipschitz est une condition suffisante de l’existence et de
l’unicit´e de la solution. Dans le cas ou` la fonction f (x, y) est diff´erentiable
born´ee ( ∂f (x, y)/∂y < L), le crit`ere de Cauchy est v´erifi´e.
|| ||
Soit l’´equation diff´erentielle suivante
y
1 + yJ = (4.8)
x
y(0) = 0 (4.9)
On v´erifie ais´ement que la condition de Cauchy n’est pas v´erifi´ee et
dans ce cas, il existe une infinit´e de solutions y = (C ln(x))x ou` C est
−
une constante quelconque qui satisfait la condition initiale.

4.3 M´ethodes d’int´egration `a pas s´epar´e

4.3.1 Introduction
Soit l’´equation diff´erentielle d´efinie par les ´equations (4.5). On suppose
que la fonction f satisfait une condition de Lipschitz afin d’ˆetre certain que la
solution existe, est unique et que le probl`eme est bien pos´e.
On cherche `a calculer une approximation de la solution y(x) en un
certain nombre de points x1 , x2 , . . . xN de l’intervalle [a, b], appel´e maillage
de l’inter- valle, avec x0 = a et xN = b
Nous supposons que la suite des points est choisie de mani`ere a` ce
que la distance entre deux points cons´ecutifs soit constante et on pose
(b − a) (4.10)
h= N
ce qui donne
xk = a + kh (4.11)
avec k = 0, 1, . . . , N
On appelle m´ethode d’int´egration `a pas s´epar´e toute formule de r
´ecurrence de la forme
yk+1 = yk + hφ(xk, yk, h)
k = 0, 1, . . . N avec y0 donn´e (4.12)

5
Equations diff´erentielles

la fonction φ est suppos´ee continue par rapport aux trois variables x, y, h.

On appelle m´ethode `a pas multiples les m´ethodes telles que yk+1 d´epend
de plusieurs valeurs pr´ec´edentes yk , yk−1 , . . . yk−r .

4.3.2 M´ethode d’Euler

Cette m´ethode est d´efinie par
yk+1 = yk + hf (yk, xk) (4.13)
avec y0 donn´e.
Cette m´ethode revient `a approximer la solution au voisinage de xk
par sa tangente et nous allons voir qu’elle est d’ordre 1. En effet, si la
solution est suffisamment d´erivable, on peut ´ecrire

y(xk + h) = y(xk) + hyJ(xk) + h

2
JJ
(xk + θh) (4.14)
avec 0 ≤ θ ≤ 1. ce qui donne 2 y

h2 JJ

+ k + h) = y(xk) + hf (yk, xk)

y(x y (xk + θh) (4.15)
2
D’apr`es la d´efinition de la m´ethode
1 y(xk + h) − y(xk)
(y(xk + h) − y(xk ) − φ(yk (xk ),
h , h)) = h −f (y(x k), xk)
xk
(4.16)
ou
1 yJJ(xk + θh)
(y(xk + h) − y(xk ) − φ(yk (xk ), , h)) = h (4.17)
h xk 2
Si la d´eriv´ee seconde de y est born´ee par une constante K dans
l’intervalle d’int´egration [a, b], on aura
1
max (y(x
k + h) − y(xk ) − φ(yk (xk ), , h))ǁ ≤ Kh (4.18)
ǁh
xk
ce qui montre que la m´ethode est d’ordre un.
La m´ethode d’Euler est une m´ethode num´erique peu couˆteuse num
´erique- ment, mais peu pr´ecise quand on int`egre sur plusieurs pas de
temps. Des am´e- liorations sont possibles d`es que l’on consid`ere des points
interm´ediaires, ce que nous allons voir ci-dessous en consid´erant des m
´ethodes dites de Runge-Kutta

4.3.3 M´ethode RK explicites `a un point

h
yk,1 = yk + f (xk, yk)
2α h
yk+1 = yk + h(1 − α)f (xk, yk) + αf (xk + , yk,1)
2α
5
 y0 donn´e
avec α un nombre r´eel compris entre 0 et 1. Les valeurs de α couramment
utilis´ees sont α = 1, α = 1/2 et α = 3/4. Ces m´ethodes sont d’ordre 2.

5
 4.3 M´ethodes d’int´egration `a pas s´epar´e

4.3.4 M´ethodes RK implicites `a un point

La formule de r´ecurrence est d´efinie par la relation

yk+1 = yk + h[(1 − θ)f (xk, yk) + θf (xk+1, yk+1)] (4.19)

ou` θ est un nombre r´eel appartenant `a l’intervalle ]0, 1] (si θ = 0, on

retrouve la m´ethode d’Euler). Si θ = 1/2, la m´ethode est d’ordre 2 et
s’appelle la m´ethode des trap`ezes. Si θ /= 1/2, la m´ethode est d’ordre 1.

4.3.5 M´ethodes RK explicites `a 2 points interm´ediaires

Ces m´ethodes sont d´efinies par les relations

yk,1 = yk + hf (xk, yk)

3

k 2h k h
y = y k+ 3f x + ,3y k

2h
yk+1 = yk , ) + 3f + ,y (4.20)
h 3
k,2
+ 4 f (x kyk xk
ou par

yk,1 = yk + hf (xk, yk)

2 h
yk,2 = yk , ) + 2f + ,
k,1
yk xk y
+h f 2
(xk
h
yk+1 = yk , ) + 4f + , +f ) (4.21)
h k,1 k+1 ,
+ f (x yk xk y (x yk,2
6 k 2

Ces deux m´ethodes sont d’ordre 3. La premi`ere est parfois appel´ee m

´ethode de Heun.

4.3.6 M´ethodes RK explicites `a 3 points interm´ediaires

La m´ethode suivante est de loin la plus connue et utilis´ee. Les
relations de r´ecurrence sont les suivantes.

yk,1 = yk + hf (xk, yk)

2

k h h
y = y k+ 2f x k+ 2, y k

k h h
y = y k+ 2f x k+ 2, y k

5
 h h
yk+1 = yk , ) + 2f + , + 2f + , +f yk,3)
h k,1 k,2 k+1 ,
+ f (x yk xk y xk y (x
6 k 2 2
(4.22)

Cette m´ethode est d’ordre 4.

5
Equations diff´erentielles

4.3.7 Formule g´en´erale des m´ethodes RK explicites

Les m´ethodes de Runge Kutta s’´ecrivent de mani`ere g´en´erale

K1 = h[f (xk + θ1 h, yk + α1,1 K1 + α1,2 K2 + . . . + α1,n Kn )]

K2 = h[f (xk + θ2 h, yk + α2,1 K1 + α2,2 K2 + . . . + α2,n Kn )]
......
Kn = h[f (xk + θn h, yk + αn,1 K1 + αn,2 K2 + . . . + αn,n Kn )]
yk+1 = yk + h[γ1 K1 + γ2 K2 + . . . γn Kn ] (4.23)

Les coefficients sont d´etermin´es afin que l’ordre soit le plus ´elev´e
possible. On note A la matrice de coefficients (αi,j), Γ le vecteur des
coefficients γi et Θ le vecteur des coefficients θi .
Quand la matrice A est triangulaire inf´erieure stricte, αij = 0 pour j
≥
i, on dit que la m´ethode est explicite. Si seule la partie triangulaire sup
´erieure est nulle, αij = 0 pour j > i, la m´ethode est dite implicite ; sinon elle
est totalement implicite.
Une repr´esentation en forme de tableau des ´equations (4.23) donne

γ1 γ2 ... γn
θ1 α1,1 α2,1α1,2 α2,2. . . α1,n α2,n
θ2 ... ... ...
... αn,1 αn,2 . . .
θn ...
αn,n
Avec cette repr´esentation, la m´ethode Runge-Kutta explicite `a deux
points qui est d’ordre 4 est repr´esent´ee par le tableau suivant

1/6 1/3 1/3 1/6

0 0 0 0 0
1/2 1/2 0 0 0
1/2 0 1/2 0 0
1 0 0 1 0

4.4 M´ethode d’int´egration `a pas variables

4.4.1 Introduction
Un int´egrateur “intelligent” poss`ede une proc´edure de controˆle de la m
´ethode de convergence, c’est `a dire un moyen d’estimer l’erreur commise
par le calcul sur un pas d´int´egration et la possibilit´e de choisir en cons
´equence un nouveau pas si le syst`eme diff´erentiel aborde une r´egion ou` la
fonction prend des valeurs plus importantes. Le calcul de cette estimation
entraˆıne un surcouˆt de calcul, qu’il convient de bien g´erer afin de
minimiser cet effort suppl´ementaire.
L’id´ee la plus simple pour estimer cette erreur consiste `a calculer la
solution donn´ee par un algorithme (de Runge-Kutta d’ordre 4 par exemple)
pour deux pas d´int´egration diff´erents, h et 2h. Soit y(x + 2h) la solution

5
exacte `a x + 2h

5
 4.5 M´ethodes de Runge-Kutta “embarqu´ees”

et y(x + h) la solution exacte `a x + h, on a

y(x + 2h) = y1 + (2h)5φ + O(h6) (4.24)

y(x + 2h) = y2 + 2(h )φ + O(h )
5 6
(4.25)

ou` φ est une fonction qui reste constante sur l’intervalle x, x + 2h `a

l’ordre h5 . La premi`ere ´equation correspond `a une int´egration avec un pas
´egal `a 2h tandis que la seconde correspond `a deux int´egrations successives
avec un pas de h. La diff´erence
∆ = y2 − y1 (4.26)

fournit une estimation de l’erreur commise avec un pas d’int´egration h.

4.5 M´ethodes de Runge-Kutta “embarqu´ees”

Une autre m´ethode pour estimer l’erreur commise par l’utilisation d’un
pas h est due `a Fehlberg. Il utilise le fait qu’en choisissant des valeurs
particuli`eres de γi (voir section 4.3.7), on peut changer l’ordre de
l’´evaluation de la solution pour un pas de temps h donn´e.

6
Σ
y1 = y(x) + γiKi + 0(h6) (4.27)
i=1
6
Σ ∗
y2 = y(x) + γi Ki + 0(h5 ) (4.28)
i=1

ce qui conduit `a une estimation de l’erreur

6
Σ ∗
∆= (γi − γi )Ki (4.29)
i=1

Pour d´eterminer la valeur du pas la plus adapt´ee, on note tout d’abord

que ∆ est calcul´e `a l’ordre h5 . Si on a un pas h1 qui donne une erreur
∆1 , le pas h0 donn´e pour une erreur ∆0 fix´ee `a l’avance, est donn´e par la
relation

..∆ 1/5
h0 = . (4.3
∆ 1
Il est donc possible, pour une valeur de ∆0 donn´ee `a l’avance de
|
diminuer h0 pour obtenir une erreur plus faible ou d’augmenter h0 de
mani`ere raisonnable si ∆1 est inf´erieur en valeur absolue `a |∆0 |.
Une difficult´e apparaˆıt pour ce type de m´ethode quand on consid`ere un
sys- t`eme diff´erentiel `a plusieurs variables. L’estimation de l’erreur est alors
donn´ee `a priori par un vecteur. La g´en´eralisation de la proc´edure ´ebauch
´ee reste possible, mais nous ne d´etaillerons pas ce type de subtilit´e dans ce

5
chapitre.

5
Equations diff´erentielles

10

1
6

2
3
4

-2

0.5 1 1.5 2 2.5 3 3.5

FIG. 4.1 – Sch´ema illustrant le principe de la m´ethode de Burlish-Stoer.

Les lignes en pointill´e correspondent aux int´egrations r´ealis´ees avec
trois pas d’in- t´egration diff´erents. La courbe en trait plein correspond `a la
solution exacte.

4.6 M´ethode de Bulirsh-Stoer

L’id´ee de cette m´ethode repose sur les trois principes suivants : la r

´esolution de l’´equation diff´erentielle pour un accroissement de ∆x est
donn´ee par une fonction qui d´epend de h, mais qui tend vers une limite
finie (ind´ependante de h) ; on cherche `a estimer la valeur exacte du
syst`eme diff´erentiel `a int´egrer en calculant pour diff´erents pas et en
prenant la limite d’un pas tendant vers z´ero.

La figure 4.1 illustre graphiquement les diff´erentes estimations de la

valeur de la fonction solution en utilisant trois pas d’int´egration diff
´erents. La limite d’un pas de temps nul correspond `a la solution exacte
donn´ee par la courbe en trait plein. Ce point est tr`es comparable `a celui
de la m´ethode de Romberg d´ecrite pour l’int´egration.

Le second principe de la m´ethode consiste `a extrapoler cette limite non

pas sur la base de d´eveloppements polynˆomiaux mais de d´eveloppements
en frac- tions rationnelles. Le troisi`eme principe consiste `a utiliser des

5
fonctions d’erreur qui sont paires en pas d´int´egration.

6
 4.7 Conclusion

4.7 Conclusion
Les m´ethodes expos´ees ci-dessus permettent d’obtenir une solution de
plus en plus pr´ecise pour des syst`emes diff´erentiels ou` le nombre de
fonctions n’est pas tr`es important. En simulation de Dynamique Mol
´eculaire, un crit`ere de qualit´e pour le choix d’une m´ethode est le
respect de la propri´et´e d’invariance par renversement du temps, qui est
une propri´et´e satisfaite pour les syst`emes hamiltoniens. Les algorithmes
satisfaisant cette propri´et´e sont appel´es symplec- tiques ; un exemple de ce
type d’algorithme est l´algorithme dit de Verlet dont le principe est d´ecrit
dans le cours Simulation num´erique en Physique Statis- tique.(http
://[Link])

6
Equations diff´erentielles

6
Chapitre 5

Transform´ee de Fourier rapide

Contenu
5.1 Introduction....................................................................... 39
5.2 Propri´et´es ..................................................................................39
5.3 Discr´etisation de la transform´ee de Fourier ...............42
5.3.1 E´ chantillonage ......................................................................42
5.3.2 Transform´ee de Fourier discr`ete .....................................43
5.4 Transform´ee de Fourier rapide ...........................................44

5.1 Introduction
Largement utilis´ee en Physique, la transform´ee de Fourier d’une
fonction d´ependant d’une variable (par exemple du temps) est devenue si
naturelle que sa repr´esentation graphique est g´en´eralement aussi utile,
voire plus, que celle de la fonction elle-mˆeme. Apr`es un rappel des propri´et
´es ´el´ementaires et fonda- mentales des transform´ees de Fourier pour la r
´esolution de probl`emes math´ema- tiques, nous allons pr´esenter le principe
de la m´ethode num´erique de l’´evaluation de cette transform´ee. Plus sp
´ecifiquement, depuis pr`es de 40 ans est apparu un algorithme performant
pour le calcul de la transform´ee de Fourier dont le temps de calcul varie
essentiellement comme N ln2 (N ) ou` N est le nombre de points ou` la
fonction f a ´et´e ´evalu´ee. Par opposition `a une approche trop na¨ıve,
ou` le nombre d’op´erations croˆıt comme N 2 cette m´ethode a re¸cu le nom
de trans- form´ee de Fourier rapide (FFT, Fast Fourier Transform, en
anglais), que toute biblioth`eque math´ematique propose g´en´eralement `a
son catalogue.

5.2 Propri´et´es
Soit une fonction f d´efinie sur R, on appelle transform´ee de Fourier
de f , la fonction fˆ
∫

6
fˆ(ν) =
+∞
2πiνt
f dt (5.1)
(t)e
−∞

6
Transform´ee de Fourier rapide

La transform´ee de Fourier inverse est d´efinie comme

+∞
f1(t) = ∫ f −2πiνt
dν (5.2)
ˆ
(ν)e
−∞

Pour une fonction f int´egrable, il y a identit´e entre la fonction f1 et la fonction

f , hormis ´eventuellement sur un support de R de mesure nulle.
Une autre d´efinition de la transform´ee de Fourier rencontr´ee dans la litt
´era- ture est celle qui correspond `a une repr´esentation en pulsation au
lieu de celle en fr´equence donn´ee ci-dessus.

fˆ(ω) = ∫ +∞
iωt
f dt (5.3)
(t)e
−∞

La transform´ee de Fourier inverse correspondante est d´efinie comme1

1
+∞
f1(t) = ∫ f (ω)e−iωt dω (5.6)
2π
−∞ ˆ
Avec le logiciel Maple, la transform´ee de Fourier est d´efinie d’une
mani`ere encore diff´erente !
∫ +∞
fˆ(ω) = f −iωt
dt (5.7)
(t)e
−∞
avec bien entendu la transform´ee de Fourier inverse correspondante
1
+∞
f1(t) = ∫ f (ω)eiωt dω (5.8)
2π
−∞ ˆ
Compte tenu des diff´erentes d´efinitions des transform´ees de Fourier
et en l’absence de documentation pr´ecise dans certaines biblioth`eques math
´ematiques, il est fortement recommand´e de tester la transform´ee de Fourier
num´erique que l’on va utiliser, en testant une fonction dont la transform´ee de
Fourier est connue analytiquement pour d´eterminer de mani`ere simple la
convention utilis´ee dans la biblioth`eque disponible. Au del`a du principe g
´en´eral de la transform´ee de Fourier rapide, il existe des optimisations
possibles compte tenu de l’architec- ture de l’ordinateur sur lequel le
programme sera ex´ecut´e. Cela renforce l’aspect totalement inutile de la r´e
´ecriture d’un code r´ealisant la dite transform´ee, ce code serait en g´en
´eral tr`es peu performant compar´e aux nombreux programmes disponibles.
Dans le cas ou` la fonction f poss`ede des sym´etries, sa transform´ee de
Fourier pr´esente des sym´etries en quelque sorte duales de la fonction initiale
1
Pour sym´etriser les expressions, la transform´ee de Fourier est parfois d´efinie comme
∫ +∞
ˆ = √ 1 iωt
f (ω)
2π et la transform´ee de Fourier inverse comme

6
f (t)e −∞
dt (5.4)
∫ +∞
1 ˆf
−iωt
f1(t) = dω
−

6
 5.2 Propri´et´es

– Si f (t) est r´eelle, la transform´ee de Fourier des fr´equences n

´egatives est complexe conjugu´ee de celle des arguments positifs fˆ(
∗ −
ν) = (fˆ(ν)) .
– Si f (t) est imaginaire, la transform´ee de Fourier des fr´equences n
´egatives est ´egale `a l’oppos´ee du complexe conjugu´e de celle des
∗
arguments positifs fˆ( ν) = (fˆ(ν)) .
—
– Si f (t) est paire, la transform´ee de Fourier est aussi paire, fˆ( ν) = fˆ(ν).
−
– Si f (t) est impaire, la transform´ee de Fourier est aussi impaire, fˆ( ν) =
−
fˆ(ν).
−
– Si f (t) est paire et r´eelle, la transform´ee de Fourier l’est aussi.
– Si f (t) est impaire et r´eelle, la transform´ee de Fourier est
imaginaire et impaire.
– Si f (t) est paire et imaginaire, la transform´ee de Fourier l’est aussi.
– Si f (t) est impaire et imaginaire, la transform´ee de Fourier est r
´eelle et impaire.
Des propri´et´es compl´ementaires sont associ´ees aux op´erations de
translation et dilatation
– Soit f (at), sa transform´ee de Fourier est donn´ee par 1 fˆ( ν ).
|a| a
– Soit 1 f ( t ), la transform´ee de Fourier est donn´ee par fˆ(bν).
|b| |b|
– Soit f (t t0 ), la transform´ee de Fourier est donn´ee par fˆ(ν)e2πiνt0 .
−
– Soit f (t)e−2πiνt0 , la transform´ee de Fourier est donn´ee par fˆ(ν ν0 ).
−
Une identit´e particuli`erement utile concernant les transform´ees de Fourier
concerne la distribution de Dirac δ : soit f une fonction continue d´efinie sur
R, La distribution δ est d´efinie `a partir de la relation suivante
∫
dxf (x)δ(x − x0) = f (x0) (5.9)
+∞
−∞
Ainsi, on en d´eduit facilement
que
∫ +∞ 2πkix 2πkix

dxδ(x − x0)e =e 0
(5.10)
−∞
et donc que
∫ +∞
dxδ(x)e 2πikx = 1 (5.11)
−∞
Par inversion de cette relation, on a la repr´esentation int´egrale de la distribution
δ
∫
δ(x) = +∞
dke −2πikx (5.12)
−∞
De mani`ere ´equivalente, on montre que
1
∫ +∞
δ(x) = dke −ikx (5.13)
2π −∞
La propri´et´e sans doute la plus importante concernant les transform´ees de
Fou- rier concerne la convolution. Soit deux fonctions f et g d´efinie sur C, on
d´efinit la convolution de f avec g not´ee f ∗ g comme

6
 ∫
+∞
( )( ) = ( ) ( ) (5.14)
f∗g t f τ g t − τ dτ
−∞

6
Transform´ee de Fourier rapide

On v´erifie imm´ediatement que f g = g f . L’op´erateur de convolution

∗ ∗
est commutatif, car la multiplication de deux nombres complexes l’est aussi.
La transform´ee de Fourier de la convolution de deux fonctions f et g est
´egale au produit de leurs transform´ees de Fourier.

(^f ∗ g)(ν) = fˆ(ν)gˆ(ν) (5.15)

Soit une fonction r´eelle f , on d´efinit la fonction d’autocorr´elation

comme

Cf (t) = ∫
duf (u)f (t + u) (5.16)
+∞
−∞

La transform´ee de Fourier de la fonction Cf est alors donn´ee par

Cˆf (ν) = fˆ(ν)fˆ(−ν) (5.17)

∗
Compte tenu du fait que f est r´eelle, on a fˆ(−ν) = (fˆ(ν)) , ce qui donne

C^f (ν) = |f ˆ(ν)|2 (5.18)

Cette relation est appel´ee th´eor`eme de Wiener-Khinchin.

Le th´eor`eme de Parseval donne que
∫ +∞ ∫ +∞
|fˆ(ν)| dν
2 2
|f (t)| dt = (5.19)
−∞ −∞

5.3 Discr´etisation de la transform´ee de Fourier

5.3.1 E´ chantillonage
Dans la situation la plus fr´equemment rencontr´ee, la fonction f est
´echan- tillonn´ee `a intervalle r´egulier. Soit ∆ le pas s´eparant deux
“mesures” cons´ecu- tives, on a

fn = f (n∆) avec n = . . . , −2, −1, 0, 1, 2, 3, . . . (5.20)

Pour le pas de l’´echantillonnage ∆, il existe une fr´equence critique, appel

´ee fr´equence de Nyquist,
1
νc = (5.21)
2∆
au del`a de laquelle il n’est pas possible d’avoir une information sur le
spectre de la fonction ´echantillonn´ee. Supposons que la fonction soit une
sinuso¨ıde de fr´equence ´egale `a celle de Nyquist. Si la fonction est
maximale pour un point, elle est minimale au point suivant et ainsi de
suite.
La cons´equence de ce r´esultat est que si on sait que le support de la
trans- form´ee de Fourier est strictement limit´e `a l’intervalle [ νc , νc ], la
−
fonction f (t) est alors donn´ee par la formule

6
 +∞
Σ sin(2πνc(t − n∆))
f (t) = ∆ f (5.22)
n=−∞
n π(t − n∆)

7
 5.3 Discr´etisation de la transform´ee de Fourier

Si la fonction f (t) est multipli´ee par la fonction e2πν1 t avec ν1 qui est
un multiple de 1/∆, les points ´echantillonn´es sont exactement les
mˆemes. Une cons´equence de cette propri´et´e est l’apparition du ph
´enom`ene de duplication (“aliasing” en anglais). Cela se manifeste de la
mani`ere suivante : supposons que la transform´ee de Fourier exacte d’une
fonction ait un spectre plus large que celui donn´e par l’intervalle de
Nyquist, le spectre situ´e `a droite se replie `a gauche et celui `a droite se
replie `a gauche.

5.3.2 Transform´ee de Fourier discr`ete

Soit un nombre fini de points ou` la fonction f a ´et´e ´echantillonn´ee

fk = f (tk) avec k = 0, 1, 2, . . . , N − 1 (5.23)

Avec N nombres fournis, il semble ´evident que l’on peut obtenir N fr

´equences pour la transform´ee de Fourier. Soit la suite

νn = n avec n = −N/2, . . . , 0, . . . , N/2 (5.24)

N∆
Le nombre de fr´equences est a priori ´egale `a N + 1, mais comme les bornes
in- f´erieure et sup´erieure de l’intervalle en fr´equence correspondent aux
bornes d´e- finies par la fr´equence critique de Nyquist, les valeurs obtenues
−
pour n = N/2 et n = N/2 sont identiques et nous avons donc bien N points
ind´ependants.
Soit un nombre fini de points not´e hk , on d´efinit la transform´ee de
Fourier discr`ete comme
hˆ n = 1 NΣ−1
N hk e2πikn/N (5.25)
k=0
La transform´ee de Fourier discr`ete inverse est donn´ee par
N −1
Σ −
hk = hne (5.26)
n=0

Ces transform´ees de Fourier sont d´efinies ind´ependamment des abscisses

cor- respondant aux valeurs originales ou` la fonction h a ´et´e ´evalu´ee. Il
est toutefois possible de relier la transform´ee de Fourier discr`ete `a la
transform´ee de Fourier de d´epart par la relation
fˆ(νn) ' ∆fˆn (5.27)
2
Notons que les propri´et´es ´etablies pour les transform´ees de Fourier au d
´ebut de ce chapitre, selon le cas ou` la fonction est paire ou impaire, r´eelle ou
imaginaire se transposent compl`etement au cas de la transform´ee de
Fourier discr`ete. Les
2
La transform´ee de Fourier continue est approch´ee de la mani`ere suivante
−1
∫+∞
fˆ(νn ) = f (t)e
−∞
N
Σ
7
 2πiνnt dt (5.28)
≈∆ fke2πiνntk (5.29)
k=0

7
Transform´ee de Fourier rapide

r`egles concernant la convolution se retrouvent bien ´evidemment et par exemple

l’expression du th´eor`eme de Parseval est alors
N −1 −1
Σ 1 NΣ
|h |hk |
2
(5.30)
k| =
2
ˆ
k=0 N k=0
Le passage de la transform´ee de Fourier continue `a la transform´ee de
Fourier discr`ete est en fait un changement de mesure, ce qui explique que les
propri´et´es
´etablies pour l’une restent vraies pour la seconde.
On peut noter que les diff´erences entre la transform´ee de Fourier
discr`ete et son inverse sont au nombre de deux : un changement de signe dans
l’exponen- tielle complexe et un facteur de normalisation en 1/N . Ces
faibles diff´erences expliquent pourquoi les proc´edures sont souvent
identiques pour le calcul de la transform´ee de Fourier et de son inverse dans
les biblioth`eques math´ematiques.

5.4 Transform´ee de Fourier rapide

L’id´ee de base de la transform´ee de Fourier rapide repose sur la
remarque suivante : posons
W = e2πi/N (5.31)
La composante de Fourier de h s’exprime alors comme
N −1
Σ
hn = nkˆW hk (5.32)
k=0

Ainsi le vecteur h de composante hk est multipli´e par une matrice de

coefficients ank = W nk . Sans utiliser d’astuces particuli`eres, le temps de
calcul pour une telle op´eration est alors en N 2 . Pour am´eliorer de
mani`ere spectaculaire la rapidit´e de ce traitement, on note que la
transform´ee de Fourier discr`ete de longueur N peut ˆetre ´ecrite comme la
somme de deux transform´ees de Fourier discr`ete chacune de longueur N/2.
En effet, en supposant que le nombre N est pair, on peut s´eparer la
contribution des termes pairs et celle des termes impairs dans l’´equation
−1
NΣ
Fk = e fj (5.33)
j=0
2πijk/N
N/2−1 N/2−1
= Σ Σ
e2πi(2j)k/N f2j + e2πi(2j+1)k/N f2j+1 (5.34)
j=0 j=0
N/2−1 N/2−1
= Σ 2πijk/(N/2) k
Σ
e f2j + W e2πijk/(N/2)f2j+1 (5.35)
j=0 j=0
=F p + W k F i (5.36)
k k

7
ou` F p repr´esente la ki`eme composante de Fourier pour les composantes
paireskde la fonction de d´epart et F i repr´esente la ki`eme composante de
Fourier pour k

7
 5.4 Transform´ee de Fourier rapide

000 000
000
001 001
001

010 010 010

011
011 011
100 100
100

101 101 101

110 110
110

111 111
111

(a) (b)

FIG. 5.1 – Sch´ema illustrant le principe de la mise en ordre du tableau

par renversement de bits. (a) Sur un tableau de 8 ´el´ements ´etiquett´es en
binaire, on note les op´erations de d´eplacements. (b) Bilan des d
´eplacements `a effectuer : deux op´erations d’´echange d’´el´ements.

les composantes impaires de la fonction de d´epart. A noter que ces deux

com- posantes sont p´eriodiques en k avec une p´eriode N/2.
Pour it´erer ce type de proc´edure, nous allons supposer que le nombre
de points de la fonction de d´epart est dor´enavant une puissance de deux. On
peut alors exprimer chaque composante de Fourier en fonction de deux
nouvelles composantes de Fourier sur un intervalle de longueur N/4 et ainsi de
suite. Une fois obtenue un intervalle de longueur 1, on a la transform´ee de
Fourier pour cet intervalle qui est ´egale au nombre fn correspondant. Le
nombre d’indices pour caract´eriser ce nombre est ´egal `a ln2 (N )

iipip...ppi
F k
= fn (5.37)

On remarque donc que pour chaque valeur de n on a un alphabet fini d’indice

iipip . . . ppi constitu´e de i et de p. Cette correspondance est biunivoque.
Si on pose que i = 1 et p = 0, on a l’expression invers´ee en binaire de
chaque valeur de n.
La proc´edure de transform´ee de Fourier rapide est alors constitu´ee d’une
pre- mi`ere op´eration de tri pour r´eordonner la matrice selon la lecture de
l’alphabet g´en´er´e, puis on calcule successivement les composantes de
Fourier de longueur 2, puis 4, jusqu’`a N . La figure 5.1 illustre le nombre
d’op´erations `a effectuer pour inverser les ´el´ements du tableau selon le
m´ecanisme du renversement de bits.

7
Transform´ee de Fourier rapide

De nombreuses variantes de cet algorithme original sont disponibles, qui

per- mettent d’obtenir la transform´ee de Fourier pour un nombre de points
diff´erent d’une puissance de deux.

7
Chapitre 6

Alg`ebre lin´eaire

Contenu
6.1 Introduction....................................................................... 47
6.2 E´ limination de Gauss-Jordan ..............................................49
6.2.1 Rappels sur les matrices.......................................................49
6.2.2 M´ethode sans pivot .............................................................49
6.2.3 M´ethode avec pivot ............................................................50
6.3 E´ limination gaussienne avec substitution .......................50
6.4 D´ecomposition LU ...................................................................51
6.4.1 Principe...................................................................................51
6.4.2 R´esolution d’un syst`eme lin´eaire ...................................52
6.5 Matrices creuses.....................................................................53
6.5.1 Introduction............................................................................53
6.5.2 Matrices tridiagonales...........................................................54
6.5.3 Formule de Sherman-Morison..............................................54
6.6 D´ecomposition de Choleski ..................................................54
6.7 Conclusion.......................................................................... 55

6.1 Introduction
Les deux chapitres qui suivent concernent le traitement num´erique des
ma- trices Ce premier chapitre est consacr´e aux m´ethodes employ´ees pour r
´esoudre les quatre taˆches les plus fr´equemment rencontr´ees pour les
syst`emes lin´eaires. Par d´efinition, on peut ´ecrire celles-ci comme

A.x = b (6.1)
ou` A est une matrice M × N
a11 a12 ... a1N
a a22 ... a2N
21

A=  . .. . 
... (6.2)
.
.. 

 . . . . 

7
aM 1 aM 2 . . . aMN

7
Alg`ebre lin´eaire

x est un vecteur colonne de M ´el´ements et b un vecteur colonne de N ´el´ements.

   
b1 x1
 b  x
. .
b= .  x= . (6.3)
. .
  
   
bN xM

Si N = M , il y autant d’´equations que d’inconnues et si aucune des ´equations

n’est une combinaison lin´eaire des N −1 autres, la solution existe et est unique. 1
Hormis pour des valeurs de N tr`es petites (typiquement N ≤ 4), ou` les
formules sont simples `a exprimer, il est g´en´eralement n´ecessaire de proc´eder
`a un calcul num´erique pour obtenir la solution d’un syst`eme
d’´equations li- n´eaires. Comme nous allons le voir ci-dessous, si la matrice
A n’est pas une matrice creuse (matrice dont la majorit´e des ´el´ements
est nulle), il est n´eces- saire d’appliquer une m´ethode g´en´erale qui
revient en quelque sorte `a inverser la matrice A (ou `a la factoriser), ce qui n
´ecessite un grand nombre d’op´erations qui augmente comme le cube de la
dimension lin´eaire de la matrice, N 3 .
Mˆeme si la proc´edure num´erique est cens´ee conduire `a une solution
dans le cas d’une matrice non singuli`ere (det(A) = 0), l’accumulation
/
d’erreurs d’ar- rondi, souvent li´ee `a des soustractions de nombres voisins
fournit un vecteur x erron´e. Dans le cas d’une matrice dont le d
´eterminant est tr`es voisin de z´ero (matrice presque singuli`ere), les
solutions obtenues peuvent ˆetre aussi fausses. L’´enorme biblioth`eque de
sous-programmes pour les probl`emes d’alg`ebre li- n´eaire montre
l’importance de ces probl`emes dans le calcul num´erique et la n´ecessit´e
de choisir une m´ethode sp´ecifique d`es que la matrice a des propri´et´es
particuli`eres.
Parmi les taˆches typiques d’alg`ebre lin´eaire hormis la r´esolution
d’un sys- t`eme lin´eaire, on peut citer
– la d´etermination des solutions d’un ensemble de syst`emes, par
exemple [Link] = bj ou` xj et bj sont des vecteurs `a N composantes
et j un indice parcourant un ensemble fini d’entiers,
– le calcul de l’inverse de A : A−1,
– le calcul du d´eterminant de A.
De mani`ere encore plus aigu¨e que dans les chapitres pr´ec´edents, une m
´ethode de force brute est bien moins efficace pour r´esoudre un probl`eme
d’alg`ebre li- n´eaire qu’une m´ethode sp´ecifique. Nous allons donc voir tout
d’abord le principe d’une m´ethode g´en´erale pour une matrice quelconque,
puis consid´erer quelques unes des techniques sp´ecifiques qui d´ependent de
la structure des matrices.

7
 1
Si N > M , les ´equations sont ind´ependantes entre elles et il n’y a pas de solution. Si
N < M , il y a une ind´etermination et il existe une infinit´e de solutions.

8
 6.2 E´ limination de Gauss-Jordan

6.2 E´ limination de Gauss-Jordan

6.2.1 Rappels sur les matrices
Il est utile de remarquer que le calcul des quatre taˆches pr´ec´edentes
peut ˆetre a priori ex´ecut´e de mani`ere similaire : en effet, si on consid`ere
par exemple un syst`eme 4 × 4, on a
       
a11 a12 a13 a14 x1 t1 y11 y12 y13 y14

 a21 a22 a23 a24 . x2  H  t2 H  y21 y22 y23 

=a31 a34 a   y y y y
132 133 x3  31 132
y41 y42 y43 y44
133 34
a41 a42 a43 x t3
  t4 
 
a44 b1   1 0 0 0
=   H  c   (6.4)
bb23b c1 H 00 0 1 10
  c 0
3 
ou` l’operateur´ H designe´ l’operateur´ reunion´ de colonnes. Ainsi on voit que
l’inversion d’une matrice est identique `a la r´esolution d’un ensemble de N
sys- t`emes lin´eaires pour lequel le vecteur situ´e dans le membre de droite a
tous les
´el´ements nuls sauf un seul qui est diff´erent pour chaque colonne.
Les trois r`egles ´el´ementaires que l’on peut r´ealiser sur les matrices
sont les suivantes.
– E´ changer les lignes de A et de b (ou c ou de la matrice identit´e) et en
gardant x (ou t ou Y ) ne change pas la solution du syst`eme lin
´eaire. Il s’agit d’une r´e´ecriture des ´equations dans un ordre diff
´erent.
– De mani`ere similaire, la solution du syst`eme lin´eaire est inchang´ee si
toute ligne de A est remplac´ee par une combinaison lin´eaire d’elle-
mˆeme et des autres lignes, en effectuant une op´eration analogue
pour b (ou c ou la matrice identit´e).
– L’´echange de deux colonnes de A avec ´echange simultan´e des
lignes cor- respondantes de x (ou t ou Y ) conduit `a la mˆeme
solution. Toutefois, si on souhaite obtenir la matrice Y dans l’ordre
initialement choisi, il est n´ecessaire de proc´eder `a l’op´eration
inverse une fois la solution obtenue.

6.2.2 M´ethode sans pivot

Pour la simplicit´e de l’expos´e, on consid`ere tout d’abord une matrice
dont les ´el´ements diagonaux sont strictement diff´erents de z´ero.
On multiplie la premi`ere ligne de A par 1/a11 (et la ligne
correspondante de b ou de la matrice identit´e). On soustrait a21 fois la 1`ere
ligne `a la deuxi`eme ligne, a31 fois la 1`ere ligne `a la troisi`eme ligne et ainsi
de suite jusqu’`a la derni`ere ligne. La matrice A a maintenant la structure
suivante
 J J J 
1 a 12 a 13 a 14

8
 J J J
0 a 22 a 23 a 24
(6.5)
0 aJ aJ J
 32 33 a 34
J J J
0 a 42 a 43 a 44


8
Alg`ebre lin´eaire

ou` les coefficients aJij s’expriment en fonction des coefficients ou` les coefficients
aij.
On multiplie alors la deuxi`eme ligne par 1/aJ22 , puis on soustrait aJ12
fois la deuxi`eme ligne `a la premi`ere ligne, aJ32 fois la deuxi`eme ligne `a
la troisi`eme ligne et ainsi de suite jusqu’`a la derni`ere ligne. La matrice
A a maintenant la
structure suivante

1 0 aJ1J3 aJ1J4
0 1 aJ2J J3 aJJ2JJ4 (6.6)
0 0 a 33 a 3 4

0 0 aJ4J3 aJ4J4


On it`ere le proc´ed´e jusqu’`a la derni`ere ligne et on obtient alors la
matrice identit´e. On peut obtenir alors facilement la solution du syst`eme
d’´equations.

6.2.3 M´ethode avec pivot

La m´ethode pr´ec´edente souffre du d´efaut de ne s’appliquer qu’aux
matrices dont tous les ´el´ements de la diagonale sont non nuls. Or, une
matrice n’est pas n´ecessairement singuli`ere si un seul ´el´ement de la
diagonale est nul. Pour permettre une plus grande g´en´eralit´e de la m
´ethode, on ajoute aux op´erations pr´ec´edentes la troisi`eme r`egle
´enonc´ee ci-dessus. On peut en effet, en utilisant l’inversion des colonnes,
placer sur la diagonale un ´el´ement non nul de la ligne (si tous les ´el
´ements d’une ligne sont nuls, cela signifie que le d´eterminant de la matrice
est nul et que la matrice est singuli`ere, ce qui a ´et´e exclu par hypoth`ese). De
plus, on peut montrer que la recherche du plus grand ´el´ement de la ligne
pour faire la permutation des colonnes rend la proc´edure bien plus stable,
en particulier en augmentant de mani`ere importante la pr´ecision num
´erique des solutions.
Ce type d’algorithme est disponible dans de nombreuses biblioth`eques
et il est parfaitement inutile de chercher `a r´e´ecrire un code qui n
´ecessite un temps important `a la fois pour sa r´ealisation et pour sa
validation.

6.3 E´ limination gaussienne avec substitution

Pour la r´esolution stricte d’un syst`eme lin´eaire, il n’est pas n´ecessaire
d’ob- tenir une matrice diagonale comme celle obtenue par la m´ethode pr
´ec´edente. Ainsi en effectuant `a chaque ligne la soustraction des lignes
situ´ees au dessous de la ligne dont le coefficient de la diagonale vient
d’ˆetre r´eduit `a l’unit´e, on obtient une matrice triangulaire sup´erieure
dont la structure est la suivante

8
 J J J J 
a 11 a 12 a 13 a 14
   
J J J
0 a 22 a 23 a 24 xx21 bJ1
 J J     (6.7)
. 0 b 20J a 133
 x = b
a 34 3
J
  3

8
 6.4 D´ecomposition LU

Il est possible de faire la r´esolution de ce syst`eme lin´eaire par une proc

´edure de substitution successive
J
b4
x= J
a4 (6.8)
4
4
1 (bJ − x aJ ) (6.9)
x =3 J 3 4 34
a 33
et de mani`ere g´en
´erale
N
1 Σ
x = (b −
i i
aJ x ) (6.10)
a
i j
j=i
Cette m´ethode est avantageuse dans le cas des syst`emes lin´eaires car le
nombre d’op´erations `a effectuer est d’ordre N 2 , contrairement `a la m
´ethode pr´ec´edente qui n´ecessite de calculer compl`etement l’inverse de
la matrice et n´ecessite un temps de calcul en N 3.
S’il faut r´esoudre un ensemble de syst`emes lin´eaires comprenant N
termes, cela est alors ´equivalent `a r´esoudre l’inversion d’une matrice et l’on
retrouve un temps de calcul de l’ordre de N 3.

6.4 D´ecomposition LU
Les m´ethodes pr´ec´edentes n´ecessitent de connaˆıtre `a l’avance le
membre de gauche de l’´equation (6.1). Les m´ethodes qui suivent consiste
`a r´e´ecrire la matrice A afin que la r´esolution du syst`eme d’´equations
soit ex´ecut´ee plus faci- lement.

6.4.1 Principe
Nous allons montrer que toute matrice N N peut se d´ecomposer de
×
la mani`ere suivante.
A = L.U (6.11)
ou` L est une matrice triangulaire inf´erieure et U une matrice triangulaire
sup´e- rieure, soit
α11 0 0 ... 0 
α21 α22 .0 . . . . . . . . . . . . 0 
 . . . ... (6.12)

L . . . 
. 
αN1. αN2. . . .. α(N .−1)(N. . −1) αNN0
et
β11 β12 β13 ... β1N 
0 β22 β. 23 . . . . . . . . . . . . β2N 
. . ... (6.13)
.

U . . . . 
. . . β β
0 0 . . . (N −1)(N0 −1) (N

8
Alg`ebre lin´eaire

En effectuant la multiplication matricielle de L par U , on obtient les relations

suivantes
i
Σ
aij = αilβlj i≤j (6.14)
l=1
j
Σ
aij = αilβlj i>j (6.15)
l=1

Ce syst`eme d’´equations lin´eaires donne N 2 ´equations et le nombre

d’inconnues est 2 N (N +1) . Ce syst`eme est donc surd´etermin´e. On peut
2
donc choisir N
´equations suppl´ementaires afin d’obtenir une et une seule solution. On fixe
donc la valeur de la diagonale de la matrice L

αii = 1 (6.16)

pour i [1, N ].
∈
L’algorithme de Crout permet de calculer simplement les N 2 coefficients
restants, en utilisant les relations suivantes
i−1
Σ
βij = aij − αikβkj i≤j (6.17)
k=1
j−1
Σ
αij 1
= βjj (aij — αik βkj ) i≥j+1 (6.18)
k=1

Il faut noter que les sous-programmes cr´eant cette d´ecomposition s’appuie

sur la recherche du meilleur pivot afin que la m´ethode soit stable.2

6.4.2 R´esolution d’un syst`eme lin´eaire

La r´esolution d’un syst`eme lin´eaire devient tr`es simple en introduisant
le vecteur y.

A.x = L.U.x (6.20)

= L.(U x) = b (6.21)

Soit
L.y = b (6.22)
2
Sachant que le nombre d’inconnues est devenu ´egale a` N 2 , on peut ´ecrire les deux
matrices sous la forme d’une seule matrice de la forme
β11 β12 β13 ... β1N 
α21 β22 β23........................................β2N
 . . . . 
 . . . ... . 
 (6.19)
 
. . . β(N −1)(N −1) 0
αN 1 αN 2 ... ... βNN

8
en se rappelant que les ´el´ements diagonaux de la matrice α ont ´et´e choisis ´egaux a` un.

8
 6.5 Matrices creuses

U.x = y (6.23)
Chaque syst`eme peut ˆetre r´esolu par une proc´edure de substitution.

b1
y1 = (6.24)
i−1
α11 Σ
1
yi = (bi — αij yj ) (6.25)
j=1
αii

(substitution `a partir du premier ´el´ement car la matrice est triangulaire

inf´e- rieure). On obtient la solution pour le vecteur x en utilisant

yN
xN = βNN (6.26)
1
x= (y — N
i
βii
i βij xj ) (6.27)
Σ

j= i+
1

car la matrice est triangulaire sup´erieure.

Une fois effectu´ee la d´ecomposition LU , le calcul du d´eterminant d’une
ma- trice devient tr`es simple. On sait que le d´eterminant d’un produit de
matrices est ´egale au produit des d´eterminants. De plus pour une matrice
triangulaire, le d´eterminant de la matrice est ´egal au produit des ´el
´ements de sa diagonale. Comme les ´el´ements de la diagonale de L ont ´et
´e choisis ´egaux `a 1, le d´etermi- nant de cette matrice est donc ´egal `a
1, et le d´eterminant de A est donc ´egal
`a
N
Y
det(A) = βjj (6.28)
j=1

6.5 Matrices creuses

6.5.1 Introduction
On appelle matrice creuse une matrice dont la plupart des ´el´ements sont
´egaux `a z´ero. Si de plus, la structure des ´el´ements non nuls est
simple, il n’est pas n´ecessaire de r´eserver une quantit´e de m´emoire
´egale `a celle de la matrice compl`ete. Des algorithmes sp´ecifiques
permettent de r´eduire le temps de calcul de mani`ere consid´erable. Parmi
les cas simples de matrices creuses, citons les matrices
– tridiagonales : les ´el´ements de matrice non nuls sont sur la
diagonale et de part et d’autre de celle-ci sur les deux lignes
| − 8
 adjacentes. On a aij = 0 pour i j > 1,
– diagonales par bande de largeur M : les ´el´ements de matrice tels que
i j > M sont nuls aij = 0,
| −
– simplement ou doublement bord´ees : par rapport `a la d´efinition pr
´ec´e- dente, des ´el´ements non nuls suppl´ementaires existent le
long des lignes ou colonnes du bord de la matrice.

8
Alg`ebre lin´eaire

6.5.2 Matrices tridiagonales

Soit le syst`eme suivant
b    
1 c1 0 ... 0 x1 y1
a2 b2 c2 ... 0   
  x2
0 . . . (6.29)
 . . . . . .
 . = .  . 
 
bN −1 c−1
N  y
0
0 . 0. . .. . . . . . x  
   . 
Il y a un vaste de choix
Noter qu’avec cette notation a1 et cN ne sont pas d´efinis.
.
de biblioth`eques disponibles pour calculer les solutions qui prennent un
temps de calcul proportionnel `a N . De mani`ere g´en´erale, on peut
obtenir aussi avoir un algorithme proportionnel `a N pour une matrice `a
bandes. Le pr´efacteur de l’algorithme est proportionnel `a M .

6.5.3 Formule de Sherman-Morison

Supposons que l’on ait une matrice A dont on a facilement calcul´e
l’inverse (cas d’une matrice tridiagonal). Si on fait un petit changement dans A
en modi- fiant par exemple un ou quelques ´el´ements de l’ensemble de la
matrice, peut-on calculer encore facilement l’inverse de cette nouvelle matrice ?

B =A+u⊗ v (6.30)

ou le symbole d´esigne le produit ext´erieur. u v repr´esente une matrice

⊗
dont l’´el´ement ij est le produit de la i`eme composante de u par la j`eme
composante de v.
La formule de Sherman-Morison donne l’inverse de B

B−1 = (1 + A−1.u ⊗ v)−1.A−1 (6.31)

−1 −1 −1 −1
= (1 − A .u ⊗ v + A .u ⊗ v.A .u ⊗ v − . . .).A (6.32)
= A−1 − A−1.u ⊗ vA−1(1 − λ + λ2 + . . .) (6.33)

= A−1 − (A−1.u) ⊗ (v.A−1)

(6.34)
1+λ
ou` λ = v.A−1 u. On a utilis´e l’associativite´ du produit exterieur´ et produit
interne.
Posons z = A−1u et w = (A−1)T v, on a λ = v.z et

−1 −1 z⊗w
= (6.35)
B A 1+λ

6.6 D´ecomposition de Choleski

Une matrice est d´efinie positive sym´etrique (aij = aji ) quand ∀x, on a

9
x.A.x > 0 (6.36)

9
 6.7 Conclusion

quel que soit x. Il existe alors une matrice L triangulaire inf´erieure telle que

A = [Link] (6.37)

Les ´el´ements de matrice sont d´etermin´es par les relations

i−1
Σ
2 1/2 ik
Lk=1= (a − L )
i−1
(6.38)ii ii
Lji 1
= Li (aij — ΣLik Ljk ) (6.39)
i k=1

avec j [i + 1, N ]. Une fois cette construction r´ealis´ee, on peut r

´esoudre∈ un syst`eme lin´eaire simplement en utilisant la proc´edure de
substitution pr´ec´edem- ment d´etaill´ee dans la section 6.4.2.

6.7 Conclusion
Avant de r´esoudre des syst`emes lin´eaires de grande dimension, il est imp
´era- tif de commencer par une analyse des propri´et´es de la matrice afin de d
´eterminer la m´ethode la plus adapt´ee afin d’obtenir une solution avec
une pr´ecision cor- recte et pour un temps de calcul qui sera minimal. Les
diff´erentes m´ethodes pr´esent´ees dans ce chapitre ne sont qu’une
introduction `a ce tr`es vaste sujet.

9
Alg`ebre lin´eaire

9
Chapitre 7

Analyse spectrale

Contenu
7.1 Introduction....................................................................... 57
7.2 Propri´et´es des matrices .......................................................58
7.3 M´ethodes directes ....................................................................60
7.3.1 M´ethode de Jacobi ..............................................................60
7.3.2 R´eduction de Householder .................................................62
7.3.3 Algorithme QL...................................................................64
7.3.4 Factorisation de Schur...........................................................65
7.4 M´ethode it´eratives .................................................................66
7.4.1 M´ethodes des puissances ....................................................66
7.4.2 M´ethode de Lanczo`s .........................................................67

7.1 Introduction
Ce chapitre est consacr´e aux op´erations que l’on peut r´ealiser sur
des ma- trices. Plus sp´ecifiquement, nous allons nous int´eresser `a la d
´etermination des valeurs propres et/ou vecteurs propres correspondants.
Tout d’abord, quelques rappels utiles pour la suite de ce chapitre : soit une
matrice carr´ee A de dimen- sion N , on appelle un vecteur propre x associ´e
`a la valeur propre λ, un vecteur qui satisfait la relation
A.x = λx (7.1)
Si x est un vecteur propre, pour tout r´eel
0, αx est aussi un vecteur propre
α
avec la mˆeme valeur propre λ.
Les valeurs propres d’une matrice peuvent ˆetre d´etermin´ees comme
les ra- cines du polynˆome caract´eristique de degr´e N

det(A − λ.1) = 0 (7.2)

ou` 1 d´esigne la matrice identit´e.

Compte tenu du fait que dans C, tout polynˆome de degr´e N a N racines,
la matrice A poss`ede N valeurs propres complexes.
9
Analyse spectrale

Quand une racine du polynˆome caract´eristique est multiple, on dit

que la valeur propre est d´eg´en´er´ee, et la dimension de l’espace associ´e
`a cette valeur propre est sup´erieure ou ´egale `a deux.
La d´etermination des valeurs propres d’une matrice `a partir de
l’´equation caract´eristique n’est pas efficace sur le plan num´erique. Des m
´ethodes plus adap- t´ees sont expos´ees dans la suite de ce chapitre.
Si une matrice poss`ede une valeur propre nulle, la matrice est dite singuli`ere.
– Une matrice est sym´etrique si elle est ´egale `a sa transpos´ee
A = AT (7.3)
aij = aji ∀i, j (7.4)
– Une matrice est hermitienne ou auto-adjointe si elle est ´egale au
complexe conjugu´e de sa transpos´ee.
A = A† (7.5)
– Une matrice est orthogonale si son inverse est ´egale `a sa transpos´ee
[Link] = AT .A = 1 (7.6)
– Une matrice est unitaire si sa matrice adjointe est ´egale a` son inverse
A.A† = A†.A = 1 (7.7)
Pour les matrices `a coefficients r´eels, il y a identit´e de d´efinition entre
matrice sym´etrique et Hermitienne, entre matrice orthogonale et unitaire.
Une matrice est dite normale si elle commute avec son adjointe.

7.2 Propri´et´es des matrices

Les valeurs propres d’une matrice Hermitienne sont toutes r´eelles. Parmi
les matrices Hermitiennes tr`es utilis´ees en Physique, il vient `a l’esprit la
repr´esen- tation matricielle de l’op´erateur de Schr¨odinger en m´ecanique
quantique1
Un corollaire de la propri´et´e pr´ec´edente est que les valeurs propres d’une
matrice r´eelle sym´etrique sont elles aussi toutes r´eelles.
Les vecteurs propres d’une matrice normale ne poss´edant que des
valeurs propres non d´eg´en´er´ees forment une base d’un espace vectoriel
de dimension N . Pour une matrice normale avec des valeurs propres d´eg´en
´er´ees, les vecteurs propres correspondant `a une valeur propre d´eg´en´er
´ee peuvent ˆetre remplac´es par une combinaison lin´eaire de ceux-ci.
Pour une matrice quelconque, l’ensemble des vecteurs propres ne constituent
pas n´ecessairement une base d’un espace de dimension N .
Pour une matrice non normale, les vecteurs propres ne sont pas orthogonaux.
On appelle vecteur `a droite les vecteurs tels que
[Link] = λixR (7.8)
i i
1
La repr´esentation de cet op´erateur correspond en g´en´eral a` une matrice de
dimension infinie, mais nous ne consid´erons ici que les syst`emes ou la repr´esentation
matricielle est possible dans un espace de dimension finie.

9
 7.2 Propri´et´es des matrices

ou` λi est la i`eme valeur propre. De mani`ere similaire, on appelle vecteurs propres
`a gauche, les vecteurs tels que
xL.A = λixL (7.9)
i i

Le transpos´e du vecteur `a gauche de A est le vecteur propre `a droite

de la transpos´ee de la mˆeme matrice. Si la matrice est sym´etrique, les
vecteurs propres
`a gauche sont les transpos´es des vecteurs propres `a droite2 .
Si la matrice est Hermitienne, les vecteurs propres `a gauche sont les
trans- pos´es des vecteurs propres conjugu´es `a droite.
Dans le cas d’une matrice non normale, on d´efinit la matrice XR comme
la matrice constitu´ee de colonnes form´ees par les vecteurs `a droite. On
introduit la matrice XL form´ee par les lignes des vecteurs `a gauche. On obtient
par d´efinition que
[Link] = [Link](λ1, . . . , λn) (7.10)
de mˆeme on
a XL.A = diag(λ1, . . . , λn).XL (7.11)
En multipliant l’´equation (7.10) par XL et l’´equation (7.11) par XR , on obtient

[Link](λ1, . . . , λn) = diag(λ1, . . . , λn)[Link] (7.12)

ce qui montre que la matrice diagonale form´ee par les valeurs propres de
A commute avec le produit [Link]. Sachant que les seules matrices qui
commutent avec une matrice diagonale constitu´ee d’´el´ements diff´erents
sont elles-mˆemes diagonales, on en d´eduit que chaque vecteur `a gauche
est orthogonal a` chaque vecteur `a droite et r´eciproquement. En
normalisant les vecteurs `a droite et `a gauche, on peut obtenir que la
matrice XL .XR soit ´egale `a l’identit´e.
Dans le cas ou` l’ensemble des vecteurs propres ne constitue une base
com- pl`ete, il est toujours possible de compl´eter cet ensemble afin d’avoir une
matrice telle que [Link] = 1.
Si la matrice A est inversible, on obtient en multipliant l’´equation (7.10) par
−1
X que
R
X−[Link] = diag(λ1, . . . , λn) (7.13)
R
Nous avons alors construit une matrice de transformation similaire de A
Rappelons la propri´et´e suivante : soit B une matrice telle que

B = P −1.A.P (7.14)

ou` P est une matrice inversible. On a

det(B − λ1) = det(P −1.A.P − λ.1) (7.15)

−1
= det(P .(A − λ.1)P ) (7.16)
= det(A − λ.1) (7.17)
2
Puisque le d´eterminant d’une matrice et de sa transpos´ee sont les mˆemes, les

9
valeurs propres de ces deux matrices sont identiques.

9
Analyse spectrale

Ainsi, on a montr´e que l’on peut construire une matrice de

transformation similaire ou` la matrice A devient diagonale dans cette
nouvelle base.
Pour une matrice r´eelle sym´etrique, la matrice de passage est une
matrice orthogonale.
La strat´egie g´en´erale pour d´eterminer les valeurs propres d’une
matrice consiste `a construire une suite de transformations de similarit´e
jusqu’`a l’ob- tention d’une matrice diagonale, ou plus simplement jusqu’`a
l’obtention d’une matrice tridiagonale `a partir de laquelle il est possible de
d´eterminer assez fa- cilement les valeurs propres.
Pour r´ealiser ces transformations, deux grandes classes de m´ethodes
sont disponibles : les m´ethodes directes et les m´ethodes it´eratives.
Les premi`eres consistent `a effectuer une suite de transformations
similaires et ne s’appliquent qu’aux matrices de taille relativement modeste, car
le temps de calcul croˆıt comme le cube de la dimension lin´eaire de la
matrice.
Pour les matrices de grande taille et g´en´eralement creuses, les m
´ethodes it´eratives sont plus adapt´ees. Dans la mesure ou` la plupart du
temps, le spectre complet d’une tr`es grande matrice n’est pas recherch´ee,
mais seulement une partie, les m´ethodes it´eratives peuvent r´ealiser plus
efficacement cette taˆche. Nous allons voir dans la suite de ce chapitre
quelques algorithmes de base, tout en ayant `a l’esprit qu’il existe une tr`es
vaste litt´erature sur ce sujet, et pour un probl`eme particulier, il est n
´ecessaire de commencer par analyser pr´ecis´ement le type de matrice dont
on souhaite obtenir le spectre, afin de choisir la m´ethode la plus adapt´ee
pour r´esoudre ce probl`eme.

7.3 M´ethodes directes

7.3.1 M´ethode de Jacobi
La m´ethode de Jacobi revient `a effectuer une suite de transformations
si- milaires orthogonales. Chaque transformation est une simple rotation
planaire qui permet d’annuler un ´el´ement de la matrice A initial.
La rotation ´el´ementaire Ppq est donn´e par la matrice
1 0 0 ... ... ... ... 0
0
0
1
...
0
1
...
...
... ... ... 0
... ... ... 0
 
0 ... c ... ... s ... 0
0 ... ... ... ... ... ... 0
Ppq =  (7.18)
0 ... ... ... 1....................0
0 ... ... ... . . . . . . . . . 0

0 ... s ... ... c ... 0
0 ... .—
.. ... ... ... 1 0
0 ... ... ... ... ... ... 1
avec la condition que
c2 + s2 = 1. (7.19)

9
Soit la matrice AJ telle que
AJ = P T .[Link] (7.20)
p

9
 7.3 M´ethodes directes

En notant les coefficients de la matrice aij , on obtient apr`es calculs que

aJrp = carp − sarq (7.21)

aJrq = carq + sarp (7.22)
a pp = c app + s aqq − 2csapq
J 2 2
(7.23)
aJqq = s2app + c2aqq + 2csapq (7.24)
a pq = (c − s )apq + cs(app − aqq)
J 2 2
(7.25)

avec r /= p et r /= q.
Si on annule le terme aJpq , en introduisant l’angle de rotation φ, (c = cos(φ),
s = sin(φ)) on a le rapport θ
c2 − s2
= (7.26)
θ 2sc
= cot(2φ) (7.27)
aqq − app
= apq (7.28)

Si on appelle t = s/c, on obtient en utilisant l’´equation

(7.26)

t2 + 2tθ − 1 = 0 (7.29)
La plus petite des racines correspond `a un angle de rotation inf´erieur `a
π/4 et donne la m´ethode la plus stable num´eriquement. Cette racine3 peut
s’exprimer sous la forme
sgn(θ)
√ t= (7.30)
θ
| | + θ2 + 1
En utilisant que c + s = 1 on obtient pour c que
2 2
”
1
c= √ (7.31)
1 + t2
et on a imm´ediatement s = tc. En imposant que le terme aJpq s’annule, on a
finalement les relations suivantes

aJpp = app − tapq (7.32)

a qq = aqq + tapq
J
(7.33)
a rp = arp − s(arq + τarq)
J
(7.34)
a rq = arq + s(arp + τarq)
J
(7.35)

avec τ d´efini par

s
τ = (7.36)
1+c
En calculant la somme S
Σ
S= |ars |2 (7.37)
3 r s de l’ordinateur, on choisit t = 1/2θ
Pour ´eviter les d´epassements de capacit´e

1
Analyse spectrale

on peut obtenir une estimation de la convergence de la m´ethode. Pour

une transformation similaire ´el´ementaire, on a

SJ = S − 2|apq|2 (7.38)

Ainsi la suite des transformations conduit `a faire d´ecroˆıtre la contribution des

´el´ements non diagonaux. Comme la transformation est orthogonale, la
somme des carr´es des ´el´ements de la matrice est conserv´ee, ce qui
revient `a ce que la somme des carr´es de la diagonale augmente de ce qui
a ´et´e perdu par les
´el´ements non diagonaux. Ainsi formellement, on peut choisir les ´el´ements
de la matrice A dans n’importe quel ordre et on obtient une m´ethode qui
converge vers une matrice diagonale. Au terme de cette proc´edure, on a

D = V T .A.V (7.39)

ou` D est une matrice diagonale contenant les diff´erentes valeurs propres et V
est une matrice contenant les vecteurs propres correspondants.

7.3.2 R´eduction de Householder

La m´ethode pr´ec´edente est tr`es couˆteuse en temps de calcul. Pour
r´eduire celui-ci, la proc´edure de Householder se propose de transformer
une matrice sym´etrique en une matrice tridiagonale par une s´erie de
transformations ortho- gonales suivantes.
Une matrice de Householder est d´efinie par la relation suivante

P = 1 − [Link] (7.40)

ou` w est un vecteur r´eel normalis´e, w T .w = w 2 = 1.

|
V´erifions que la matrice P est une matrice orthogonale

P 2 = (1 − [Link] ).(1 − [Link] ) (7.41)

= 1 − 4w.w + 4w.(w .w).w
T T T
(7.42)
=1 (7.43)

Cela implique que P = P −1 . En utilisant la d´efinition de P , on v´erifie

facilement que P T = P , et on a donc bien construit une transformation
orthogonale.
Nous allons maintenant appliquer `a P le vecteur x constitu´e de la
premi`ere colonne de A. Pour cela, on exprime la matrice P sous la forme
suivante
[Link]
P =1− (7.44)
H
avec H = |u|2/2. Si on choisit le vecteur u tel que

u = x ∓ |x|e1 (7.45)

ou` e1 est le vecteur colonne unitaire tel que seule la premi`ere composante
est non nulle et ´egale `a 1. On obtient facilement la valeur de H
1
H = 2(|x2| ± |x|x1) (7.46)

1
 7.3 M´ethodes directes

En appliquant P `a
x u
P.x =x − .(x ∓ |x|e 1)T .x (7.47)
H
=x − 2u.(|x|2 ∓ |x|x1)
2|x|2 ∓ 2|x|x (7.48)
1
=x − u (7.49)
= ± |x|e1 (7.50)
Cela montre que la matrice P annule tous les ´el´ements du vecteur x
hormis le premier.
La strat´egie pour construire les matrices de Householder est la suivante :
on choisit un vecteur x constitu´e des n 1 derniers ´el´ements de la premi`ere
−
colonne pour construire la matrice P1 . En cons´equence, on obtient la structure
suivante pour P1
 
1 0 0 ... ... 0
0
0 
P1 = 0  (7.51)
 (n−1)
.
. P1
 
0
0 

En appliquant la transformation orthogonale `a la matrice A, on obtient

AJ =P.A.P (7.52)
 
a11 k 0.....................0
k
0 
= 0 (7.53)
 
.
. 
0 

 T
ou` le nombre k est au signe pr`es 0 la norme du vecteur (a21 , . . . , an1 ) .
On choisit la seconde matrice de Householder avec un vecteur x qui constitu
´e avec les (n − 2) derniers ´el´ements de la seconde colonne
1 0 0 ... ... 0
0 1 0 ... ... 0
0 0 
P1 =  0 0

(7.54)
. .. . P2 (n

 
0 0
0 0
La tridiagonalisation de la partie sup´erieure de la matrice pr´ec´edente est pr
´eser- v´ee, on construit ainsi colonne par colonne une tridiagonalisation de la
matrice

1
Analyse spectrale

A. La proc´edure est compl`ete apr`es (n 2) transformations similaires de

−
Hou- seholder.
Pratiquement, pour ´eviter la multiplication de matrices, tr`es
couˆteuse nu- m´eriquement, on peut exprimer le produit P.A.P en
introduisant la notation suivante
A.u
p (7.55)
=H
En cons´equence, la premi`ere multiplication matricielle peut ˆetre ´ecrite comme
[Link]
A.P = A.(1 − ) (7.56)
H
= A − [Link] (7.57)

de mˆeme pour la
seconde

AJ = P.A.P = A − [Link] − [Link] + [Link] (7.58)

avec
uT .p
K= (7.59)
2H
En posant
q = p − Ku (7.60)
on a la matrice A qui s’exprime alors simplement
J

AJ = A − [Link] − [Link] (7.61)

7.3.3 Algorithme QL
En utilisant une suite de transformations de Householder, on peut ´ecrire
tout matrice r´eelle sous la forme

A = Q.R (7.62)

ou` Q est une matrice orthogonale et R une matrice triangulaire sup´erieure.

Pour obtenir une telle d´ecomposition, les matrices de Householder sont
construites dans ce cas de mani`ere `a ce que la premi`ere colonne ne
poss`ede que le premier
´el´ement non nul une fois la transformation effectu´ee, pour la deuxi`eme
colonne, on choisit le vecteur de la matrice de Householder pour que tous les ´el
´ements de la matrice transform´ee soient nuls sous la diagonale et ainsi de
suite jusqu’`a for- mer une matrice triangulaire sup´erieure. Le nombre de
matrices de Householder pour obtenir ce r´esultat est ´egal `a (n 1).
−
De mani`ere analogue, on peut montrer qu’il existe une d´ecomposition de
la forme
A = Q.L (7.63)
. ou` Q est une matrice orthogonale et L une matrice triangulaire inf´erieure.
Les transformations de Householder correspondantes consistent `a annuler,
colonne par colonne, les ´el´ements de la matrice transform´ee qui sont situ´es

1
au dessus de la diagonale.

1
 7.3 M´ethodes directes

Si on d´efinit la matrice AJ comme

AJ = L.Q (7.64)

Puisque Q est orthogonal, on en d´eduit de l’´equation (7.63) que

L = Q−1.A (7.65)

ce qui donne
AJ = QT .A.Q (7.66)

ce qui montre que AJ est une transformation orthogonale de A.

Pour des raisons de minimisation d’erreurs d’arrondi, il est pr´ef´erable
d’uti- liser la d´ecomposition QL au lieu de la d´ecomposition QR.
L’algorithme QL est d´efini par la suite suivante

As = [Link] (7.67)
As+1 = [Link] (7.68)

La m´ethode repose sur les bases suivantes : (i) Si A a des valeurs

propres toutes distinctes de valeur absolue λi , alors As tend vers une matrice
|
triangu- laire inf´erieure quand s . Les valeurs propres apparaissent sur la
→
diagonale par valeur absolue croissante (ii) Si A a une valeur propre d´eg´en´er
´ee de multipli- cit´e p,quand s , As tend vers une matrice triangulaire inf
→
´erieure, except´ee pour un bloc d’ordre p correspondant `a la valeur propre
d´eg´en´er´ee. Pour une matrice quelconque, une it´eration a un couˆt de
calcul proportionnel `a n3 , mais pour une matrice tridiagonale, ce couˆt est
lin´eaire avec n. On peut montrer que la convergence d´epend de la diff
´erence entre deux valeurs propres succes- sives. Quand deux valeurs propres
sont trop proches, il est quand mˆeme possible d’am´eliorer la convergence de
l’algorithme en d´eplac¸ant ces valeurs propres suc- cessives.

7.3.4 Factorisation de Schur

La factorization de Schur consiste `a r´e´ecrire une matrice carr´ee A
sous la forme suivante
– Si la matrice A est complexe

A = ZTZ † (7.69)

ou` Z est unitaire et T est une matrice triangulaire sup´erieure.

– Si la matrice A est r´eelle
A = ZTZ T (7.70)

ou` Z est orthogonale et T est une matrice quasi-triangulaire sup

´erieure, ce qui signifie que la diagonale est constitu´ee soit de blocs 1
× 1 soit de blocs 2 × 2.

1
Analyse spectrale

Les colonnes de Z sont appel´ees les vecteurs de Schur. Les valeurs propres
de A apparaissent sur la diagonale de T ; les valeurs propres complexes conjugu
´ees d’une matrice A r´eelle correspondent aux blocs 2x2 de la diagonale.
L’algorithme utilis´ee dans la biblioth`eque LAPACK, on commence par
trans- former la matrice A en la transformant en une matrice de Hessenberg, qui
est une matrice triangulaire sup´erieure bord´ee une ligne d’´el´em´ents
nuls sous la diagonale.
– Si la matrice A est complexe

A = QHQ† (7.71)

ou` Q est unitaire et H est une matrice de Hessenberg.

– Si la matrice A est r´eelle
A = QT QT (7.72)
ou` Q est orthogonale et H est une matrice de Hessenberg.
Dans une deuxi`eme ´etape, on transforme la matrice de Hessenberg en une ma-
trice de Schur.

7.4 M´ethode it´eratives

Ces m´ethodes sont principalement appliqu´ees `a des matrices de grande
taille et largement creuses. Les calculs croissent dans ce cas lin´eairement
avec n. Pour toute m´ethode it´erative, il convient de s’assurer que la
convergence est suffisamment rapide pour que le temps de calcul ne soit pas
consomm´e sans que la recherche d’une solution ne soit r´eellement effectu´ee.

7.4.1 M´ethodes des puissances

Le principe de cette m´ethode est tr`es simple, car il repose sur le fait
qu’en appliquant un grand nombre de fois la matrice sur un vecteur initial
quelconque, les vecteurs successifs vont prendre une direction qui se
rapproche du vecteur propre de la plus grande valeur propre (en valeur
absolue). Le principe it´eratif de cette m´ethode est la suivante : soit x0 un
vecteur initial quelconque, et A la matrice dont on cherche `a d´eterminer la
plus grande valeur propre. On effectue l’op´eration suivante
x1 = A.x0 (7.73)
Si on d´esigne α comme l’angle entre ces deux vecteurs, on a
x1.x0
cos(α) =
|x1|.|x0| (7.74)

Si x0 n’est pas perpendiculaire au vecteur recherch´e (ce qui est rarement

le cas pour un vecteur choisi au d´epart al´eatoirement), le cosinus de l’angle
entre x0 et x1 est diff´erent de z´ero. En appliquant `a nouveau la matrice
A sur le vecteur x1 on cr´ee un vecteur x2 et on calcule l’angle entre x1 et
x2 qui est inf´erieur au pr´ec´edent en valeur absolue. On continue jusqu’`a ce
que l’angle entre deux vecteurs successifs devienne plus petit qu’une nombre ϵ
choisi initialement.
1
 7.4 M´ethode it´eratives

On en d´eduit alors le vecteur propre recherch´e et donn´e par xn , ainsi

que la valeur propre correspondante. La convergence de cette m´ethode
varie comme (λ1/λ2) ce qui peut devenir assez lent quand les valeurs propres
deviennent quasi-d´eg´en´er´ees. Cette m´ethode n’est pas tr`es efficace,
mais poss`ede le m´erite de s’´ecrire rapidement.

7.4.2 M´ethode de Lanczo`s

La m´ethode de Lanczo`s consiste `a la fois `a calculer les puissances
successives de A, en s’inspirant de la m´ethode pr´ec´edente, mais de
mani`ere bien plus effi- cace en construisant un ensemble de vecteurs
orthogonaux. Par simplicit´e, nous allons voir la m´ethode pour des matrices
hermitiennes, mais la m´ethode peut ˆetre ´etendue pour des matrices plus g
´en´erales, en particulier quand les vecteurs
`a gauche diff`erent des vecteurs `a droite.
Soit un vecteur de d´epart normalis´e u0 : ce vecteur est convenablement
choisi c´est `a dire que sa projection sur la base de la valeur propre `a d
´eterminer est non nulle. On construit les vecteurs successifs de la base dite de
Krylov `a partir de la relation

β2u1 = A.u0 − α1u0 (7.75)

o β2 et α1 sont des constantes d´etermin´ees de la mani`ere suivante : α1 est

u
`
choisi de mani`ere `a ce que le vecteur u1 soit orthogonal au vecteur u0

α1 = uT .A.u0 (7.76)
0

et β2 est d´etermin´e de mani`ere `a ce que le vecteur u1 soit normalis´e.

q
β2 = uT .A2.u0 − α2
(7.77)
0 1

Pour le vecteur suivant u2, on utilise la relation :

β3u2 = A.u1 − α2u1 − β2u0 (7.78)

On note tout d’abord qu’avec cette construction, le vecteur u2 est orthogonal au

vecteur u0 . De mani`ere similaire, on impose que u2 soit un vecteur
orthogonal
`a u1 ce qui conduit `a la relation

α2 = uT .A.u1 (7.79)
1

et la normalisation de u2 est impos´ee en choisissant β3 comme

q
1 2
β3 = uT .A2.u1 − α2 − β2

1
Par it´eration, on obtient pour le i`eme vecteur

βi+1ui = [Link]−1 − αiui−1 − βiui−2 (7.81)

1
Analyse spectrale

On v´erifie facilement que tous les vecteurs uj avec j < i 2 sont

−
orthogonaux avec le vecteur ui+1 en raison de la construction dans un sous
espace ortho- norm´e sur ses vecteurs. Les valeurs αi et βi+1 sont d´etermin
´ees par les relations suivantes
αi =uTi−[Link]+1 (7.82)
q
βi+1 = uT .[Link]−1 − α2 −i βi (7.83)
i

Dans cette base, la matrice AJ est tridiagonale et est donn´ee

α1 
β2 β2 0 ... ... ... 0
0 α2 β 3 0 ... ... 0
A =
J
0 β 3 α3 β 4 0 ... 0 (7.84)
 0 ... ... ... 0
... ... ... ... ... 0
n−1 αn−1 β 
0 ... ... 0
β
0 
... ... ... 0 βn α
Il est possible d´utiliser alors une m´ethode de type QL pour obtenir
rapidement les valeurs propres de la matrice AJ , c´est `a dire aussi la matrice
A.
Quelques remarques : en construisant progressivement une base de vecteurs
de plus en plus grande, on voit que la valeur propre recherch´ee peut ˆetre estim
´ee en utilisant les sous espaces successifs. On peut donc arrˆeter l´it
´eration quand la diff´erence entre deux estimations successives de la valeur
propre est devenue suffisamment petite.
Dans le cas ou` la matrice A est une repr´esentation tronqu´ee d´un op
´erateur dont la base propre est infinie (par exemple, un op´erateur de
Schro¨ringer), le proc´ed´e it´eratif conduit progressivement `a des valeurs
de βi de plus en en plus petites. On peut consid´erer que la base de
Krylov est compl`ete quand βn est devenu inf´erieure en valeur absolue `a
une valeur ϵ choisie `a l´avance (g´en´eralement 1012 ).

1
Chapitre 8

Equations int´egrales

Contenu
8.1 Introduction....................................................................... 69
8.2 Equation de Fredholm......................................................69
8.2.1 Equation de premi`ere esp`ece ............................................69
8.2.2 Equation de seconde esp`ece ................................................70
8.3 Equation de Volterra.........................................................71
8.3.1 Equation de premi`ere esp`ece ............................................71
8.3.2 Equation de seconde esp`ece ................................................71
8.4 Conclusion.......................................................................... 71

8.1 Introduction
Les ´equations int´egrales sont a priori moins simples `a r´esoudre que les
´equa- tions alg´ebriques ou les ´equations diff´erentielles. Nous allons voir
dans ce cha- pitre que pour des ´equations int´egrales lin´eaires, une fois r
´ealis´ee la discr´etisa- tion de ces ´equations, on se ram`ene au probl`eme
de la recherche de solutions d’un syst`eme lin´eaire que nous avons vu
chapitre 6.

8.2 Equation de Fredholm

8.2.1 Equation de premi`ere esp`ece
L’´equation de Fredholm inhomog`ene de premi`ere esp`ece est d´efinie
par la relation suivante
∫b
K(t, s)f (s)ds = g(t) (8.1)
a

ou` f (t) est la fonction inconnue que l’on souhaite d´eterminer. g(t) est le
terme de source et K(t, s) est appel´e le noyau.
En notant gi = g(ti ), Kij = K(si , tj ) et fj = f (tj ) ou` i est un indice
variant de 1 `a N et j un indice variant de 1 `a M (M peut ˆetre diff´erent de N ).
1
Equations int´egrales

L’´equation (8.1) se r´e´ecrit alors comme

M
Σ
Kijfj = gi (8.2)
j=1

Soit encore sous une forme matricielle

K.f = g (8.3)

Formellement, si le noyau K n’est pas singulier, la solution existe, est unique

et est donn´ee par la relation
f = K−1.g (8.4)

8.2.2 Equation de seconde esp`ece

L’´equation de Fredholm inhomog`ene de deuxi`eme esp`ece est d´efinie
par la relation suivante
∫
λf (t) = b
K(t, s)f (s)ds + g(t) (8.5)
a

ou` f (t) est la fonction inconnue que l’on souhaite d´eterminer. g(t) est le terme
de source, K(t, s) est appel´e le noyau et λ est un scalaire introduit par
commodit´e pour la suite de cette ´equation.
Suivant le mˆeme principe que celui d´efini ci-dessus, une fois discr
´etis´ee, l’´equation (8.5) se r´eexprime comme

M
Σ
λfi = Kijfj + gi (8.6)
j=1

Soit encore sous une forme matricielle

(K − λ1).f = −g (8.7)

Si la fonction g est nulle, le probl`eme se ram`ene `a la d´etermination des

valeurs propres de la matrice K, et on parle d’´equation de Fredholm
homog`ene. Dans le cas ou` g est diff´erent de z´ero, la solution existe et
est unique si λ n’est pas l’une des valeurs propres du noyau, sinon la matrice
−
`a inverser K λ.1 devient singuli`ere. Si cette derni`ere est inversible, on a
formellement la solution comme

f = (λ.1 − K)−1.g (8.8)

La r´esolution num´erique des ´equations de Fredholm homog`ene de

premi`ere esp`ece est g´en´eralement d´elicate car le noyau correspond
souvent `a une matrice presque non inversible (mal conditionn´ee), c’est-a`-
dire avec un d´eterminant voi- sin de z´ero. Corr´elativement, si λ est

1
suffisament diff´erent de z´ero, la solution des ´equations de Fredholm de
seconde esp`ece est relativement simple `a obtenir.

1
 8.3 Equation de Volterra

8.3 Equation de Volterra

Les ´equations de Volterra sont des cas particuliers de ceux de Fredholm
dans lesquelles le noyau K est tel que

K(t, s) = 0 pour s > t (8.9)

8.3.1 Equation de premi`ere esp`ece

L’´equation de Volterra homog`ene de premi`ere esp`ece est d´efinie par
la rela- tion suivante
∫t
g(t) =
K(t, s)f (s)ds (8.10)
a
ou` f (t) est la fonction inconnue que l’on souhaite d´eterminer. g(t) est le
terme de source et K(t, s) est appel´e le noyau.
La r´e´ecriture matricielle de l’´equation de Volterra est identique formellement
`a celle de Fredholm, mais avec la particularit´e que la matrice associ´ee au
noyau K est une matrice triangulaire inf´erieure. Ce type d’´equations lin
´eaires, comme nous l’avons vu au chapitre 6, peut ˆetre r´esolu par une m
´ethode de substitution. Alors que les ´equations de Fredholm de premi`ere
esp`ece sont g´en´eralement mal conditionn´ees, les ´equations de Volterra ne
le sont pas.

8.3.2 Equation de seconde esp`ece

De mani`ere similaire, l’´equation de Volterra de premi`ere esp`ece
inhomog`ene s’´ecrit comme
∫t
λf (t) =
K(t, s)f (s)ds + g(t) (8.11)
a
De mani`ere identique, la repr´esentation matricielle de l’´equation de Volterra
est identique `a celle correspondante de Fredholm, seule la structure du
noyau K est diff´erente, puisque comme pour toutes les ´equations de Volterra
lin´eaires, la matrice K triangulaire inf´erieure.

8.4 Conclusion
L’existence de structure de noyaux presque singuliers n´ecessite d’utiliser
des m´ethodes plus complexes qui d´epassent largement le cadre de ce
cours intro- ductif aux m´ethodes num´eriques.

1
Equations int´egrales

1
Chapitre 9

Equations aux d´eriv´ees

partielles

Contenu
9.1 Introduction....................................................................... 73
9.2 Equations avec conditions aux fronti`eres .......................76
9.2.1 Introduction............................................................................76
9.2.2 Diff´erences finies ..................................................................76
9.2.3 M´ethodes matricielles ..........................................................77
9.2.4 M´ethodes de relaxation .......................................................77
9.2.5 M´ethodes de Fourier ...........................................................78
9.3 Equations avec conditions initiales.................................80
9.3.1 Equations `a flux conservatif ..............................................80
9.3.2 Une approche na¨ıve..............................................................81
9.3.3 Crit`ere de Stabilit´e de Von Neumann ............................81
9.3.4 M´ethode de Lax ...................................................................82
9.4 Conclusion.......................................................................... 83

9.1 Introduction
Les ´equations aux d´eriv´ees partielles interviennent dans de nombreux
do- maines de physique, qui comprennent les probl`emes de diffusion, les ph
´enom`enes de propagation, ainsi que le domaine de la m´ecanique des
fluides d´ecrite par les ´equations hydrodynamiques comme celles de Navier-
Stokes et l’´equation de Schr¨odinger d´ependante du temps pour la m
´ecanique quantique. Ces ´equations diff´erentielles n’ont g´en´eralement pas
de solutions analytiques et une r´esolution num´erique de ces ´equations est
alors n´ecessaire.
Une ´equation aux d´eriv´ees partielles est une relation liant une fonction
de n variables `a ses d´eriv´ees partielles. L’ordre de l’´equation est donn´e
par l’ordre le plus ´el´ev´e des d´eriv´ees partielles apparaissant dans
l’´equation. La forme g´en´erale d’une ´equation aux d´eriv´ees partielles lin
´eaires est

1
L[f (x)] = g(x) (9.1)

1
Equations aux d´eriv´ees partielles

4
Conditions
aux
limites
3
t

Conditions initiales
0
0 1 2 3 4 5 6 7 8 9
-
x

FIG. 9.1 – Sch´ema illustrant le principe de la r´esolution d’une ´equation

aux d´eriv´ees partielles avec des conditions initiales et des conditions aux
limites : La propagation de la solution le long de l’axe temporel (axe vertical) se
fait par “tranche horizontale”.

ou` un x est un vecteur de composante (x1 , x2 , . . . , xn ) et ou` L est un op

´erateur d´efini par la relation
n n
Σ Σ
L = p0 (x) + pi (x)∂i + pij (x)∂i ∂j + . . . (9.2)

Si g(x) = 0, on dit que l’´equation

i est homog`ene.
i,j

Dans ce cours d’introduction, nous allons nous limiter aux ´equations

aux d´eriv´ees partielles lin´eaires. Une premi`ere classification non
exhaustive fait ap- paraˆıtre trois types d’´equations : ´equations
hyperboliques, paraboliques et ellip- tiques. Cette classification provient de
l’analyse des d´eriv´ees partielles d’ordre le plus ´el´ev´e. Pour les
´equations d’ordre deux, une solution unique existe si la relation appel´e´e
´equation aux caract´eristiques
√ 2
— p11p22
12
p /= 0 (9.3)
p11
±
est vraie. Si la valeur du membre de gauche de l’´equation (9.3) est r´eelle, on
dit que l’on a une ´equation hyperbolique. Si la valeur est complexe, on dit que
l’on a une ´equation elliptique et si elle est nulle, on a une ´equation
parabolique.
On peut illustrer ces trois cat´egories par les exemples physiques suivants :
1
 9.1 Introduction

6
Conditions aux bords

y
3

0
0 1 2 3 4 5 6 7 8 9
- x

FIG. 9.2 – Sch´ema illustrant le principe de la r´esolution d’une ´equation

aux d´eriv´ees partielles avec des conditions aux bords : l´ensemble des
points de la grille doit ˆetre gard´e tout au long du calcul.

Equation hyperbolique. La propagation des ondes dans un milieu continu

(unidimensionnel) ob´eit `a l’´equation suivante

∂2u(x, t)
∂t2 = 2 ∂2u(x, t)
c (9.4)
∂x2
ou` u(x, t) peut designer´ une variable scalaire comme la densit´e locale
ρ(x, t) ou une variable vectorielle comme la vitesse locale v(x, t).
Equation parabolique. La diffusion de particules dans un milieu continu
ob´eit `a l’´equation suivante

∂ρ(x, t)
∂t = D ∂2ρ(x, t)
(9.5)
∂x2
ou` D est la constante de diffusion et ρ(x, t) la densit´e locale
instantan´ee des particules diffusantes.
Equation elliptique. En dimensions deux, le potentiel coulombien induit
par une densit´e de charges ´electriques satisfait l’´equation suivante
(avec des unit´es choisies de mani`ere appropri´ee)

−ρ(x, y) =
∂2V (x, y, t) ∂2V (x, y, t)
∂x2 + ∂y2 (9.6)

1
ou` ρ(x, y) est la densit´e locale de charge et V (x, y) le potentiel
´electrosta- tique `a d´eterminer.

1
Equations aux d´eriv´ees partielles

La r´esolution num´erique repose sur une distinction importante entre

ces type d’´equations : Les deux premi`eres ´equations correspondent `a une
´evolution temporelle tandis que la troisi`eme ´equation correspond `a un
probl`eme statique. Dans le premier cas, la solution est donn´ee `a partir d’une
fonction (ou condition initiale) que l’on doit propager le long de l’axe du
temps ; en compl´ement on doit tenir compte des conditions au limites du
probl`eme. En discr´etisant la r´esolution, on obtient un sch´ema du type de
celui illustr´e par la figure 9.1.
En effet pour une ´equation d´ependante du temps du type ´equation
de dif- fusion, la r´esolution de l’´equation se fait pour un temps donn´e en
fonction des valeurs de la fonction `a l’instant ant´erieur pr´ec´edent (ou
de quelques instants ant´erieurs pr´ec´edents). Il n’est donc pas n´ecessaire
de garder en m´emoire la solution `a tous les temps.
Pour l’´equation de Poisson, la solution est unique une fois donn´ee la
fonc- tion ρ(x, y). La solution ne pouvant ˆetre obtenue num´eriquement
directement
`a partir de l’int´egration `a partir d’une fronti`ere du syst`eme, il est n
´ecessaire de proc´eder `a des r´eajustements, ce qui implique de garder en m
´emoire la totalit´e de valeurs de la fonction V `a d´eterminer tout au long de
la r´esolution num´erique (voir figure 9.2).

9.2 Equations avec conditions aux fronti`eres

9.2.1 Introduction
Le prototype de ces ´equations est la r´esolution de l’´equation de Poisson. Ces
´equations sont a priori plus simples `a r´esoudre car le probl`eme de la
stabilit´e de l’algorithme se pose moins fr´equemment que dans pour des
´equations avec conditions initiales, que nous allons discuter dans la section
suivante.
Avant de proc´eder `a la r´esolution num´erique, il est important de bien
distin- guer le type de conditions aux fronti`eres du probl`eme : condition
de Dirichlet, qui sp´ecifie les valeurs des points `a la fronti`ere ; condition
de Neumann qui sp´ecifie les valeurs des gradients normaux `a la
fronti`ere, voire des conditions mixtes.

9.2.2 Diff´erences finies

La premi`ere ´etape consiste `a choisir un maillage r´egulier de l’espace
(que nous choisissons bidimensionnel dans notre exemple). Soit la fonction de
deux variables V (x, y)

xi = x0 + j∆ j = 0, 1, ..., J (9.7)
yk = x0 + k∆ k = 0, 1, ..., K (9.8)

ou` ∆ est le pas de la grille. On note Vjk = V (xj , yk ). On approxime le

laplacien par la formule de diff´erence finie suivante

1
Vj+1,k + Vj−1,k − 2Vj,k Vj,k+1 + Vj,k−1 − 2Vj,k
+ = −ρ (9.9)
j
∆ ∆

1
 9.2 Equations avec conditions aux fronti`eres

que l’on peut ´ecrire simplement

Vj+1,k + Vj−1,k + Vj,k+1 + Vj,k−1 − 4Vj,k = −∆2 ρjk (9.10)

Ce syst`eme d’´equations lin´eaires peut s’exprimer sous une forme matricielle

en utilisant l’indexation suivante

i = j(K + 1) + k j = 0, 1, ..., J et k = 0, 1, ..., K (9.11)

L’´equation (9.9) est correcte pour l’ensemble des points int´erieurs au

rectangle. Pour les points aux fronti`eres, les valeurs de V ou de ces d´eriv
´ees sont don- n´ees par le probl`eme lui-mˆeme. On obtient finalement la
structure matricielle suivante.
A.Y = b (9.12)
La matrice A est une matrice tridiagonale avec des franges, le vecteur Y est un
vecteur `a JK composantes (Yi Vjk ) et le vecteur b `a JK composantes
≡
(bi
ρ ) Ce type de repr´esentation est similaire pour des ´equations elliptiques
− jk
du second ordre
∂V
∂2V 2
+ c(x, y) ∂ V
a(x, y) 2 +b(x,
∂x y) ∂x ∂y2 2
∂V ∂V
+ d(x, y) + e(x, y) + f (x, y)V = g(x, y) (9.13)
∂y ∂x∂y
En choisissant une classification tr`es rapide, on peut dire qu’il existe
trois types de m´ethodes pour r´esoudre ce type d’´equations : les m´ethodes
de relaxa- tion, les m´ethodes de Fourier et les m´ethodes matricielles directes.

9.2.3 M´ethodes matricielles

Les m´ethodes matricielles ou m´ethodes directes consistent `a d
´eterminer les valeurs de la fonction V en calculant la matrice inverse de
A. La difficult´e de cette m´ethode repose sur la taille consid´erable de la
matrice que l’on doit inverser. En effet, si on consid`ere un r´eseau
bidimensionnel dont la taille est
100 100, et compte tenu de l’´equation (9.11), la matrice A est alors une
×
matrice carr´ee de taille
× 10000 10000, contenant 108 ´el´ements. Il est
donc n´ecessaire que cette matrice soit tr`es creuse et de structure simple
afin que le calcul de son inverse soit effectu´e avec un algorithme rapide.
Pratiquement, la matrice inverse est facile `a obtenir quand les ´el
´ements de matrice ne d´ependent pas des abscisses du r´eseau (cas d’une
´equation elliptique
`a coefficients constants).

9.2.4 M´ethodes de relaxation

1
 Les m´ethodes de relaxation reposent sur le sch´ema it´eratif suivant.
On d´e- compose la matrice A comme

A=E−F (9.14)

1
Equations aux d´eriv´ees partielles

ou` E est une matrice facilement inversible et F le reste. On peut ´ecrire

alors que
E.u = F.u + b (9.15)
En choisissant un vecteur u(0) , on proc`ede `a la succession de calculs suivants
E.u(r) = F.u(r−1) + b (9.16)
La convergence de la m´ethode est obtenue quand la diff´erence entre deux
valeurs de la fonction est inf´erieure `a ϵ, ǁu(r) − u(r−1) ǁ < ϵ. La norme ǁuǁ
Σ
d´esigne par
ij |uij |

9.2.5 M´ethodes de Fourier

Principe
On exprime la transform´ee de Fourier discr`ete de V dans les deux directions
Ox et Oy.
K−1 L−1
Σ Σ 1V −2iπkm/K
V =
k me e−2iπln/L (9.17)
K
m=0
Une expression analogue pour la fonction
ˆ ρ peut ˆetre obtenue
K−1 L−1
1 ΣΣ
ρ = −2iπkm/K −2iπln/L
(9.18)
me
k
e
K
ρˆ
m=0
En prenant la transform´ee de Fourier discr`ete de l’´equation (9.10), les
compo- santes de Fourier de V sont reli´ees `a celles de ρ par la relation

Vˆmn e2iπm/K + e−2iπm/K + e2iπn/L + e−2iπn/L − 4 = −ρˆmn ∆2 (9.19)

ce qui peut se r´e´ecrire plus simplement pour si n et m diff´erents de z´ero

comme
— mn
=
Vˆmn ρˆ ∆2
(9.20)
2πm 2πn
2 cos + cos −2
K L
La relation alg´ebrique entre les composantes de Fourier de V et de ρ
indique que le calcul dans l’espace de Fourier est ´el´ementaire et revient
`a multiplier individuellement chaque composante de ρ par un scalaire.
La m´ethode de r´esolution num´erique est donc la suivante :
– Calculer les composantes de Fourier de ρ : ρˆmn
– Calculer les composantes de Fourier de V en utilisant la relation (9.20) :
Vˆmn
– Calculer la transform´ee de Fourier inverse de Vˆ pour obtenir V .
Noter que la transform´ee de Fourier obtenue donne les relations suivantes
ujk = uj+J,k = uj,k+K (9.21)
ce qui en d’autres termes correspond `a une solution satisfaisant des
1
conditions aux limites p´eriodiques. A partir de l’´equation (9.20), on voit
que l’on peut choisir Vˆ00 = 0, mais qu’il est n´ec´essaire que ρˆ00 = 0. Si
cette derni`ere relation
n’est pas satisfaite, il faut modifier la solution choisie pour V comme nous allons
le voir ci-dessous.

1
 9.2 Equations avec conditions aux fronti`eres

Condition de Dirichlet
Si on impose des conditions de Dirichlet aux bords du rectangle, c’est-a`-
dire u = 0 pour j = 0 et j = J et pour k = 0 et k = K, la transform´ee
de Fourier adapt´ee est alors celle en sinus.
Σ K−1
Σ L−1
2 πkm πln
Vkl = ˆV sin (9.22)
mn
KL m=0 n=0 K L
sin
Une expression similaire pour la fonction ρ peut ˆetre ´ecrite. Un calcul
simple montre que la relation (9.20) relie les composantes de Fourier en
sinus des deux fonctions. et la m´ethode de calcul donn´ee dans le
paragraphe pr´ec´edent s’applique `a nouveau `a condition de consid´erer
les transform´ees de Fourier en sinus.

Condition de Neumann
Si on impose des conditions de Neumann aux bords du rectangle, c’est-
a`- dire u = 0 pour j = 0 et j = J et pour k = 0 et k = K, la
∇
transform´ee de Fourier adapt´ee est alors celle en cosinus.
Σ K−1
Σ L−1
2 πkm πln
Vkl = ˆV cos (9.23)
KL K L
mncos
m=0 n=0

ainsi que pour la fonction ρ. A nouveau, les composantes de Fourier sont reli
´ees par la relation (9.20). Le mˆeme sch´ema donn´e ci-dessus s’applique
alors.

Conditions complexes
Ceci repr´esente le cas g´en´eral ou` par exemple la fonction u(x, y)
s’annule sur l’une des fronti`eres (j = 0), mais vaut u = f (y) sur la fronti`ere
oppos´ee (j = J). La m´ethode consiste alors `a ´ecrire que la solution de ce
probl`eme est la somme de deux contributions : celle de la solution de
l’´equation pour une condition de Dirichlet et celle de la solution qui est
nulle `a l’int´erieur du rectangle. En
´ecrivant la solution sous la forme
u = uJ + uB (9.24)
avec uB = 0 sauf uB = f `a la fronti`ere, on obtient

∇2uJ = −∇2uB + ρ (9.25)

ce qui montre que la fonction satisfait alors une condition de Dirichlet. Le calcul
de uJ est effectu´e en ajoutant un terme de source `a ρ. Pratiquement, cela
ajoute une contribution sur la ligne situ´ee juste avant la fronti`ere.
De mani`ere analogue, si la condition aux fronti`eres est donn´ee par
un gra- dient normal non nul, on d´ecompose la solution d’une premi`ere
contribution correspondant `a la solution du probl`eme de Neumann et

1
d’une seconde contri- bution nulle en dehors de la fronti`ere. En utilisant une
d´ecomposition analogue
`a celle introduite ci-dessus, cela revient `a r´esoudre le probl`eme de
Neumann en modifiant le terme de source.

1
Equations aux d´eriv´ees partielles

9.3 Equations avec conditions initiales

9.3.1 Equations `a flux conservatif
Les ´equations `a flux conservatif dans un espace unidimensionel peuvent
s’´ecrire sous la forme
∂u ∂F (u) (9.26)
=−
∂ ∂
ou` u est une fonction scalaire et F un vecteur, appel´e flux conserv´e, on
peut citer l´exemple de l’´equation de la conservation de la masse
∂ρ
= .(ρv) (9.27)
∂t −
Si de plus on a une relation constitutive, comme la loi de Fick, F = D (ρ),
—
on obtient l’´equation parabolique de la diffusion.
Si u est un vecteur, l’´equation de propagation des ondes (sonores ou
´elec- tromagn´etiques, par exemple) se ram`ene `a ce type d’´equation.
En effet soit l’´equation des ondes suivante,

∂ 2v 2
∂t2 = c2 ∂ v (9.28)
∂x2
on introduit les fonctions auxiliaires r et s
∂v
r c (9.29)
= ∂x
∂v
s (9.30)
= ∂t
et l’´equation de propagation s’exprime alors par un ensemble de deux ´equations
aux d´eriv´ees partielles du premier ordre

∂r ∂s
=c (9.31)
∂t ∂x
∂s ∂r (9.32)
∂t = c ∂x
Si on consid`ere que r et s sont les deux composantes du vecteur u, l’´equation
de propagation est alors donn´ee par l’´equation (9.26) avec la relation
matricielle suivante
0 −v
F (u) = .u (9.33)
−v 0
Ainsi les ´equations `a flux conserv´e comprennent `a la fois les
´equations aux d´eriv´ees partielles paraboliques et hyperboliques.
Par souci de simplicit´e, nous allons consid´erer par la suite l’´equation `a
flux conservatif suivante
∂u ∂u (9.34)
= −c
∂ ∂
La solution de cette ´equation est bien ´evidemment connue u = f (x ct),
mais nous allons chercher `a construire une proc´edure num´erique − pour d
´eterminer la solution.
1
 9.3 Equations avec conditions initiales

9.3.2 Une approche na¨ıve

Une premi`ere ´etape consiste `a discr´etiser l’espace et le temps sur un r
´eseau r´egulier

xj = x0 + j∆x j = 0, 1, . . . J (9.35)
tn = t0 + n∆t n = 0, 1, . . . N (9.36)

ou` ∆x et ∆t sont les pas d’espace et de temps respectivement.

La difficult´e r´eside maintenant dans le choix de l’algorithme d’int´egration.
La m´ethode la plus simple consiste `a utiliser un algorithme de type Euler

∂u .
un+1
j − un
= j
+ O(∆ )
∂
j
(9.37)
t
Mˆeme si la pr´ecision de. l’algorithme n’est pas tr`es importante, cela
permet d’exprimer la solution `a l’instant n + 1 uniquement en fonction de
l’instant n.
Pour int´egrer spatialement, on utilise un algorithme du second ordre

∂u
. = n n
j− 1
∂x j,n uj+1 —
. + O(∆x2) (9.38)
2
ce qui donne finalement

un+1 − un un − un
j j+1 j−1 (9.39)
= −c 2∆
∆t
soit encore x
n n
un+1 = un − c∆t uj+1 j−1
— (9.40)
j j
La d´etermination de u n+1
se fait en fonction des
2 trois points : u , u
n n
et un .
j j j−1
j+1
On parle de m´ethode explicite. Malheureusement, ce sch´ema tr`es simple
ne fonctionne pas en ce sens qu’il ne permet pas d’obtenir une solution correcte
de l’´equation aux d´eriv´ees partielles. En d’autres termes, on dit qu’il est
instable.

9.3.3 Crit`ere de Stabilit´e de Von Neumann

Pour faire l’analyse de la stabilit´e des algorithmes, on suppose que les
coef- ficients de l’´equation aux diff´erences varient tr`es peu dans l’espace et
le temps. Les modes propres solutions de cette ´equation peuvent alors
s’´ecrire sous la forme
un = ξ(k)n eikj∆x (9.41)
j
ou` k est un vecteur d’onde r´eel et ξ(k) une fonction complexe qui d´epend
de k. On voit facilement que l’int´egration temporelle conduit `a une
progression g´eo- m´etrique pour les valeurs de un . Afin que la solution soit
j

1
stable, il est n´ecessaire que les modes exponentiellement divergents
n’existent pas. Cela revient `a dire que l’algorithme est instable d`es qu’il
existe un vecteur d’onde k0 tel que

|ξ(k0)| > 1 (9.42)

1
Equations aux d´eriv´ees partielles

5 5
(a) (b)
Instable Stable
4.5 4.5

4 4

3.5 3.5

3 3
t

t
2.5 2.5

2 2

1.5 1.5

1
1 2 3 4 5 6 7 8 1
x 1 2 3 4 5 6 7 8
x

FIG. 9.3 – Sch´ema illustrant le principe de la r´esolution d’une ´equation

aux d´eriv´ees partielles avec la m´ethode de Lax : La zone hachur´ee
indique la zone de d´ependance (a) M´ethode instable (b) M´ethode stable

En ins´erant l’´equation (9.41) dans l’´equation (9.39) et on obtient

c∆t
ξ(k) = 1 − i
sin(k∆x) (9.43)
∆x
Le module de ξ(k) est strictement sup´erieur `a 1 pour tous les vecteurs
d’onde, hormis ceux ou` k∆x = pπ ou` p est un entier. Ceci permet de
comprendre la nature catastrophique de la m´ethode pr´ec´edente.

9.3.4 M´ethode de Lax

Pour corriger le grave probl`eme d’instabilit´e de la m´ethode pr´ec´edente,
Lax a propos´e la modification suivante : le terme un provenant de la d´eriv´ee
j
temporelle est remplac´e par la moyenne sur les deux points adjacents :

1
un → (un
+ un ) (9.44)
j j
2 j

ce qui donne la relation it´erative

suivante
un − un
1
n + un (9.45)
un+1 = j+1 j−1 ) +
j− 1 ∆t
j+1
2 (
j 2∆x
En ins´erant la solution (9.41) dans l’´equation (9.45), on obtient une expression
pour ξ(k)
1
 c∆t
ξ(k) = cos(k∆x) − i sin(k∆x) (9.46)
∆x

1
 9.4 Conclusion

La condition de stabilit´e |ξ(k)| < 1 ´equivaut `a choisir

c∆t
≤1 (9.47)
∆x
Cette condition, dite de Courant est appel´ee aussi condition de stabilit´e
de Courant-Friedrichs-Lewy. Il est facile de comprendre que la discr´etisation
spatio- temporelle doit se faire `a une vitesse qui est donn´ee par ∆x inf´erieur
`a la vitesse de propagation du ph´enom`ene c.∆ Pour comprendre
graphiquement la raison pour laquelle le rapport du pas en temps sur celui de
la position ne doit pas exc´eder la vitesse de propagation, consid´erons les
figures suivantes : les zones hachur´ees correspondent `a la solution donn´ee
par la solution exacte de l’´equation de propagation : si la valeur du point `a
l’´etape n + 1 d´epend des valeurs des points situ´ees dans la zone hachur´ee,
l’algorithme est instable. Inversement dans la figure de droite, le sommet de
cette zone d´epend de points situ´es `a l’ext´erieur de la zone hachur´ee et
l’algorithme est stable.

9.4 Conclusion
Les difficult´es des m´ethodes num´eriques pour la r´esolution des
´equations d´eriv´ees partielles ont des origines diverses : elles proviennent soit
des subtilit´es li´ees au choix de l´algorithme comme dans le cas des
´equations avec conditions initiales, soit des probl`emes de stockage en m
´emoire dans le cas des ´equations avec conditions aux fronti`eres. Des
difficult´es suppl´ementaires apparaissent pour des ´equations aux d´eriv´ees
partielles non-lin´eaires (cas de l’hydrodynamique) ou` il est n´ecessaire
d’introduire des nouveaux crit`eres pour l’analyse de la stabilit´e des
algorithmes.
De tels d´eveloppements d´epassent largement le cadre de ce cours
introductif, mais cela illustre la grande complexit´e de ce sujet. En effet, Le
domaine de la r´esolution des ´equations aux d´eriv´ees partielles reste un
domaine tr`es actif de recherche en analyse num´erique, soulignant `a la
fois l’int´erˆet de nouvelles m´ethodes et la n´ecessite d’obtenir des m
´ethodes de plus en plus performantes.

1
Equations aux d´eriv´ees partielles

1
Annexe A

Coordonn´ees hypersph´eriques

Il convient de distinguer le cas pair du cas impair pour le syst`eme de

co- ordonn´ees. Les calculs analytiques et num´eriques sont plus simples
dans le cas impair, car en particulier le calcul de la transform´ee de Fourier
d’une fonction
`a sym´etrie sph´erique se ram`ene au calcul d’une transform´ee de
Fourier `a une dimension sur l’axe r´eel positif.
Consid´erons tout d’abord le cas ou` n est impair : A n dimensions, on d
´efinit les coordonn´ees hypersph´eriques de la mani`ere suivante

x1 = ρ sin(φ1) sin(φ2) sin(φ3) . . . sin(φn−3) sin(φn−2) sin(φn−1) x2

= ρ sin(φ1) sin(φ2) sin(φ3) . . . sin(φn−3) sin(φn−2) cos(φn−1) x3 =
ρ sin(φ1) sin(φ2) sin(φ3) . . . sin(φn−3) cos(φn−2)
x4 = ρ sin(φ1) sin(φ2) sin(φ3) . . . cos(φn−3)
...
xn−1 = ρ sin(φ1) cos(φ2)
xn = ρ cos(φ1) (A.1)

ou` les variables φi varient de 0 `a π pour j compris entre 1 et n − 2, la variable

φn−1 varie entre 0 et 2π, et la variable ρ est positive ou nulle.
Pour n pair, on d´efinit les coordonn´ees hypersph´eriques de la mani`ere sui-
vante

x1 = ρ sin(φ1) sin(φ2) sin(φ3) . . . sin(φn−3) sin(φn−2) sin(φn−1) x2

= ρ sin(φ1) sin(φ2) sin(φ3) . . . sin(φn−3) sin(φn−2) cos(φn−1) x3 =
ρ sin(φ1) sin(φ2) sin(φ3) . . . sin(φn−3) cos(φn−2)
x4 = ρ sin(φ1) sin(φ2) sin(φ3) . . . cos(φn−3) cos(φn−2)
...
xn−1 = ρ sin(φ2) cos(φ1)
xn = ρ cos(φ2) cos(φ1) (A.2)

On v´erifie facilement que quel que soit n

−ρ ≤ xi ≤ ρ (A.3)
1
Coordonn´ees hypersph´eriques

pour toute valeur de i. De plus, on a

n
Σ 2 2
ix =ρ
i=1
Pour ´evaluer une int´egrale `a n dimensions, il est n´ecessaire de calculer le
Jaco- bien de la transformation des coordonn´ees cart´esiennes aux coordonn
´ees hyper- sph´eriques pour n impair
∂(x1x2 . . . xn)
n−1
J =∂(ρφ1φ2 . . . φn) = ρ ×
sin(φ1) sin(φ2) . . . sin(φn−1) cos(φ1) . . . sin(φn−1) ... sin(φ1) . . . cos(φn−1)
sin(φ1) sin(φ2) . . . cos(φn−1 cos(φ1) . . . cos(φn−1) . . . − sin(φ1) . . . sin(φn−1).
. ... ... ... .
. .. .
cos(φ1) − sin(φ1) 0 (A.5)
Par r´ecurrence, on montre que le jacobien J se r´eexprime sous la forme
n−2 J
n−2
Y
− − − Y
J =ρ n 1 sin (φj ) = sin(φn−1−j )) (A.6)
j=1
n 1 j j=1

On peut v´erifier que la formule convient aussi pour n pair.

Quand l’integrant est une fonction qui ne d´epend que de ρ, on peut int
´egrer sur tous les angles et le calcul se ram`ene `a celui d’une int´egrale
simple sur l’axe r´eel positif. Par exemple on peut calculer le volume d’une
hypersph`ere de rayon a (qui correspond `a un int´egrant ´egal `a 1)
nY−2 ∫ 2π
∫a
Vn = ρn−1dρ ∫ (sin(φj )) n − 1 − j
dφj dφn−1
π
0 j=1 0
0

πn/2an
Γ(n/2
= 1) +
 pour pair
an
π n/2
(n/2)! n
= (n−1)/2 n (A.7)
π (2a) ((n − 1)/2)!
n! pour n impair
L’´el´ement de surface de l’hypersph`ere de rayon a est donn´e par
Yn−2
−
dSn = a n 1 sin − − (φi)dφi (A.8)
i=1 n 1 i
On peut en d´eduire la surface d’une hypersph`ere dans un espace `a n dimensions.
dVn
Sn = nπn/2an−1
da = Γ(n/2 + 1)!

pour n pair
n/2 an−1
nπ
(n/2)! (A.9)
= nπ
(n−1)/2
(2a)n−1(( n 1)/2)!
2
 −
1
 pour n impair
(n − 1)!

1
Annexe B

Les biblioth`eques BLAS

et Lapack

Contenu
B.1 Introduction....................................................................... 87
B.2 Terminologie...........................................................................88

B.1 Introduction
La biblioth`eque BLAS est un biblioth`eque qui r´ealise les op´erations
d’alg`ebre lin´eaire de bas niveau. Elle est disponible en Fortran et en C par
l’interm´ediaire de la biblioth`eque GSL, par exemple. Ces op´erations sont
class´ees par ordre de “difficult´e” croissante :
1. op´erations vectorielles
y = αx + y (B.1)

2. op´erations matrices-
vecteurs
y = αAx + βy (B.2)

3. op´erations matrices-matrices

C = αAB + C (B.3)

LAPACK dont le nom vient de Linear Algebra PACKage est une biblio-
th`eque ´ecrite en fortran 77 et qui comprend un ensemble de sous-
programmes pour r´esoudre les probl`emes suivants : syst`eme d’´equations lin
´eaires : probl`emes aux valeurs propres. Diverses factorisations matricielles
sont disponibles comme la d´ecomposition LU, QR Cholesky, SVD, Schur et
Schur generalis´ee. Les ma- trices consid´er´ees sont des matrices soit pleines
soit `a bande. Le cas des matrices creuses n’est pas sp´ecifiquement trait´e. Les
sous-programmes sont fournis pour des probl`emes dans lesquelles les
matrices sont soit r´eelles, soit complexes, en simple ou en double pr´ecision.

1
Les biblioth`eques BLAS et Lapack

La biblioth`eque LAPACK s´appuie sur des appels `a des sous-programmes

de plus bas niveau qui sont ceux du BLAS. Cela permet une meilleure
optimisation car des biblioth`eques BLAS sp´ecifiques existent pour chaque
type de processeur et utilisent leurs caract´eristiques. En l’absence de
biblioth`eque BLAS optimis´ee, il est possible de disposer de la biblioth`eque
BLAS de base disponible sur le site de netlib[4].
La version install´ee sur les stations de travail de l’´ecole doctorale est
la derni`ere version disponible (version 3.0) qui date de 2000. Elle est pr´evue
pour des programmes Fortran. Une version C de Lapack (clapack) est
disponible sur le site de netlib, mais n’est pas install´ee par d´efaut sur les
stations.

B.2 Terminologie
Une terminologie mn´emotechnique a ´et´e choisie pour retrouver
“facilement” les noms des sous-programmes `a appeler. Le principe est le
mˆeme pour Lapack et le BLAS.
Les sous-programmes ont un suffixe qui indique le type d’objets sur
lesquels on effectue les op´erations.
– S pour type REAL
– D pour type DOUBLE PRECISION
– C pour type COMPLEX
– Z pour type COMPLEX*16
1

Pour le BLAS de niveau 2, 3 et Lapack, un groupe de deux lettres d´esigne

le type de matrices que consid`ere le sous-programme.
– GE pour une matrice g´en´erale
– SY pour une matrice sym´etrique
– HE pour une matrice hermitienne
– TR pour une matrice triangulaire
– GB pour une matrice `a bande g´en´erale
– SB pour une matrice sym´etrique `a bande
– HB pour une matrice hermitienne `a bande
– TB pour une matrice triangulaire `a bande
Pour la recherche du nom pr´ecis d’un sous-programme, il existe un
petit index r´ef´erenc¸ant (disponible sur le site netlib ou sur les stations de
travail) de tous les sous-programmes avec la liste des arguments `a fournir. Il
existe aussi un site internet tr`es utile qui propose un moyen de rechercher
le sous-programme adapt´e `a partir de la s´election de la taˆche `a
effectuer. L’adresse de ce site est http
://[Link]/∼lapack/[Link]

1
Pour le Blas de Niveau 2, il y a un ensemble de sous-programmes pour des calculs en pr

1
´ecision ´etendue. Les pr´efixes sont ES, ED, EC, EZ

1
Annexe C

La biblioth`eque GSL

Contenu
C.1 Introduction....................................................................... 90
C.2 Sous programmes...................................................................91
C.2.1 BLAS...................................................................................91
C.2.2 fonctions simples...................................................................92
C.2.3 interfaces entre GSL et le BLAS......................................92
C.2.4 Blocs.........................................................................................92
C.2.5 S´eries de Chebyshev. ...........................................................92
C.2.6 Combinatoire..........................................................................92
C.2.7 Complexes...............................................................................92
C.2.8 Hankel......................................................................................93
C.2.9 D´eriv´ees ...............................................................................93
C.2.10 Valeurs et vecteurs propres..................................................93
C.2.11 Transform´ees de Fourier ....................................................93
C.2.12 Ajustements............................................................................93
C.2.13 Histogrammes.........................................................................94
C.2.14 IEEE....................................................................................94
C.2.15 Int´egration ............................................................................94
C.2.16 Interpolation...........................................................................94
C.2.17 Alg`ebre lin´eaire ..................................................................95
C.2.18 Matrices...................................................................................95
C.2.19 Minimisation...........................................................................95
C.2.20 Monte Carlo............................................................................95
C.2.21 Ajustements non-lin´eaires ...................................................96
C.2.22 Equations diff´erentielles ......................................................96
C.2.23 Permutation............................................................................96
C.2.24 Polynˆomes .............................................................................97
C.2.25 Puissances...............................................................................97
C.2.26 Nombres al´eatoires ..............................................................97
C.2.27 Racines unidimensionnelles..................................................97
C.2.28 Fonctions sp´eciales ..............................................................98
C.2.29 Tris.........................................................................................100
C.2.30 Lissage...................................................................................100

1
La biblioth`eque GSL

C.2.31 Statistique.............................................................................100
C.2.32 Transform´ees de Levin .....................................................101
C.2.33 Vecteurs.................................................................................101

C.1 Introduction
Parmi les nombreuses biblioth`eques math´ematiques que l’on peut
utiliser pour r´esoudre des probl`emes num´eriques, il existe une version
´evolutive bas´ee sur le concept de logiciel libre qui s’appelle la GSL (acronyme
de Gnu Scientific Library). La version actuelle est une version stable 1.2.
Cette biblioth`eque est
´ecrite en C et est install´ee sur les stations de travail de l’´ecole doctorale
(avec la version 1.1.1), ce qui permet `a tout utilisateur d’en disposer
pleinement.
Cette biblioth`eque offre des sous-programmes pour les probl`emes suivants :

Nombres complexes Racines de polynˆomes

Fonctions sp´eciales Vecteurs et Matrices
Permutations Combinatoire
Tri BLAS
Alg`ebre lin´eaire FFT
Valeurs propres G´en´erateurs de nombres al
´eatoires
Histogrammes Statistiques
Int´egration de Monte Carlo Equations diff´erentielles
Recuit simul´e Diff´erentiation num´erique
Interpolation S´eries acc´el´eratrices
Approximations de Chebyshev Recherche de racines
Transform´ees de Hankel Moindre carr´e
discr`etes
Minimisation Constantes physiques

La diversit´e des probl`emes num´eriques trait´es par cette

biblioth`eque rend tr`es difficile une pr´esentation rapide. Le manuel
proprement dit comprend ac- tuellement plus de 450 pages ! Le but de cet
appendice est donc de donner quelques cl´es pour une utilisation aussi
simple que possible de cette biblio- th`eque et par effet d’apprentissage
d’autres biblioth`eques dont les noms de sous-programmes sont g´en
´eralement moins ´eloquents.
Les fichiers d’entˆetes n´ecessaires `a la compilation d’un programme sont
pla- c´es dans le r´epertoire /usr/include/gsl.
Pour que les ex´ecutables g´en´er´es soient les plus petits possibles, les
fichiers d’entˆetes sont tr`es nombreux (206 dans la version 1.1.1), et leur
nom corres- pondent `a une classe de sous-programes que l’on souhaite utilis
´es selon la r`egle suivante :
1
– Le nom commence toujours par le pr´efixe gsl . (Moyen simple
d’identifi- cation par rapport `a d’autres biblioth`eques qui seraient `a
utiliser dans le mˆeme programme).
– La seconde partie correspond `a la classe de probl`eme. Par exemple sf
cor- respond `a un fonction sp´eciale, matrix `a un probl`eme sur les
matrices,...

1
 C.2 Sous programmes

– La troisi`eme correspond au probl`eme sp´ecifique. Par exemple, gsl sf

expint.h est le fichier pour les fonctions exponentielles int´egrales, gsl sf
bessel.h aux fonctions de bessel, ...
– Si l’on veut disposer d’une classe de fonctions plus directement au lieu
de charger un `a un les fichiers d’entˆetes, on peut utiliser le fichier
qui regroupe. Pour l’exemple pr´ec´edent, on peut utiliser gsl sf.h qui
contient les d´efinitions de toutes les fonctions sp´eciales.

C.2 Sous programmes

La Liste des fonctions disponibles est la suivante

C.2.1 BLAS
Les fonctions du BLAS ´ecrites. Elles r´ealisent les trois niveaux qui
sont les suivants :
1. op´erations vectorielles
y = αx + y (C.1)
2. op´erations matrices-
vecteurs
y = αAx + βy (C.2)

3. op´erations matrices-
matrices C = αAB + C (C.3)

cblas caxpy, cblas ccopy, cblas cdotc sub, cblas cdotu sub, cblas cgbmv, cblas cgemm,
cblas cgemv, cblas cgerc, cblas cgeru, cblas chbmv, cblas chemm, cblas chemv, cblas cher,
cblas cher2, cblas cher2k, cblas cherk, cblas chpmv, cblas chpr, cblas chpr2, cblas cscal,
cblas csscal, cblas cswap, cblas csymm, cblas csyr2k, cblas csyrk, cblas ctbmv, cblas ctbsv,
cblas ctpmv, cblas ctpsv, cblas ctrmm, cblas ctrmv, cblas ctrsm, cblas ctrsv, cblas dasum,
cblas daxpy, cblas dcopy, cblas ddot, cblas dgbmv, cblas dgemm, cblas dgemv, cblas dger,
cblas dnrm2, cblas drot, cblas drotg, cblas drotm, cblas drotmg, cblas dsbmv, cblas dscal,
cblas dsdot, cblas dspmv, cblas dspr, cblas dspr2, cblas dswap, cblas dsymm, cblas dsymv,
cblas dsyr, cblas dsyr2, cblas dsyr2k, cblas dsyrk, cblas dtbmv, cblas dtbsv, cblas dtpmv,
cblas dtpsv, cblas dtrmm, cblas dtrmv, cblas dtrsm, cblas dtrsv, cblas dzasum, cblas dznrm2,
cblas icamax, cblas idamax, cblas isamax, cblas izamax, cblas sasum, cblas saxpy, cblas scasum,
cblas scnrm2, cblas scopy, cblas sdot, cblas sdsdot, cblas sgbmv, cblas sgemm, cblas sgemv,
cblas sger, cblas snrm2, cblas srot, cblas srotg, cblas srotm, cblas srotmg, cblas ssbmv,
cblas sscal, cblas sspmv, cblas sspr, cblas sspr2, cblas sswap, cblas ssymm, cblas ssymv,
cblas ssyr, cblas ssyr2, cblas ssyr2k, cblas ssyrk, cblas stbmv, cblas stbsv, cblas stpmv,
cblas stpsv, cblas strmm, cblas strmv, cblas strsm, cblas strsv, cblas xerbla, cblas zaxpy,
cblas zcopy, cblas zdotc sub, cblas zdotu sub, cblas zdscal, cblas zgbmv, cblas zgemm,
cblas zgemv, cblas zgerc, cblas zgeru, cblas zhbmv, cblas zhemm, cblas zhemv, cblas zher,
cblas zher2, cblas zher2k, cblas zherk, cblas zhpmv, cblas zhpr, cblas zhpr2, cblas zscal,
cblas zswap, cblas zsymm, cblas zsyr2k, cblas zsyrk, cblas ztbmv, cblas ztbsv, cblas ztpmv,
cblas ztpsv, cblas ztrmm, cblas ztrmv, cblas ztrsm, cblas ztrsv

1
La biblioth`eque GSL

C.2.2 fonctions simples

Les fonctions simples avec une grande pr´ecision
gsl acosh, gsl asinh, gsl atanh, gsl expm1, gsl log1p

C.2.3 interfaces entre GSL et le BLAS

gsl blas caxpy, gsl blas ccopy, gsl blas cdotc, gsl blas cdotu, gsl blas cgemm, gsl blas cgemv,
gsl blas cgerc, gsl blas cgeru, gsl blas chemm, gsl blas chemv, gsl blas cher, gsl blas cher2,
gsl blas cher2k, gsl blas cherk, gsl blas cscal, gsl blas csscal, gsl blas cswap, gsl blas csymm,
gsl blas csyr2k, gsl blas csyrk, gsl blas ctrmm, gsl blas ctrmv, gsl blas ctrsm, gsl blas ctrsv,
gsl blas dasum, gsl blas daxpy, gsl blas dcopy, gsl blas ddot, gsl blas dgemm, gsl blas dgemv,
gsl blas dger, gsl blas dnrm2, gsl blas drot, gsl blas drotg, gsl blas drotm, gsl blas drotmg,
gsl blas dscal, gsl blas dsdot, gsl blas dswap, gsl blas dsymm, gsl blas dsymv, gsl blas dsyr,
gsl blas dsyr2, gsl blas dsyr2k, gsl blas dsyrk, gsl blas dtrmm, gsl blas dtrmv, gsl blas dtrsm,
gsl blas dtrsv, gsl blas dzasum, gsl blas dznrm2, gsl blas icamax, gsl blas idamax, gsl blas isamax,
gsl blas izamax, gsl blas sasum, gsl blas saxpy, gsl blas scasum, gsl blas scnrm2, gsl blas scopy,
gsl blas sdot, gsl blas sdsdot, gsl blas sgemm, gsl blas sgemv, gsl blas sger, gsl blas snrm2,
gsl blas srot, gsl blas srotg, gsl blas srotm, gsl blas srotmg, gsl blas sscal, gsl blas sswap,
gsl blas ssymm, gsl blas ssymv, gsl blas ssyr, gsl blas ssyr2, gsl blas ssyr2k, gsl blas ssyrk,
gsl blas strmm, gsl blas strmv, gsl blas strsm, gsl blas strsv, gsl blas zaxpy, gsl blas zcopy,
gsl blas zdotc, gsl blas zdotu, gsl blas zdscal, gsl blas zgemm, gsl blas zgemv, gsl blas zgerc,
gsl blas zgeru, gsl blas zhemm, gsl blas zhemv, gsl blas zher, gsl blas zher2, gsl blas zher2k,
gsl blas zherk, gsl blas zscal, gsl blas zswap, gsl blas zsymm, gsl blas zsyr2k, gsl blas zsyrk,
gsl blas ztrmm, gsl blas ztrmv, gsl blas ztrsm, gsl blas ztrsv

C.2.4 Blocs
Les fonctions pour les blocs de m´emoire , gsl block alloc, gsl block
calloc, gsl block fprintf, gsl block fread, gsl block free, gsl block fscanf, gsl block
fwrite

C.2.5 S´eries de Chebyshev.

gsl cheb alloc, gsl cheb calc deriv, gsl cheb calc integ, gsl cheb eval, gsl cheb eval err,
gsl cheb eval n, gsl cheb eval n err, gsl cheb free, gsl cheb init

C.2.6 Combinatoire
gsl combination alloc, gsl combination calloc, gsl combination data, gsl combination fprintf,
gsl combination fread, gsl combination free , gsl combination fscanf, gsl combination fwrite,
gsl combination get, gsl combination init first, gsl combination init last, gsl combination k,
gsl combination n, gsl combination next, gsl combination prev , gsl combination valid

C.2.7 Complexes
Les fonctions d’op´erations sur les complexes gsl complex abs, gsl complex abs2,
gsl complex add, gsl complex add imag, gsl complex add real, gsl complex arccos,

1
gsl complex arccos real, gsl complex arccosh, gsl complex arccosh real, gsl complex arccot,

1
 C.2 Sous programmes

gsl complex arccoth, gsl complex arccsc, gsl complex arccsc real, gsl complex arccsch,
gsl complex arcsec, gsl complex arcsec real, gsl complex arcsech, gsl complex arcsin ,
gsl complex arcsin real, gsl complex arcsinh, gsl complex arctan, gsl complex arctanh,
gsl complex arctanh real, gsl complex arg, gsl complex conjugate, gsl complex cos, gsl complex cosh,
gsl complex cot, gsl complex coth, gsl complex csc, gsl complex csch, gsl complex div,
gsl complex div imag, gsl complex div real, gsl complex exp, gsl complex inverse, gsl complex linalg LU invert,
gsl complex log, gsl complex log10, gsl complex log b, gsl complex logabs, gsl complex mul,
gsl complex mul imag, gsl complex mul real, gsl complex negative, gsl complex polar,
gsl complex pow, gsl complex pow real, gsl complex rect, gsl complex sec, gsl complex sech,
gsl complex sin, gsl complex sinh, gsl complex sqrt, gsl complex sqrt real, gsl complex sub,
gsl complex sub imag, gsl complex sub real, gsl complex tan, gsl complex tanh

C.2.8 Hankel
Les fonctions pour les transform´ees de Hankel
gsl dht alloc, gsl dht apply, gsl dht free, gsl dht init, gsl dht k sample, gsl dht new,
gsl dht x sample

C.2.9 D´eriv´ees
gsl diff backward, gsl diff central, gsl diff forward

C.2.10 Valeurs et vecteurs propres

gsl eigen herm, gsl eigen herm alloc, gsl eigen herm free, gsl eigen hermv,
gsl eigen hermv alloc, gsl eigen hermv free, gsl eigen hermv sort, gsl eigen symm,
gsl eigen symm alloc, gsl eigen symm free, gsl eigen symmv, gsl eigen symmv alloc,
gsl eigen symmv free, gsl eigen symmv sort

C.2.11 Transform´ees de Fourier

gsl fft complex backward, gsl fft complex forward, gsl fft complex inverse,
gsl fft complex radix2 backward, gsl fft complex radix2 dif backward, gsl fft complex radix2 dif forward,
gsl fft complex radix2 dif inverse, gsl fft complex radix2 dif transform, gsl fft complex radix2 forward,
gsl fft complex radix2 inverse, gsl fft complex radix2 transform, gsl fft complex transform,
gsl fft complex wavetable alloc, gsl fft complex wavetable free, gsl fft complex workspace alloc,
gsl fft complex workspace free, gsl fft halfcomplex radix2 backward, gsl fft halfcomplex radix2 inverse,
gsl fft halfcomplex transform, gsl fft halfcomplex unpack, gsl fft halfcomplex wavetable alloc,
gsl fft halfcomplex wavetable free, gsl fft real radix2 transform, gsl fft real transform,
gsl fft real unpack, gsl fft real wavetable alloc, gsl fft real wavetable free, gsl fft real workspace alloc,
gsl fft real workspace free
gsl finite

C.2.12 Ajustements
gsl fit linear, gsl fit linear est, gsl fit mul, gsl fit mul est, gsl fit wlinear, gsl fit wmul

1
La biblioth`eque GSL

C.2.13 Histogrammes
gsl histogram2d accumulate, gsl histogram2d add, gsl histogram2d alloc, gsl histogram2d clone,
gsl histogram2d cov, gsl histogram2d div, gsl histogram2d equal bins p, gsl histogram2d find,
gsl histogram2d fprintf, gsl histogram2d fread , gsl histogram2d free, gsl histogram2d fscanf,
gsl histogram2d fwrite , gsl histogram2d get, gsl histogram2d get xrange, gsl histogram2d get yrange,
gsl histogram2d increment, gsl histogram2d max bin, gsl histogram2d max val, gsl histogram2d memcpy,
gsl histogram2d min bin, gsl histogram2d min val, gsl histogram2d mul, gsl histogram2d nx,
gsl histogram2d ny, gsl histogram2d pdf alloc, gsl histogram2d pdf free, gsl histogram2d pdf init,
gsl histogram2d pdf sample, gsl histogram2d reset, gsl histogram2d scale, gsl histogram2d set ranges,
gsl histogram2d set ranges uniform, gsl histogram2d shift, gsl histogram2d sub, gsl histogram2d sum,
gsl histogram2d xmax, gsl histogram2d xmean, gsl histogram2d xmin, gsl histogram2d xsigma,
gsl histogram2d ymax, gsl histogram2d ymean, gsl histogram2d ymin, gsl histogram2d ysigma,
gsl histogram accumulate, gsl histogram add, gsl histogram alloc, gsl histogram bins,
gsl histogram clone, gsl histogram div, gsl histogram equal bins p, gsl histogram find,
gsl histogram fprintf, gsl histogram fread, gsl histogram free, gsl histogram fscanf, gsl histogram fwrite,
gsl histogram get, gsl histogram get range, gsl histogram increment, gsl histogram max,
gsl histogram max bin, gsl histogram max val, gsl histogram mean, gsl histogram memcpy,
gsl histogram min, gsl histogram min bin, gsl histogram min val, gsl histogram mul,
gsl histogram pdf alloc, gsl histogram pdf free , gsl histogram pdf init, gsl histogram pdf sample,
gsl histogram reset, gsl histogram scale, gsl histogram set ranges, gsl histogram set ranges uniform,
gsl histogram shift, gsl histogram sigma, gsl histogram sub, gsl histogram sum
gsl hypot

C.2.14 IEEE
gsl ieee env setup, gsl ieee fprintf double, gsl ieee fprintf float, gsl ieee printf double,
gsl ieee printf float

C.2.15 Int´egration
gsl integration qag, gsl integration qagi, gsl integration qagil, gsl integration qagiu,
gsl integration qagp, gsl integration qags, gsl integration qawc, gsl integration qawf,
gsl integration qawo, gsl integration qawo table alloc, gsl integration qawo table free,
gsl integration qawo table set, gsl integration qawo table set length, gsl integration qaws,
gsl integration qaws table alloc, gsl integration qaws table free, gsl integration qaws table set,
gsl integration qng, gsl integration workspace alloc, gsl integration workspace free

C.2.16 Interpolation
gsl interp accel alloc, gsl interp accel find , gsl interp accel free, gsl interp
akima, gsl interp akima periodic, gsl interp alloc, gsl interp bsearch, gsl interp
cspline,
gsl interp cspline periodic, gsl interp eval, gsl interp eval deriv, gsl interp eval deriv2,
gsl interp eval deriv2 e, gsl interp eval deriv e, gsl interp eval e, gsl interp eval integ,
gsl interp eval integ e, gsl interp free, gsl interp init, gsl interp linear, gsl interp min size,
gsl interp name, gsl interp polynomial

1
gsl isinf, gsl isnan

1
 C.2 Sous programmes

C.2.17 Alg`ebre lin´eaire

gsl linalg bidiag decomp, gsl linalg bidiag unpack, gsl linalg bidiag unpack2, gsl linalg bidiag unpack B,
gsl linalg cholesky decomp, gsl linalg cholesky solve, gsl linalg cholesky svx, gsl linalg complex LU decomp,
gsl linalg complex LU det, gsl linalg complex LU lndet, gsl linalg complex LU refine,
gsl linalg complex LU sgndet, gsl linalg complex LU solve, gsl linalg complex LU svx,
gsl linalg hermtd decomp, gsl linalg hermtd unpack, gsl linalg hermtd unpack T, gsl linalg HH solve,
gsl linalg HH svx, gsl linalg LU decomp, gsl linalg LU det, gsl linalg LU invert, gsl linalg LU lndet,
gsl linalg LU refine, gsl linalg LU sgndet, gsl linalg LU solve, gsl linalg LU svx, gsl linalg QR decomp,
gsl linalg QR lssolve, gsl linalg QR QRsolve, gsl linalg QR QTvec, gsl linalg QR Qvec,
gsl linalg QR Rsolve, gsl linalg QR Rsvx, gsl linalg QR solve, gsl linalg QR svx, gsl linalg QR unpack,
gsl linalg QR update, gsl linalg QRPT decomp , gsl linalg QRPT decomp2, gsl linalg QRPT QRsolve,
gsl linalg QRPT Rsolve, gsl linalg QRPT Rsvx, gsl linalg QRPT solve, gsl linalg QRPT svx,
gsl linalg QRPT update, gsl linalg R solve, gsl linalg R svx, gsl linalg solve symm cyc tridiag,
gsl linalg solve symm tridiag, gsl linalg SV decomp, gsl linalg SV decomp jacobi, gsl linalg SV decomp mod,
gsl linalg SV solve, gsl linalg symmtd decomp, gsl linalg symmtd unpack, gsl linalg symmtd unpack T

C.2.18 Matrices
gsl matrix add, gsl matrix add constant, gsl matrix alloc, gsl matrix calloc, gsl matrix column,
gsl matrix const column, gsl matrix const diagonal , gsl matrix const row, gsl matrix const subdiagonal,
gsl matrix const submatrix, gsl matrix const superdiagonal, gsl matrix const view array,
gsl matrix const view array with tda, gsl matrix const view vector, gsl matrix const view vector with tda,
gsl matrix diagonal, gsl matrix div elements, gsl matrix fprintf, gsl matrix fread, gsl matrix free,
gsl matrix fscanf, gsl matrix fwrite, gsl matrix get, gsl matrix get col, gsl matrix get row,
gsl matrix isnull, gsl matrix max, gsl matrix max index, gsl matrix memcpy, gsl matrix min,
gsl matrix min index, gsl matrix minmax, gsl matrix minmax index, gsl matrix mul elements,
gsl matrix ptr, gsl matrix ptr, gsl matrix row, gsl matrix scale, gsl matrix set, gsl matrix set all,
gsl matrix set col, gsl matrix set identity, gsl matrix set row, gsl matrix set zero, gsl matrix sub,
gsl matrix subdiagonal, gsl matrix submatrix, gsl matrix superdiagonal, gsl matrix swap,
gsl matrix swap columns, gsl matrix swap rowcol, gsl matrix swap rows, gsl matrix transpose,
gsl matrix transpose memcpy , gsl matrix view array, gsl matrix view array with tda,
gsl matrix view vector, gsl matrix view vector with tda

C.2.19 Minimisation
gsl min fminimizer alloc, gsl min fminimizer brent, gsl min fminimizer f lower, gsl min fminimizer f minimum,
gsl min fminimizer f upper, gsl min fminimizer free, gsl min fminimizer goldensection,
gsl min fminimizer iterate, gsl min fminimizer name, gsl min fminimizer set , gsl min fminimizer set with values,
gsl min fminimizer x lower, gsl min fminimizer x minimum, gsl min fminimizer x upper,
gsl min test interval

C.2.20 Monte Carlo

gsl monte miser alloc, gsl monte miser free, gsl monte miser init , gsl monte miser integrate,
gsl monte plain alloc, gsl monte plain free, gsl monte plain init, gsl monte plain integrate,

1
La biblioth`eque GSL

gsl monte vegas alloc, gsl monte vegas free, gsl monte vegas init , gsl monte vegas integrate

C.2.21 Ajustements non-lin´eaires

gsl multifit covar, gsl multifit fdfsolver alloc, gsl multifit fdfsolver free, gsl multifit fdfsolver iterate,
gsl multifit fdfsolver lmder, gsl multifit fdfsolver lmsder, gsl multifit fdfsolver name,
gsl multifit fdfsolver position, gsl multifit fdfsolver set, gsl multifit fsolver alloc, gsl multifit fsolver free,
gsl multifit fsolver iterate, gsl multifit fsolver name, gsl multifit fsolver position, gsl multifit fsolver set,
gsl multifit gradient, gsl multifit linear, gsl multifit linear alloc, gsl multifit linear free,
gsl multifit test delta, gsl multifit test gradient, gsl multifit wlinear, gsl multimin fdfminimizer alloc,
gsl multimin fdfminimizer conjugate fr, gsl multimin fdfminimizer conjugate pr,
gsl multimin fdfminimizer free, gsl multimin fdfminimizer gradient, gsl multimin fdfminimizer iterate,
gsl multimin fdfminimizer minimum, gsl multimin fdfminimizer name, gsl multimin fdfminimizer restart,
gsl multimin fdfminimizer set, gsl multimin fdfminimizer steepest descent,
gsl multimin fdfminimizer vector bfgs, gsl multimin fdfminimizer x, gsl multimin test gradient,
gsl multiroot fdfsolver alloc, gsl multiroot fdfsolver dx, gsl multiroot fdfsolver f,
gsl multiroot fdfsolver free, gsl multiroot fdfsolver gnewton, gsl multiroot fdfsolver hybridj,
gsl multiroot fdfsolver hybridsj, gsl multiroot fdfsolver iterate, gsl multiroot fdfsolver name,
gsl multiroot fdfsolver newton, gsl multiroot fdfsolver root, gsl multiroot fdfsolver set,
gsl multiroot fsolver alloc, gsl multiroot fsolver broyden, gsl multiroot fsolver dnewton,
gsl multiroot fsolver dx , gsl multiroot fsolver f, gsl multiroot fsolver free, gsl multiroot fsolver hybrid,
gsl multiroot fsolver hybrids, gsl multiroot fsolver iterate, gsl multiroot fsolver name,
gsl multiroot fsolver root, gsl multiroot fsolver set, gsl multiroot test delta, gsl multiroot test residual
gsl ntuple bookdata, gsl ntuple close, gsl ntuple create, gsl ntuple open, gsl ntuple project,
gsl ntuple read, gsl ntuple write.

C.2.22 Equations diff´erentielles

gsl odeiv control alloc, gsl odeiv control free, gsl odeiv control hadjust, gsl odeiv control init,
gsl odeiv control name, gsl odeiv control standard new, gsl odeiv control y new, gsl odeiv control yp new,
gsl odeiv evolve alloc, gsl odeiv evolve apply, gsl odeiv evolve free, gsl odeiv evolve reset,
gsl odeiv step alloc, gsl odeiv step apply, gsl odeiv step bsimp, gsl odeiv step free, gsl odeiv step gear1,
gsl odeiv step gear2, gsl odeiv step name, gsl odeiv step order, gsl odeiv step reset,
gsl odeiv step rk2, gsl odeiv step rk2imp, gsl odeiv step rk4, gsl odeiv step rk4imp,
gsl odeiv step rk8pd, gsl odeiv step rkck, gsl odeiv step rkf45

C.2.23 Permutation
gsl permutation alloc, gsl permutation calloc, gsl permutation canonical cycles, gsl permutation canonical to
gsl permutation data, gsl permutation fprintf, gsl permutation fread, gsl permutation free,
gsl permutation fscanf, gsl permutation fwrite, gsl permutation get, gsl permutation init,
gsl permutation inverse, gsl permutation inversions, gsl permutation linear cycles, gsl permutation linear to cano
gsl permutation memcpy, gsl permutation mul, gsl permutation next, gsl permutation prev,
gsl permutation reverse, gsl permutation size, gsl permutation swap, gsl permutation valid,
gsl permute, gsl permute inverse, gsl permute vector, gsl permute vector inverse

1
 C.2 Sous programmes

C.2.24 Polynˆomes
gsl poly complex solve, gsl poly complex solve cubic, gsl poly complex solve quadratic,
gsl poly complex workspace alloc, gsl poly complex workspace free, gsl poly dd eval,
gsl poly dd init, gsl poly dd taylor, gsl poly eval, gsl poly solve cubic, gsl poly solve quadratic

C.2.25 Puissances
gsl pow 2, gsl pow 3, gsl pow 4, gsl pow 5, gsl pow 6, gsl pow 7, gsl pow 8, gsl pow 9,
gsl pow int

C.2.26 Nombres al´eatoires

gsl qrng alloc, gsl qrng clone, gsl qrng free, gsl qrng get , gsl qrng init, gsl qrng memcpy,
gsl qrng name, gsl qrng niederreiter 2, gsl qrng size, gsl qrng sobol, gsl qrng state, gsl ran bernoulli,
gsl ran bernoulli pdf, gsl ran beta, gsl ran beta pdf, gsl ran binomial, gsl ran binomial pdf,
gsl ran bivariate gaussian, gsl ran bivariate gaussian pdf, gsl ran cauchy, gsl ran cauchy pdf,
gsl ran chisq, gsl ran chisq pdf, gsl ran choose, gsl ran dir 2d, gsl ran dir 2d trig method,
gsl ran dir 3d, gsl ran dir nd, gsl ran discrete, gsl ran discrete free, gsl ran discrete pdf,
gsl ran discrete preproc, gsl ran exponential, gsl ran exponential pdf, gsl ran exppow,
gsl ran exppow pdf, gsl ran fdist, gsl ran fdist pdf, gsl ran flat, gsl ran flat pdf, gsl ran gamma,
gsl ran gamma pdf, gsl ran gaussian, gsl ran gaussian pdf, gsl ran gaussian ratio method,
gsl ran gaussian tail, gsl ran gaussian tail pdf, gsl ran geometric , gsl ran geometric pdf,
gsl ran gumbel1, gsl ran gumbel1 pdf, gsl ran gumbel2, gsl ran gumbel2 pdf, gsl ran hypergeometric,
gsl ran hypergeometric pdf, gsl ran landau, gsl ran landau pdf, gsl ran laplace, gsl ran laplace pdf,
gsl ran levy, gsl ran levy skew, gsl ran logarithmic, gsl ran logarithmic pdf, gsl ran logistic,
gsl ran logistic pdf, gsl ran lognormal, gsl ran lognormal pdf, gsl ran negative binomial,
gsl ran negative binomial pdf, gsl ran pareto, gsl ran pareto pdf, gsl ran pascal, gsl ran pascal pdf,
gsl ran poisson, gsl ran poisson pdf, gsl ran rayleigh, gsl ran rayleigh pdf, gsl ran rayleigh tail,
gsl ran rayleigh tail pdf, gsl ran sample, gsl ran shuffle, gsl ran tdist, gsl ran tdist pdf,
gsl ran ugaussian, gsl ran ugaussian pdf, gsl ran ugaussian ratio method, gsl ran ugaussian tail,
gsl ran ugaussian tail pdf, gsl ran weibull, gsl ran weibull pdf, GSL REAL, gsl rng alloc,
gsl rng borosh13, gsl rng clone, gsl rng cmrg, gsl rng coveyou, gsl rng env setup, gsl rng fishman18,
gsl rng fishman20, gsl rng fishman2x, gsl rng free, gsl rng get, gsl rng gfsr4, gsl rng knuthran,
gsl rng knuthran2, gsl rng lecuyer21, gsl rng max, gsl rng memcpy, gsl rng min, gsl rng minstd,
gsl rng mrg, gsl rng mt19937, gsl rng name, gsl rng print state, gsl rng r250, gsl rng ran0,
gsl rng ran1, gsl rng ran2, gsl rng ran3, gsl rng rand, gsl rng rand48, gsl rng random bsd,
gsl rng random glibc2, gsl rng random libc5, gsl rng randu, gsl rng ranf, gsl rng ranlux,
gsl rng ranlux389, gsl rng ranlxd1, gsl rng ranlxd2, gsl rng ranlxs0, gsl rng ranlxs1,
gsl rng ranlxs2, gsl rng ranmar, gsl rng set, gsl rng size, gsl rng slatec, gsl rng state,
gsl rng taus, gsl rng taus2, gsl rng transputer, gsl rng tt800, gsl rng types setup, gsl rng uni,
gsl rng uni32, gsl rng uniform, gsl rng uniform int, gsl rng uniform pos, gsl rng vax,
gsl rng waterman14, gsl rng zuf

C.2.27 Racines unidimensionnelles

gsl root fdfsolver alloc, gsl root fdfsolver free, gsl root fdfsolver iterate, gsl root fdfsolver name,

1
La biblioth`eque GSL

gsl root fdfsolver newton, gsl root fdfsolver root, gsl root fdfsolver secant, gsl root fdfsolver set,
gsl root fdfsolver steffenson, gsl root fsolver alloc, gsl root fsolver bisection, gsl root fsolver brent,
gsl root fsolver falsepos, gsl root fsolver free, gsl root fsolver iterate, gsl root fsolver name,
gsl root fsolver root, gsl root fsolver set, gsl root fsolver x lower, gsl root fsolver x upper,
gsl root test delta, gsl root test interval , gsl root test residual
gsl set error handler, gsl set error handler off

C.2.28 Fonctions sp´eciales

Airy
gsl sf airy Ai, gsl sf airy Ai deriv, gsl sf airy Ai deriv e, gsl sf airy Ai deriv scaled,
gsl sf airy Ai deriv scaled e, gsl sf airy Ai e, gsl sf airy Ai scaled, gsl sf airy Ai scaled e,
gsl sf airy Bi, gsl sf airy Bi deriv, gsl sf airy Bi deriv e, gsl sf airy Bi deriv scaled,
gsl sf airy Bi deriv scaled e, gsl sf airy Bi e, gsl sf airy Bi scaled, gsl sf airy Bi scaled e,
gsl sf airy zero Ai, gsl sf airy zero Ai deriv, gsl sf airy zero Ai deriv e, gsl sf airy zero Ai e,
gsl sf airy zero Bi, gsl sf airy zero Bi deriv, gsl sf airy zero Bi deriv e, gsl sf airy zero Bi e

Divers
gsl sf angle restrict pos, gsl sf angle restrict pos e, gsl sf angle restrict symm,
gsl sf angle restrict symm e, gsl sf atanint, gsl sf atanint e

Bessel
gsl sf bessel I0, gsl sf bessel I0 e, gsl sf bessel I0 scaled, gsl sf bessel i0 scaled, gsl sf bessel I0 scaled e,
gsl sf bessel i0 scaled e, gsl sf bessel I1, gsl sf bessel I1 e, gsl sf bessel I1 scaled, gsl sf bessel i1 scaled,
gsl sf bessel I1 scaled e, gsl sf bessel i1 scaled e, gsl sf bessel i2 scaled, gsl sf bessel i2 scaled
e, gsl sf bessel il scaled, gsl sf bessel il scaled array, gsl sf bessel il scaled e, gsl sf bessel In,
gsl sf bessel In array, gsl sf bessel In e, gsl sf bessel In scaled, gsl sf bessel In scaled array,
gsl sf bessel In scaled e, gsl sf bessel Inu, gsl sf bessel Inu e, gsl sf bessel Inu scaled,
gsl sf bessel Inu scaled e, gsl sf bessel J0, gsl sf bessel j0, gsl sf bessel j0 e, gsl sf bessel J0 e,
gsl sf bessel J1, gsl sf bessel j1, gsl sf bessel J1 e, gsl sf bessel j1 e, gsl sf bessel j2 ,
gsl sf bessel j2 e, gsl sf bessel jl, gsl sf bessel jl array, gsl sf bessel jl e, gsl sf bessel jl steed array,
gsl sf bessel Jn, gsl sf bessel Jn array, gsl sf bessel Jn e, gsl sf bessel Jnu, gsl sf bessel Jnu e,
gsl sf bessel K0, gsl sf bessel K0 e, gsl sf bessel K0 scaled, gsl sf bessel k0 scaled,
gsl sf bessel K0 scaled e, gsl sf bessel k0 scaled e, gsl sf bessel K1, gsl sf bessel K1 e,
gsl sf bessel K1 scaled, gsl sf bessel k1 scaled, gsl sf bessel K1 scaled e, gsl sf bessel k1 scaled e,
gsl sf bessel k2 scaled, gsl sf bessel k2 scaled e, gsl sf bessel kl scaled, gsl sf bessel kl scaled array,
gsl sf bessel kl scaled e, gsl sf bessel Kn, gsl sf bessel Kn array, gsl sf bessel Kn e, gsl sf bessel Kn scaled,
gsl sf bessel Kn scaled array, gsl sf bessel Kn scaled e, gsl sf bessel Knu, gsl sf bessel Knu e,
gsl sf bessel Knu scaled, gsl sf bessel Knu scaled e, gsl sf bessel lnKnu, gsl sf bessel lnKnu e,
gsl sf bessel sequence Jnu e, gsl sf bessel y0, gsl sf bessel Y0, gsl sf bessel y0 e, gsl sf bessel Y0 e,
gsl sf bessel y1, gsl sf bessel Y1, gsl sf bessel Y1 e, gsl sf bessel y1 e, gsl sf bessel y2,
gsl sf bessel y2 e, gsl sf bessel yl, gsl sf bessel yl array, gsl sf bessel yl e, gsl sf bessel Yn,
gsl sf bessel Yn array, gsl sf bessel Yn e, gsl sf bessel Ynu, gsl sf bessel Ynu e, gsl sf bessel zero J0,
gsl sf bessel zero J0 e, gsl sf bessel zero J1, gsl sf bessel zero J1 e, gsl sf bessel zero Jnu,
gsl sf bessel zero Jnu e

1
 C.2 Sous programmes

Beta, Coniques,Elliptiques, ...

gsl sf beta, gsl sf beta e, gsl sf beta inc, gsl sf beta inc e, gsl sf Chi, gsl sf Chi e,
gsl sf choose, gsl sf choose e, gsl sf Ci, gsl sf Ci e, gsl sf clausen, gsl sf clausen e, gsl sf complex cos e,
gsl sf complex dilog e, gsl sf complex log e, gsl sf complex logsin e, gsl sf complex sin e,
gsl sf conicalP 0, gsl sf conicalP 0 e, gsl sf conicalP 1 , gsl sf conicalP 1 e, gsl sf conicalP cyl reg,
gsl sf conicalP cyl reg e, gsl sf conicalP half, gsl sf conicalP half e, gsl sf conicalP mhalf,
gsl sf conicalP mhalf e, gsl sf conicalP sph reg, gsl sf conicalP sph reg e, gsl sf cos,
gsl sf cos e, gsl sf cos err, gsl sf cos err e, gsl sf coulomb CL array , gsl sf coulomb CL e,
gsl sf coulomb wave F array, gsl sf coulomb wave FG array, gsl sf coulomb wave FG e,
gsl sf coulomb wave FGp array , gsl sf coulomb wave sphF array, gsl sf coupling 3j,
gsl sf coupling 3j e, gsl sf coupling 6j, gsl sf coupling 6j e, gsl sf coupling 9j, gsl sf coupling 9j e,
gsl sf dawson, gsl sf dawson e, gsl sf debye 1, gsl sf debye 1 e, gsl sf debye 2, gsl sf debye 2 e,
gsl sf debye 3, gsl sf debye 3 e, gsl sf debye 4, gsl sf debye 4 e, gsl sf dilog, gsl sf dilog e,
gsl sf doublefact, gsl sf doublefact e, gsl sf ellint D, gsl sf ellint D e, gsl sf ellint E,
gsl sf ellint E e, gsl sf ellint Ecomp, gsl sf ellint Ecomp e, gsl sf ellint F, gsl sf ellint F
e, gsl sf ellint Kcomp, gsl sf ellint Kcomp e, gsl sf ellint P, gsl sf ellint P e, gsl sf ellint
RC, gsl sf ellint RC e, gsl sf ellint RD, gsl sf ellint RD e, gsl sf ellint RF, gsl sf ellint RF
e,
gsl sf ellint RJ, gsl sf ellint RJ e, gsl sf elljac e, gsl sf erf, gsl sf erf e, gsl sf erf Q, gsl sf erf Q
e, gsl sf erf Z, gsl sf erf Z e, gsl sf erfc, gsl sf erfc e, gsl sf eta, gsl sf eta e, gsl sf eta int,
gsl sf eta int e, gsl sf exp, gsl sf exp e, gsl sf exp e10 e , gsl sf exp err e, gsl sf exp err e10 e,
gsl sf exp mult, gsl sf exp mult e, gsl sf exp mult e10 e, gsl sf exp mult err e, gsl sf exp mult err e10 e,
gsl sf expint 3, gsl sf expint 3 e, gsl sf expint E1, gsl sf expint E1 e, gsl sf expint E2,
gsl sf expint E2 e, gsl sf expint Ei, gsl sf expint Ei e, gsl sf expm1, gsl sf expm1 e,
gsl sf exprel, gsl sf exprel 2, gsl sf exprel 2 e, gsl sf exprel e, gsl sf exprel n, gsl sf exprel n e,
gsl sf fact, gsl sf fact e, gsl sf fermi dirac 0, gsl sf fermi dirac 0 e, gsl sf fermi dirac 1,
gsl sf fermi dirac 1 e, gsl sf fermi dirac 2, gsl sf fermi dirac 2 e, gsl sf fermi dirac 3half,
gsl sf fermi dirac 3half e, gsl sf fermi dirac half, gsl sf fermi dirac half e, gsl sf fermi dirac inc 0,
gsl sf fermi dirac inc 0 e, gsl sf fermi dirac int, gsl sf fermi dirac int e, gsl sf fermi dirac m1,
gsl sf fermi dirac m1 e, gsl sf fermi dirac mhalf, gsl sf fermi dirac mhalf e, gsl sf gamma,
gsl sf gamma e, gsl sf gamma inc P, gsl sf gamma inc P e, gsl sf gamma inc Q, gsl sf gamma inc Q e,
gsl sf gammainv , gsl sf gammainv e, gsl sf gammastar, gsl sf gammastar e, gsl sf gegenpoly 1,
gsl sf gegenpoly 1 e, gsl sf gegenpoly 2, gsl sf gegenpoly 2 e, gsl sf gegenpoly 3, gsl sf gegenpoly 3 e,
gsl sf gegenpoly array, gsl sf gegenpoly n, gsl sf gegenpoly n e, gsl sf hydrogenicR, gsl sf hydrogenicR 1,
gsl sf hydrogenicR 1 e, gsl sf hydrogenicR e, gsl sf hyperg 0F1, gsl sf hyperg 0F1 e,
gsl sf hyperg 1F1, gsl sf hyperg 1F1 e, gsl sf hyperg 1F1 int, gsl sf hyperg 1F1 int e,
gsl sf hyperg 2F0, gsl sf hyperg 2F0 e, gsl sf hyperg 2F1, gsl sf hyperg 2F1 conj, gsl sf hyperg 2F1 conj e,
gsl sf hyperg 2F1 conj renorm , gsl sf hyperg 2F1 conj renorm e, gsl sf hyperg 2F1 e,
gsl sf hyperg 2F1 renorm, gsl sf hyperg 2F1 renorm e, gsl sf hyperg U , gsl sf hyperg U e,
gsl sf hyperg U e10 e, gsl sf hyperg U int, gsl sf hyperg U int e, gsl sf hyperg U int e10 e,
gsl sf hypot, gsl sf hypot e, gsl sf hzeta, gsl sf hzeta e, gsl sf laguerre 1, gsl sf laguerre 1 e,
gsl sf laguerre 2, gsl sf laguerre 2 e, gsl sf laguerre 3 , gsl sf laguerre 3 e, gsl sf laguerre n,
gsl sf laguerre n e, gsl sf lambert W0, gsl sf lambert W0 e, gsl sf lambert Wm1, gsl sf lambert Wm1 e,
gsl sf legendre array size, gsl sf legendre H3d, gsl sf legendre H3d 0, gsl sf legendre H3d 0 e,
gsl sf legendre H3d 1, gsl sf legendre H3d 1 e, gsl sf legendre H3d array, gsl sf legendre H3d e,
gsl sf legendre P1, gsl sf legendre P1 e, gsl sf legendre P2, gsl sf legendre P2 e, gsl sf legendre P3,

1
gsl sf legendre P3 e, gsl sf legendre Pl, gsl sf legendre Pl array, gsl sf legendre Pl e,

1
La biblioth`eque GSL

gsl sf legendre Plm, gsl sf legendre Plm array, gsl sf legendre Plm e, gsl sf legendre Q0,
gsl sf legendre Q0 e, gsl sf legendre Q1, gsl sf legendre Q1 e, gsl sf legendre Ql, gsl sf legendre Ql e,
gsl sf legendre sphPlm, gsl sf legendre sphPlm array, gsl sf legendre sphPlm e, gsl sf lnbeta,
gsl sf lnbeta e, gsl sf lnchoose, gsl sf lnchoose e, gsl sf lncosh, gsl sf lncosh e, gsl sf lndoublefact,
gsl sf lndoublefact e, gsl sf lnfact, gsl sf lnfact e, gsl sf lngamma, gsl sf lngamma complex e,
gsl sf lngamma e, gsl sf lngamma sgn e, gsl sf lnpoch, gsl sf lnpoch e, gsl sf lnpoch sgn e,
gsl sf lnsinh, gsl sf lnsinh e, gsl sf log, gsl sf log 1plusx, gsl sf log 1plusx e, gsl sf log 1plusx mx,
gsl sf log 1plusx mx e, gsl sf log abs, gsl sf log abs e , gsl sf log e, gsl sf log erfc, gsl sf log erfc e,
gsl sf multiply e, gsl sf multiply err e, gsl sf poch, gsl sf poch e, gsl sf pochrel, gsl sf pochrel e,
gsl sf polar to rect, gsl sf pow int, gsl sf pow int e, gsl sf psi, gsl sf psi 1 int, gsl sf psi 1 int e,
gsl sf psi 1piy, gsl sf psi 1piy e, gsl sf psi e, gsl sf psi int, gsl sf psi int e, gsl sf psi n,
gsl sf psi n e, gsl sf rect to polar, gsl sf Shi, gsl sf Shi e, gsl sf Si, gsl sf Si e, gsl sf sin,
gsl sf sin e, gsl sf sin err, gsl sf sin err e, gsl sf sinc, gsl sf sinc e, gsl sf synchrotron 1,
gsl sf synchrotron 1 e, gsl sf synchrotron 2, gsl sf synchrotron 2 e, gsl sf taylorcoeff,
gsl sf taylorcoeff e, gsl sf transport 2, gsl sf transport 2 e, gsl sf transport 3, gsl sf transport 3 e,
gsl sf transport 4, gsl sf transport 4 e, gsl sf transport 5, gsl sf transport 5 e, gsl sf zeta,
gsl sf zeta e, gsl sf zeta int, gsl sf zeta int e,
gsl siman solve

C.2.29 Tris
gsl sort, gsl sort index, gsl sort largest, gsl sort largest index, gsl sort smallest, gsl sort smallest index,
gsl sort vector, gsl sort vector index, gsl sort vector largest, gsl sort vector largest index,
gsl sort vector smallest, gsl sort vector smallest index, gsl heapsort, gsl heapsort index

C.2.30 Lissage
gsl spline alloc, gsl spline eval, gsl spline eval deriv, gsl spline eval deriv2, gsl spline eval deriv2 e,
gsl spline eval deriv e, gsl spline eval e, gsl spline eval integ, gsl spline eval integ e,
gsl spline free, gsl spline init

C.2.31 Statistique
gsl stats absdev, gsl stats absdev m, gsl stats covariance, gsl stats covariance m,
gsl stats kurtosis, gsl stats kurtosis m sd, gsl stats lag1 autocorrelation ,
gsl stats lag1 autocorrelation m, gsl stats max, gsl stats max index, gsl stats mean,
gsl stats median from sorted data, gsl stats min, gsl stats min index, gsl stats minmax,
gsl stats minmax index,
gsl stats quantile from sorted data, gsl stats sd, gsl stats sd m, gsl stats sd with fixed mean,
gsl stats skew, gsl stats skew m sd, gsl stats variance, gsl stats variance m
, gsl stats variance with fixed mean, gsl stats wabsdev, gsl stats wabsdev m, gsl stats wkurtosis,
gsl stats wkurtosis m sd, gsl stats wmean, gsl stats wsd, gsl stats wsd m, gsl stats wsd with fixed mean,
gsl stats wskew, gsl stats wskew m sd, gsl stats wvariance, gsl stats wvariance m,
gsl stats wvariance with fixed mean
gsl strerror

1
 C.2 Sous programmes

C.2.32 Transform´ees de Levin

gsl sum levin u accel, gsl sum levin u alloc, gsl sum levin u free, gsl sum levin utrunc accel,
gsl sum levin utrunc alloc, gsl sum levin utrunc free

C.2.33 Vecteurs
gsl vector add, gsl vector add constant, gsl vector alloc, gsl vector calloc,
gsl vector complex const imag, gsl vector complex const real , gsl vector complex imag,
gsl vector complex real, gsl vector const subvector , gsl vector const subvector with stride,
gsl vector const view array, gsl vector const view array with stride, gsl vector div, gsl vector fprintf,
gsl vector fread, gsl vector free, gsl vector fscanf, gsl vector fwrite, gsl vector get, gsl vector isnull,
gsl vector max, gsl vector max index, gsl vector memcpy, gsl vector min, gsl vector min index,
gsl vector minmax, gsl vector minmax index, gsl vector mul, gsl vector ptr, gsl vector ptr,
gsl vector reverse, gsl vector scale, gsl vector set, gsl vector set all, gsl vector set basis,
gsl vector set zero, gsl vector sub, gsl vector subvector, gsl vector subvector with stride,
gsl vector swap, gsl vector swap elements, gsl vector view array gsl vector view array with stride

Nous renvoyons au manuel pour une pr´esentation d´etaill´ee de

chacune de ces fonctions.

1
La biblioth`eque GSL

1
Bibliographie

[1] W. H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical

Recipes in C, Cambridge University Press (1992).

[2] http ://[Link]/gsl/

[3] http ://[Link]/lapack/

[4] http ://[Link]/blas/

[5] http ://[Link]/B673/

[6] http ://[Link]/[Link]/[Link]

1
BIBLIOGRAPHIE

1
Table des mati`eres

1 Int´egration et sommes discr`etes 3

1.1 Introduction.................................................................................................3
1.2 Les m´ethodes de Coˆtes ..........................................................................4
1.2.1 Trap`eze .........................................................................................5
1.2.2 Simpson...........................................................................................5
1.3 M´ethode de Romberg ..............................................................................6
1.4 M´ethodes de Gauss ..................................................................................6
1.5 Int´egrales multiples ..................................................................................9

2 Fonctions sp´eciales et ´evaluation de fonctions 11

2.1 Introduction..............................................................................................11
2.2 Fonctions transcendantes simples.........................................................12
2.3 Fonction Gamma.......................................................................................12
2.3.1 D´efinition et propri´et´es .........................................................12
2.3.2 Fonctions reli´ees : Ψ, B ............................................................14
2.4 Fonctions de Bessel..................................................................................15
2.5 Fonctions Hyperg´eom´etriques .............................................................17
2.5.1 Fonction Hyperg´eom´etrique Gaussienne ...............................17
2.5.2 Fonctions Hyperg´eom´etriques g´en´eralis´ees .....................18
2.6 Fonction erreur, exponentielle int´egrale ..............................................18
2.7 Conclusion.................................................................................................20

3 Racines d’´equations 21
3.1 Introduction..............................................................................................21
3.2 Dichotomie................................................................................................22
3.3 M´ethode de Ridder ................................................................................23
3.3.1 M´ethode de la position fausse ..................................................23
3.3.2 M´ethode de Ridder ....................................................................24
3.4 M´ethode de Brent ..................................................................................25
3.5 Newton-Raphson......................................................................................25
3.6 Racines de Polynˆomes ...........................................................................26
3.6.1 R´eduction polynomiale ..............................................................26
3.6.2 M´ethode de Laguerre ................................................................27

1
 TABLE DES MATIE` RES

4 Equations diff´erentielles 29
4.1 Introduction..............................................................................................29
4.2 D´efinitions ...............................................................................................30
4.3 M´ethodes d’int´egration `a pas s´epar´e ...........................................31
4.3.1 Introduction..................................................................................31
4.3.2 M´ethode d’Euler ........................................................................32
4.3.3 M´ethode RK explicites `a un point ........................................32
4.3.4 M´ethodes RK implicites `a un point ......................................33
4.3.5 M´ethodes RK explicites `a 2 points interm´ediaires ............33
4.3.6 M´ethodes RK explicites `a 3 points interm´ediaires ............33
4.3.7 Formule g´en´erale des m´ethodes RK explicites ..................34
4.4 M´ethode d’int´egration `a pas variables ............................................34
4.4.1 Introduction..................................................................................34
4.5 M´ethodes de Runge-Kutta “embarqu´ees” .........................................35
4.6 M´ethode de Bulirsh-Stoer .....................................................................36
4.7 Conclusion.................................................................................................37

5 Transform´ee de Fourier rapide 39

5.1 Introduction..............................................................................................39
5.2 Propri´et´es ..............................................................................................39
5.3 Discr´etisation de la transform´ee de Fourier ......................................42
5.3.1 E´ chantillonage ...........................................................................42
5.3.2 Transform´ee de Fourier discr`ete ............................................43
5.4 Transform´ee de Fourier rapide .............................................................44

6 Alg`ebre lin´eaire 47
6.1 Introduction..............................................................................................47
6.2 E´ limination de Gauss-Jordan ...............................................................49
6.2.1 Rappels sur les matrices.............................................................49
6.2.2 M´ethode sans pivot ...................................................................49
6.2.3 M´ethode avec pivot ...................................................................50
6.3 E´ limination gaussienne avec substitution ...........................................50
6.4 D´ecomposition LU ..................................................................................51
6.4.1 Principe.........................................................................................51
6.4.2 R´esolution d’un syst`eme lin´eaire .........................................52
6.5 Matrices creuses.......................................................................................53
6.5.1 Introduction..................................................................................53
6.5.2 Matrices tridiagonales.................................................................54
6.5.3 Formule de Sherman-Morison...................................................54
6.6 D´ecomposition de Choleski ...................................................................54
6.7 Conclusion.................................................................................................55

7 Analyse spectrale 57
7.1 Introduction..............................................................................................57
7.2 Propri´et´es des matrices .......................................................................58
7.3 M´ethodes directes ...................................................................................60
7.3.1 M´ethode de Jacobi ....................................................................60

1
 TABLE DES MATIE` RES

7.3.2 R´eduction de Householder ........................................................62

7.3.3 Algorithme QL.........................................................................64
7.3.4 Factorisation de Schur................................................................65
7.4 M´ethode it´eratives ...............................................................................66
7.4.1 M´ethodes des puissances ...........................................................66
7.4.2 M´ethode de Lanczo`s ................................................................67

8 Equations int´egrales 69
8.1 Introduction..............................................................................................69
8.2 Equation de Fredholm.............................................................................69
8.2.1 Equation de premi`ere esp`ece ..................................................69
8.2.2 Equation de seconde esp`ece ......................................................70
8.3 Equation de Volterra................................................................................71
8.3.1 Equation de premi`ere esp`ece ..................................................71
8.3.2 Equation de seconde esp`ece ......................................................71
8.4 Conclusion.................................................................................................71

9 Equations aux d´eriv´ees partielles 73

9.1 Introduction..............................................................................................73
9.2 Equations avec conditions aux fronti`eres ............................................76
9.2.1 Introduction..................................................................................76
9.2.2 Diff´erences finies ........................................................................76
9.2.3 M´ethodes matricielles ................................................................77
9.2.4 M´ethodes de relaxation .............................................................77
9.2.5 M´ethodes de Fourier .................................................................78
9.3 Equations avec conditions initiales........................................................80
9.3.1 Equations `a flux conservatif .....................................................80
9.3.2 Une approche na¨ıve....................................................................81
9.3.3 Crit`ere de Stabilit´e de Von Neumann ..................................81
9.3.4 M´ethode de Lax .........................................................................82
9.4 Conclusion.................................................................................................83

A Coordonn´ees hypersph´eriques 85

B Les biblioth`eques BLAS et Lapack 87

B.1 Introduction..............................................................................................87
B.2 Terminologie.............................................................................................88

C La biblioth`eque GSL 89
C.1 Introduction..............................................................................................90
C.2 Sous programmes.....................................................................................91
C.2.1 BLAS......................................................................................... 91
C.2.2 fonctions simples.........................................................................92
C.2.3 interfaces entre GSL et le BLAS.................................................92
C.2.4 Blocs..............................................................................................92
C.2.5 S´eries de Chebyshev. .................................................................92
C.2.6 Combinatoire................................................................................92

1
TABLE DES MATIE` RES

C.2.7 Complexes.....................................................................................92
C.2.8 Hankel...........................................................................................93
C.2.9 D´eriv´ees ....................................................................................93
C.2.10 Valeurs et vecteurs propres.......................................................93
C.2.11 Transform´ees de Fourier ...........................................................93
C.2.12 Ajustements..................................................................................93
C.2.13 Histogrammes..............................................................................94
C.2.14 IEEE.......................................................................................... 94
C.2.15 Int´egration ..................................................................................94
C.2.16 Interpolation.................................................................................94
C.2.17 Alg`ebre lin´eaire ........................................................................95
C.2.18 Matrices.........................................................................................95
C.2.19 Minimisation.................................................................................95
C.2.20 Monte Carlo..................................................................................95
C.2.21 Ajustements non-lin´eaires .........................................................96
C.2.22 Equations diff´erentielles ............................................................96
C.2.23 Permutation..................................................................................96
C.2.24 Polynˆomes ...................................................................................97
C.2.25 Puissances.....................................................................................97
C.2.26 Nombres al´eatoires ....................................................................97
C.2.27 Racines unidimensionnelles.......................................................97
C.2.28 Fonctions sp´eciales ....................................................................98
C.2.29 Tris..............................................................................................100
C.2.30 Lissage.........................................................................................100
C.2.31 Statistique...................................................................................100
C.2.32 Transform´ees de Levin ............................................................101
C.2.33 Vecteurs......................................................................................101