Méthodes Numériques et Optimisation .pdf



Nom original: Méthodes Numériques et Optimisation.pdf

Ce document au format PDF 1.5 a été généré par TeX / MiKTeX pdfTeX-1.40.11, et a été envoyé sur fichier-pdf.fr le 17/01/2012 à 23:21, depuis l'adresse IP 160.228.x.x. La présente page de téléchargement du fichier a été vue 5294 fois.
Taille du document: 406 Ko (44 pages).
Confidentialité: fichier public


Aperçu du document


Méthodes Numériques
et Optimisation
Supélec
17 janvier 2012

Cours de Jérôme Juillard édité par Daniel Barral

2

Table des matières
1 Présentation générale

5

2 Problèmes linéaires
2.1 Résolution des systèmes linéaires . . . . . . . . . . .
2.1.1 Exemples d’applications . . . . . . . . . . . .
2.1.2 Résolution des SL vs. inversion de matrice . .
2.1.3 Méthodes de résolution . . . . . . . . . . . .
2.2 Calcul du déterminant d’une matrice . . . . . . . . .
2.3 Problèmes de valeurs propres . . . . . . . . . . . . .
2.3.1 Exemples d’applications . . . . . . . . . . . .
2.3.2 Méthodes de calcul d’une ou de plusieurs v.p.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

3 Problèmes non linéaires
3.1 Illustrations . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Généralités . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 Différences entre problèmes linéaires et non linéaires
3.2.2 Vitesse et critères de convergence . . . . . . . . . . .
3.2.3 Applications contractantes et itérations du point fixe
3.2.4 D’où viens-je ? Qui suis-je ? Où vais-je ? . . . . . . .
3.3 Problèmes en dimension 1 . . . . . . . . . . . . . . . . . . .
3.3.1 Méthodes à convergence garantie . . . . . . . . . . .
3.3.2 Méthodes à convergence non garantie . . . . . . . .
3.4 Problèmes en dimension N . . . . . . . . . . . . . . . . . . .
3.4.1 Méthode de Newton-Raphson . . . . . . . . . . . . .
3.4.2 Méthodes de quasi-Newton (Broyden) . . . . . . . .
3.5 Accélération de la convergence . . . . . . . . . . . . . . . .
3.5.1 ∆2 d’Aitken . . . . . . . . . . . . . . . . . . . . . . .
3.5.2 Méthode de Steffensen . . . . . . . . . . . . . . . . .
3.5.3 Méthode du passage à la limite de Richardson . . .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

7
7
7
7
8
10
10
10
11

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

14
14
14
14
14
15
15
15
16
16
17
17
17
18
18
18
18

4 Optimisation
19
4.1 Optimisation à N dimensions . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.2 Optimisation à une dimension . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.3 Méthode du point fixe "simple" ou du gradient . . . . . . . . . . . . . . . . 20
3

4

TABLE DES MATIÈRES
4.4
4.5
4.6

4.7

Méthode de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Problème des moindres carrés linéaires . . . . . . . . . . . . . . . . . . . . .
4.5.1 Régularisation de Tikhonov . . . . . . . . . . . . . . . . . . . . . . .
Moindres carrés non linéaires . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6.1 Méthode de Gauss-Newton . . . . . . . . . . . . . . . . . . . . . . .
4.6.2 Méthode de Levenberg-Marquardt . . . . . . . . . . . . . . . . . . .
Problèmes généraux d’optimisation . . . . . . . . . . . . . . . . . . . . . . .
4.7.1 Méthode de Broyden . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.7.2 Formule de Sherman-Morrison . . . . . . . . . . . . . . . . . . . . .
4.7.3 Méthode du gradient conjugué (ou méthode des directions conjuguées)

5 Interpolation et extrapolation - Dérivation et intégration
5.1 Généralités et applications . . . . . . . . . . . . . . . . . . .
5.2 Problèmes en une dimension - Interpolation polynômiale . .
5.2.1 Polynômes de Lagrange . . . . . . . . . . . . . . . .
5.2.2 Système de Van der Monde . . . . . . . . . . . . . .
5.2.3 Formule barycentrique de Lagrange . . . . . . . . . .
5.2.4 Limites de l’interpolation polynômiale . . . . . . . .
5.2.5 Interpolation par les splines cubiques . . . . . . . . .
5.2.6 Méthodes à noyaux . . . . . . . . . . . . . . . . . . .
5.3 Interpolation à D dimensions . . . . . . . . . . . . . . . . .
5.4 Dérivation numérique . . . . . . . . . . . . . . . . . . . . .
5.4.1 Formules des différences finies . . . . . . . . . . . . .
5.4.2 Choix d’un pas optimal . . . . . . . . . . . . . . . .
5.5 Intégration numérique . . . . . . . . . . . . . . . . . . . . .
5.5.1 Intégration à une dimension . . . . . . . . . . . . . .
5.5.2 Intégration en D dimensions . . . . . . . . . . . . . .
6 Problèmes aux valeurs initiales et
6.1 Problèmes aux valeurs initiales .
6.1.1 Contexte de l’étude . . .
6.1.2 Notions élémentaires . . .
6.1.3 Méthodes de résolution .
6.2 Problèmes aux limites . . . . . .
6.2.1 Notations . . . . . . . . .
6.2.2 Cadre de l’étude . . . . .
6.2.3 Démarche générale . . . .
6.2.4 Problèmes stationnaires .
6.2.5 Problèmes transitoires . .

aux
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .
. . .

limites
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

numérique
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

21
21
22
22
22
23
23
23
24
24

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

26
26
27
27
27
28
28
28
29
29
30
30
31
31
31
33

.
.
.
.
.
.
.
.
.
.

35
35
35
35
37
38
38
39
40
41
44

Chapitre 1

Présentation générale
Les méthodes numériques sont d’abord des mathématiques appliquées :
– à la résolution d’équations (systèmes linéaires, non linéaires, EDO, EDP)
– à l’approximation de certaines quantités (intégrales, dérivées)
– à l’optimisation
Ces calculs doivent être menés avec un certain nombre de contraintes : il faut respecter
un certain budget en opérations et utiliser les opérations les moins coûteuses. La précision
doit être finie.
Exemple 1.1


On cherche à trouver la racine positive de x2 = x + 1. La racine s’écrit −1+2 5 : c’est la
forme "transcendante". Comment en donner une estimation numérique ? Par la suite :
u0 = a
1
a
un+1 = (un +
)
2
un
car on a lim un =
n→∞


−1+ 5
2

L’expression de la solution est différente de la solution (idéal platonicien) ainsi que de
la valeur approchée.
Le but des méthodes numériques est de trouver des approximations de solutions à des
vrais problèmes.
Exemple 1.2
Résolution de l’équation de la chaleur.
∇(k∇T ) = Q
On a des conditions aux limites simples (k uniforme ne dépend pas de la température, la source de chaleur est ponctuelle, T = T0 à l’infini.
R
Q(M )
Q
On aura T (r) = 4πkr
+ T0 . Le calcul de l’intégrale T (M 0 ) = Sources k4πd(M,M
0 ) dM
se fera avec des méthodes numériques.
5

6

CHAPITRE 1. PRÉSENTATION GÉNÉRALE

Exemple 1.3
Optimisation du rendement d’un moteur.
Quelles sont les dimensions (d1 ,...,dn ) qui maximisent P ?
On ne dispose pas d’expression analytique de P : chaque évaluation de P nécessite
une solution numérique. La question est de trouver un optimum avec un minimum
d’évaluations et de calculs possibles. On cherche à trouver une solution respectant des
contraintes d’encombrement également.
On peut faire tout (ou presque) avec des méthodes numériques, mais chaque méthode a
coût calculatoire. un bon numéricien est avare avec ses opérations. Par exemple, l’inversion
via Cramer se fait en O(N !) ou en O(N 4 ). L’inversion via LU se fait en O(N 3 ).
Une méthode ne peut pas être plus précise que l’ordinateur sur lequel on travaille
16
(10 + 1 = 1016 )
Une petite erreur peut avoir des conséquences importantes. Par exemple, on souhaite
(t)
calculer lim f (t+h)−f
. Sur l’ordinateur, on pense calculer le taux d’accroissement de f
h
h→0

ˆ
fˆ(t)
entre t et t+h, alors qu’on calcule f (t+h)−
où fˆ comporte une erreur d’arrondi par
h
rapport à f. Les deux erreurs en t+h et t ne sont pas nécessairement égales et on peut
arriver à des erreurs grossières.
L’objectif de ce cours est de devenir un utilisateur éclairé des méthodes numériques,
pas des informaticiens ni des mathématiciens. Il s’agit d’être capable de formaliser un
problème réel de manière standard, i.e. de reconnaître des classes de problèmes.
Il faut être conscient des implications numériques (coûts calculatoire, mémoire), et du
champ d’application des différentes méthodes. Il faut donc être capable de choisir avec
discernement l’outil le plus adapté à la résolution d’un problème donné.

Chapitre 2

Problèmes linéaires
2.1
2.1.1

Résolution des systèmes linéaires
Exemples d’applications

Exemple 2.1
Micro-capteur de pression
Sous l’effet d’un différentiel de pression, la membrane se déforme, et les résistances
R2 et R3 varient. Le signal mesuré est lié à la forme prise par la membrane. Connaissant
dP comment déterminer cette forme ?
Lorsque le pont à 4 résistances est déséquilibré, le voltmètre détecte un différentiel
de tension dV que l’on espère proportionnel à dP.
Discrétisation du problème :
La membrane est assimilée à un réseau de ressorts connue. Chaque noeud du réseau
ne voit que ses plus proches voisins, les déplacements des noeuds du bord sont fixés à
zéro, chaque noeud est soumis à une force proportionnelle dP.
pi = K

X

(xi − xj )

j∈Vi

2.1.2

Résolution des SL vs. inversion de matrice

Le problème que l’on cherche à résoudre est l’équation Ax = b où A est à coefficients
réels, carrée et inversible.
Formellement, la solution va s’écrire x = A−1 b. Comment fait-on pour calculer A−1 ?
La meilleure méthode est la suivante : on pose v1 ,...,vn les vecteurs colonnes inconnus
de A−1 . On pose e1 ,...,en les vecteurs unités (colonnes de In ). Puis on résout le système :



 Av1 = e1

...



 Av = e
n
n

7

8

CHAPITRE 2. PROBLÈMES LINÉAIRES

En général, l’inversion de matrice n’est pas économique. On préfère résoudre le système
linéaire associé.

2.1.3

Méthodes de résolution

Les méthodes directes de résolution de systèmes linéaires permettent d’arriver à la
solution du problème en un nombre fini d’opérations.
On cherche x tel que Ax = b. On résout en fait (A + δA)ˆ
x = b + δb. La solution est
entachée d’erreur : x
ˆ = x + δx.
On a Ax+Aδx+δAx = b+δb d’où δx = A−1 (δb−δAˆ
x). Ainsi, ||δx|| = ||A−1 (δb−δAˆ
x)||.
−1
−1
On a donc ||δx|| ≤ ||A ||.||δb − δAˆ
x|| ≤ ||A ||(||δb|| + ||δA||.||ˆ
x||). L’erreur relative est
donc majorée par :
||δx||
||δb||
||δb||
||δA||
≤ ||A−1 ||(
+ ||δA||) = ||A||.||A−1 ||.(
+
)
||x||
||ˆ
x||
||A||.||ˆ
x||
||A||
Définition 2.1
On note Cond(A) = ||A||.||A−1 || le conditionnement de la matrice A, où ||A||2 =
σmax (A) =plus grande valeur singulière de A=Racine carrée de la plus grande vp de
AT A, et ||A−1 || = σmin1 (A) =plus petite valeur singulière de A.
On a l’inégalité de norme d’algèbre : ||BC|| ≤ ||B||.||C|| avec égalité si l’une des deux
matrices est orthogonale.
(A)
Ainsi, Cond(A) = σσmax
. Un système est dit mal conditionné si Cond(A) ' 1016 : A
min (A)
est "presque non inversible".
Factorisation LU
L’objectif est de décomposer A sous la forme A = L.U où L est triangulaire inférieure
et U est triangulaire supérieure.
Connaissant L et U, résoudre Ax = b se fait en 2 étapes :
– On résout Ly = b par substitution de proche en proche : O(N 2 ) opérations
– On résout U x = y : O(N 2 ) opérations
La factorisation LU s’obtient par l’algorithme de Crout. Celle-ci est unique si l’on
impose les valeurs de L sur la diagonale, par exemple : lii = 1. Cet algorithme est une
boucle sur les colonnes
Pour la colonne 1 :
u11 .1 = a11
1
u11 .l21 = a21 ⇒ l21 =
.a21
u11
...
Tous les li 1 (i=1..N) sont déterminés.
Pour la colonne 2 :
...
Par cet algorithme on obtient la décomposition LU en O( 32 N 3 ).

2.1. RÉSOLUTION DES SYSTÈMES LINÉAIRES

9

Propriété 2.1
D’après l’inégalité de norme d’algèbre, on a Cond(A) ≤ Cond(L)Cond(U ).
Par conséquent, si A est mal conditionnée, le problème via LU l’est encore plus.
Factorisation de Cholesky
Il s’agit de la décomposition LU adaptée pour les matrices symétriques définies positives.
Si A est symétrique définie positive, ∃!C ∈ Mn ( )/A = C T C avec C triangulaire
supérieure.

R

Factorisation QR
Cette méthode s’utilise si A est mal conditionnée.

R

∃!(Q, R) ∈ Mn ( ), A = QR
avec Q orthogonale et R triangulaire supérieure.
La résolution de Ax=b se fait alors en deux étapes :
– Trouver y tel que Qy = b ⇐⇒ y = QT b (en O(N 2 ))
– Trouver x tel que Rx = y (en O(N 2 )).
La factorisation sous forme QR nécessite O( 43 N 3 ) opérations (et est donc plus couteuse
que LU). Cependant, on a :
Cond(A) = Cond(R)
Le conditionnement est donc amélioré.
Méthodes itératives
L’objectif est de construire une suite xn convergeant vers la solution de Ax=b.
On écrit A = A1 + A2 où A1 est facile à inverser. On a alors :
Ax = b ⇔ A1 x + A2 x = b ⇔ x = A−1
1 (−A2 x + b)
Ainsi, si xn+1 = A−1
1 (−A2 xn + b) converge, c’est vers la solution de Ax = b.
La méthode de Jacobi suggère de poser A1 = diag(A) = diag(a1 , ..., an ).
On obtient alors l’algorithme itératif suivant :
for n=1 to Nmax
for i=1 to N
P
n
xn+1
= a1ii (− N
j=1,j6=i aij xj + bi )
i
next i
next n
La méthode de Gauss-Seidel est la suivante, et converge deux fois plus vite environ
que Jacobi :

10

CHAPITRE 2. PROBLÈMES LINÉAIRES

for n=1 to Nmax
for i=1 to N
P
P
n+1
n
xn+1
= a1ii (− i−1
− N
j=i+1 aij xj + bi )
i
j=1 aij xj
next i
next n
En pratique, on préférera la méthode de sur-relaxation. On pose :
xn+1 = (1 − ω)xn + ωxn+1
Gauss−Seidel
où ω ∈]1; 2[.
Celle-ci permet d’améliorer la vitesse den convergence de Gauss-Seidel.
Si l’on souhaite atteindre un résidu ||Ax||b||−b|| de l’ordre de 10−p , sur le problème de la
membrane,
O(pN) itérations sont nécéssaires. Avec Jacobi, 2 fois moins que Gauss-Seidel,

O(p N ) avec sur-relaxation.

2.2

Calcul du déterminant d’une matrice

Pour calculer le déterminant de la matrice, on effectue la décomposition de cette
matrice (par exemple sous factorisation LU), puis on calcule par produit : Det(A) =
Q
det(L)det(U ) = ni=1 uii

2.3

Problèmes de valeurs propres

Le but est de résoudre l’équation Ax = λx d’inconnues x et λ.

2.3.1

Exemples d’applications

Exemple 2.2
Les pages Google sont classées selon leur pagerank : PageRank(i)=probabilité d’arriver
sur la page i en surfant au hasard, au bout d’un temps infini.
Matrice d’interconnexion du www : A est définie par Aij = L1j si la page i pointe
vers j, où L(j) est le nombre de liens pointant vers j.
On appelle matrice de transition la matrice T = αU + (1 − α)A, avec α facteur de
lassitude (probabilité de cliquer au hasard). La ième coordonnée du 1er vecteur propre
de T correspond au PageRank de la page i.
Exemple 2.3
Modes propres des structures mécaniques
En régime harmonique, de nombreuses structures mécaniques peuvent être modélisées par :

Ku − ω 2 M u = f

2.3. PROBLÈMES DE VALEURS PROPRES

11

où K est la matrice de raideur de taille N, u le vecteur déplacement, ω la pulsation,
M la matrice de masse et f le vecteur des efforts.
Il existe N modes propres de vibration correspondant à des formes u1 , ..., uN et des
pulsations ω1 , ..., ωN :
Kui − ωi2 M ui = 0
c’est à dire que chaque mode de vibration est un vecteur propre de M −1 K.

2.3.2

Méthodes de calcul d’une ou de plusieurs v.p.

Méthode de la puissance
La méthode est la suivante : on choisit un vecteur unitaire aléatoire bO .
Ab0
Abn
, et on itère le processus : bn+1 = ||Ab
On calcule b1 = ||Ab
n ||
0 ||
Si toutes les valeurs propres de A sont positives, réelles, distinctes, on note λ1 , ...λN
avec λi > λi+1 . On note v1 , ..., vN les vecteurs propres unitaires associés.
P
On a Ab0 = N
k=1 ak λk vk donc par une récurrence immédiate :
N
X

An b0 =

ak λnk vk = λ1

k=1

N
X
k=1

ak

λk n
vk
λ1

Ainsi, An b0 ≈ λn1 a1 v1 d’où ||An b0 || ≈ |a1 |.
Cette méthode converge d’autant plus vite que le rapport | λλ12 | est petit. λˆ1 est donné

n ||
(si λ1 > 0) par ||Ab
||bn || .
Pour trouver les valeurs propres de moindre magnitude, on utilise la technique dite de

déflation : on forme A0 = A − λ1

vˆ1 vˆ1T
vˆT vˆ
1

. A’ a alors les mêmes vecteurs propres que A mais

1

v1 est associé à la valeur propre 0.
Itérations QR

On pose A0 = A.
On cherche la factorisation QR de A (cela se fait en O( 43 N 3 )), et on peut noter A0 =
Q0 R0 . On construit A1 = R0 Q0 .
Lorsqu’elle converge, la suite An tend vers une matrice triangulaire par blocs (forme
de Schur) :
– Triangulaire supérieure si toutes les valeurs propres sont réelles
– Sinon des blocs, dont la taille dépend de la multiplicité des valeurs propres complexes,
"dépassent" sous la diagonale.
– Les vecteurs propres se déduisent par des opérations simples (O(N 2 )) sur les colonnes
de pn ;
– Converge au bout de O(N ) itérations ;
– Le coût total du calcul des valeurs propres devrait-être en O(N 4 ). En fait, on peut
se ramener à du O(N 3 ).

12

CHAPITRE 2. PROBLÈMES LINÉAIRES

Idée : mettre A sous la forme de Hessenberg :
– Trangulaire supérieure avec des éléments non-nuls sur la première sous-diagonale.
A = QT HQ

O(N 3 )

– On applique ensuite la méthode QR → RQ à H. La factorisation QR d’une matrice
de Hessenberg coûte O(N 2 ) opérations et la forme de Hessenberg est conservée par
la transformation.
Calcul des racines d’un polynôme
On veut calculer les racines de P (λ) = c0 + c1 λ + ... + cn−1 λn−1 + λn .
On constate que P (λ) est le polynôme caractéristique de la matrice compagnon :


0 0 ···
1 0 · · ·


..

.
0 1
C=
..
 ..
.
. 0

 .. .. . .
. .
.
0 0 ···

0
0

0
0









0 cn−2 

0 0
. . ..
. .
1
0



c0
c1 

c2
..
.

1 cn−1

Méthode d’Arnoldi
Cette méthode permet de trouver des approximations de quelques valeurs propres et
les vecteurs propres associés.
L’idée est de ne faire qu’une factorisation partielle sous forme de Hessenberg : A =
P HP T ⇐⇒ AP = P H.
On ne s’intéresse qu’aux n << N premières colonnes de AP .
˜ n Pn+1 où Pn est la matrice constituée des n premières colonnes de
On note APn = H
˜
P et Hn est une matrice de n+1 lignes et n colonnes extraite de H.
On note Pn = [p1 , ..., pn ]. On part de p1 unitaire aléatoire, et on calcule v = Ap1
(O(N 2 ) opérations). On cherche alors à exprimer v dans une base (p1 , p2 ) orthonormée :
– On déduit h11 = v T p1
– On déduit h21 = ||v − h11 p1 ||
11 p1
– On déduit p2 = v−h
h21
Le processus (orthogonalisation de Gram-Schmidt est répété n fois).
Le coût de la nième itération est celui du produit matrice-vecteur, soit O(N 2 ) opérations. Le calcul d’une nouvelle colonne de P et de H se fait en O(n).
Propriété 2.2
Les valeurs propres de Hn tendent vers les valeurs propres de plus grande magnitude
de A.

2.3. PROBLÈMES DE VALEURS PROPRES

13

Itération inverse de la puissance
Cette méthode permet de raffiner les estimations des valeurs et vecteurs propres trouvés.
ˆ 1 de λ1 .
On suppose qu’on dispose d’une approximation λ
ˆ 1 IN )−1 :
On applique alors la méthode de la puissance à (A − λ

ˆ 1 IN )v = bn

 (A − λ



bn+1 =

v
||v||

On doit donc résoudre un système linéaire à chaque itération. Cela est peu économique.
ˆ 1 est valeur propre de A − λ
ˆ 1 et 1 est valeur
Si λ1 est valeur propre de A, λ1 − λ
ˆ1
λ1 − λ
−1
ˆ
propre de (A − λ1 ) .

Chapitre 3

Problèmes non linéaires
3.1

Illustrations

R

R

On considère une application f : N −→ N .
On cherche à trouver x ∈ N tel que f (x) = 0.

R

3.2
3.2.1

Généralités
Différences entre problèmes linéaires et non linéaires

Existence de la solution
Pour le problème linéaire, l’existence de la solution de Ax = b est toujours vérifiée si
A est inversible.
Cependant, pour le problème non linéaire, f (x) = 0 n’admet pas toujours de solution.
Multiplicité de la solution
Dans le cas linéaire, la solution est unique si elle existe.
Dans le cas non linéaires, les solutions peuvent être multiples, voire infinies en nombre.
Méthode directe de calcul de la solution
Une telle méthode existe dans le cas linéaire mais n’existe pas en général dans le cas
non linéaire.
Les problèmes non linéaires ne peuvent en général être résolus que de manière itérative.
Le but est donc de choisir la méthode appropriée.

3.2.2

Vitesse et critères de convergence

Soit une suite xn tendant vers x∞ . On parle de convergence :
−x∞ ||
−→ C
– linéaire si ∃C ∈]0, 1[ t.q. ||x||xn+1
n −x∞ ||
– Super-linéaire si C = 0.
14

3.3. PROBLÈMES EN DIMENSION 1

15

n+1 −x∞ ||
0
– super-linéaire d’ordre µ si ∃µ, C 0 > 0 t.q. ||x
||xn −x∞ ||µ −→ C . La convergence est dite
quadratique si µ = 2 et cubique si µ = 3.
Un compromis doit être souvent trouvé entre vitesse de convergence et robustesse.

3.2.3

Applications contractantes et itérations du point fixe

Définition 3.1
(y)||
<
On dit qu’une application g est contractante si et seulement si ∃C ∈]0, 1[ tq ∀(x, y), ||f (x)−f
||x−y||
C.

Théorème 3.1
Si g est contractante, elle admet un point fixe unique vers lequel la suite xn+1 = g(xn )
converge. La vitesse de convergence d’une intération du point fixe dépend des dérivées
de g.
Théorème 3.2

R

R

Soit g : −→ contractante et x0 son unique point fixe. On définit la suite (xn ) par
xn+1 = g(xn ).
Si |g 0 (x0 )| < 1, alors la convergence de (xn )n est linéaire.
Si |g 0 (x0 )| = 0, alors la convergence de (xn )n est quadratique.

3.2.4

D’où viens-je ? Qui suis-je ? Où vais-je ?

D’où viens-je ?
On initialise la méthode : "Posons x1 = ... et n = 1."
Qui suis-je ?
On développe le corps de la méthode : "faire xn+1 = g(xn , ...). Incrémenter n."
Où vais-je ?
On termine la métode : "Continuer tant que..."
On peut se fixer des critères de non échec : n < Niter_max , ainsi que des critères de
réussite, par exemple :
−xn−1 ||
– ||x||xnn||+||x
>
n−1 ||



3.3

||xn −xn−1 ||
>
||x1 |||
||f (xn )||
||f (x1 )|| >



Problèmes en dimension 1

Il existe des méthodes, dans le cas où la solution existe, où la convergence est garantie :
la méthode de dichotomie (ou bissection), la méthode de la fausse position, ou la méthode

16

CHAPITRE 3. PROBLÈMES NON LINÉAIRES

de Brent. Ces méthodes nécéssitent un encadrement de la solution, à savoir des points a
et b tels que f (a)f (b) < 0.

3.3.1

Méthodes à convergence garantie

Méthode de dichotomie
On part de a,b tels que f (a)f (b) < 0.
a+b
On calcule c = a+b
2 et f (c) = f ( 2 ).
Si f (a)f (c) < 0 alors on remplace (a,b) par (a,c) et on itère. Sinon, on remplace (a,b)
par (c,b) et on itère.
Méthode de la fausse position
On part de a,b tels que f (a)f (b) < 0.
b−a
On calcule c = a − f (a) f (b)−f
(a) et f (c).
Si f (a)f (c) < 0 alors on remplace (a,b) par (a,c) et on itère. Sinon, on remplace (a,b)
par (c,b) et on itère.

3.3.2

Méthodes à convergence non garantie

Méthode ne nécessitant pas la dérivée : méthode du point fixe simple
Elle consiste à construire g une fonction telle que x∞ = g(x∞ ) ⇐⇒ f (x∞ ) = 0.
Première idée : g(x) = x − λf (x). Le paramètre λ doit être choisi pour que g soit
localement contractante. La convergence sera linéaire sauf si f 0 (x1∞ ) = λ.
Méthode de Newton
Deuxième idée : on pose g(x) = x − ff0(x)
(x)
La convergence est alors quadratique. La méthode demande deux évaluations de fonctions par itération. Cependant, on doit pouvoir calculer la dérivée. Il est possible de relaxer
la méthode mais les coûts augmentent et la vitesse n’est plus quadratique.
Le rayon de convergence est souvent moins bon que pour la première idée.
Cependant, cette méthode ne marche pas toujours bien, par exemple lorsque la dérivée
de f s’annule au voisinage du point recherché.
Méthode de Newton relaxée
A chaque itération, on choisit λn ∈ [ , 1], et on fait :
xn+1 = xn − λn

f (xn )
f 0 (xn )

λn sera choisi par exemple en partant de λ1 = 1 puis λk+1 =
|f (xn+1 )| < |f (xn )|.

λk
2

de façon à avoir

3.4. PROBLÈMES EN DIMENSION N

17

Cette méthode peut être une bonne méthode dans le cas où f n’est pas très couteuse
à évaluer. Elle n’offre plus de vitesse quadratique.

Méthode de la sécante
Cette méthode consiste à construire une approximation λˆn de
des itérations.
λˆn =

1
f 0 (x∞ )

au fur et à mesure

xn − xn−1
xn − xn−1
⇒ xn+1 = xn −
f (xn )
f (xn ) − f (xn−1 )
f (xn ) − f (xn−1 )

L’ordre de convergence super-linéaire de cette méthode est le nombre d’or (µ = 1, 618).
Elle a pour avantage sur la méthode de Newton de ne pas évaluer la dérivée, et ne nécessite
qu’une seule évaluation de f par itération.
Deux itérations de sécante donnent un ordre de convergence µ > 2 pour autant d’évaluations de fonctions que Newton. Cependant, le problème d’arrondi près de la solution
est plus complexe à régler.
La méthode de Müller est également appliquable (méthode d’approximation de la
courbe par une parabole).

3.4

Problèmes en dimension N

Idée intuitive : à l’instar de la dimension 1, on pose g(x) = x − λf (x).

3.4.1

Méthode de Newton-Raphson

k
Cette méthode consiste à résoudre : J(xn )δxn = −f (xn ) où (Jkl )n = ∂f
∂xl (x = xn )
(matrice Jacobienne), puis à mettre à jour le point courant : xn+1 = xn + δn .
Cette méthode présente les mêmes défauts/qualités que Newton mais demande N 2 + 1
évaluations de fonction par itération et la résolution d’un système linéaire, car xn+1 =
xn − Jn−1 f (xn ).
De plus, si cond(Jn ) >> 1, on aura un problème de convergence.

3.4.2

Méthodes de quasi-Newton (Broyden)

Le principe est d’étendre la méthode de la sécante à la dimension quelconque. On
construit au fur et à mesure une approximation de la Jacobienne (ou mieux, de son inverse)
se raffinant au fur et à mesure de la convergence.
Cette méthode présente les mêmes qualités/défauts que la méthode de la sécante,
elle nécéssite N évaluations de fonction par itération (si l’on approxime l’inverse de la
Jacobienne) ainsi qu’une multiplication matrice-vecteur.

18

CHAPITRE 3. PROBLÈMES NON LINÉAIRES

3.5
3.5.1

Accélération de la convergence
∆2 d’Aitken

Soit une suite scalaire un convergeant linéairement vers u∞ .
Propriété 3.1
La suite un se comporte comme u∞ + λ1 C1n + o(C1n ) pour n suffisamment grand
On a alors λ1 et C1 inconnus à priori. On note ∆n = un+1 − un .
Propriété 3.2
La suite u
˜n = un −

∆2n
∆n+1 −∆n

converge plus vite vers u∞ que un .

Si u
˜n converge linéairement, la méthode du ∆2 peut être appliquée de manière itérée :
on élimine successivement les termes en C2 , C3 ...
Cette méthode présente des risques d’instabilité numérique car on divise par des très
petites quantités.

3.5.2

Méthode de Steffensen

Il s’agit d’une méthode de Aitken "améliorée" appliquée aux méthodes de point fixe.
Propriété 3.3
2

(g(vn )−vn )
Si un+1 = g(un ) converge linéairement, alors la suite vn+1 = vn − (g(g(vn ))−g(v
n ))−(g(vn )−vn )
converge quadratiquement.

3.5.3

Méthode du passage à la limite de Richardson

Cette méthode nécessite plus d’a priori (connaissance des Ci . Elle est cepedant plus
simple et tout autant efficace que la méthode d’Aitken.
On suppose un = u∞ + λ1 C1n + o(C1n ).
On a alors un+1 = u∞ + λ1 C1n+1 + o(C1n+1 ).
Si on connait C1 , on a C1 un −un+1 = (C1 −1)u∞ +o(C1n ), d’où la proposition suivante :
Propriété 3.4
La suite vn =

C1 un −un+1
C1 −1

tend vers u∞ plus vite de un .

Cette méthode est utilisée pour l’intégration, la dérivation et la résolution d’EDO.

Chapitre 4

Optimisation
4.1

Optimisation à N dimensions

Le problème standard d’optimisation est le suivant : on souhaite minimiser une fonction

J:

RN −→ R. Il s’agit donc de trouver xˆ = argminJ(x).

Maximiser J, c’est minimiser -J, on se ramène toujours au même problème au final.
On appelle x
ˆ le minimiseur de J, J(ˆ
x) le minimum de J, J le critère ou le coût.
Si ces quantités existent, on définit les notions suivantes :
Définition 4.1
On appelle gradient de J la quantité :
grad(J) = gJ =

∂J
∂x

Définition 4.2
On appelle matrice Hessienne de J, ou Hessien de J, la matrice
H = (H)ij
définie par :

H=

∂J
∂x

(H)ij =

∂J
∂xi ∂xj

Cette matrice est souvent symétrique car les fonctions étudiées sont au moins de
classe C2 (théorème de Schwartz).
L’objectif est de trouver des méthodes permettant de parvenir à un minimiseur local
de J.
Un exemple d’application classique est l’identification de système : on cherche à identifier les paramètres réels yi d’un système par une fonction de transfert yˆi . Par la méthode
des moindres carrés, on cherchera à minimiser la fonction de coût suivante :
19

20

CHAPITRE 4. OPTIMISATION

J(p) =

n
X

(yi − yˆi (p))2

i=0

On a les conditions nécéssaires d’optimalité suivantes :
Propriété 4.1
Si l’optimalité est atteinte en x
ˆ, on a nécéssairement :
(1) grad(J(ˆ
x))=0
(2) H est définie non négative en x
ˆ
On a la condition suffisante d’optimalité suivante :
Propriété 4.2
On a l’optimalité si (1) et (2) sont vérifiées simultanément :
(1) grad(J(ˆ
x))=0
(2) H est définie positive (strictement) en x
ˆ
Les méthodes d’optimisation dans
sations à une dimension.

4.2

RN se ramènent souvent à une succession d’optimi-

Optimisation à une dimension

R

R

Le problème est de chercher un minimum d’une fonction f : −→ . Nécéssairement,
x
ˆ vérifie f 0 (ˆ
x) = 0.
– Si on peut calculer f 0 , le problème revient à trouver le ou les zéros de f 0 . On utilise
donc la dichotomie, la fausse position, la méthode de la sécante ou du point fixe.
– Si on ne peut pas calculer f 0 , on doit faire des hypothèses supplémentaires sur f et
disposer d’un encadrement du minimum.
On dispose de deux points a,b tels que x
ˆ ∈]a, b[ et f soit minimale sur ]a, b[.
Si f est minimale sur ]a, b[, elle est décroissante à gauche de x
ˆ et croissante à droite.
L’algorithme est le suivant :
– On choisit deux points x− et x+ sur ]a, b[.
– Si f (x− ) < f (x+ ), alors x
ˆ ∈]a, x+ [ : on remplace (a, b) par (a, x+ )
– Sinon, on remplace (a, b) par (x− , b)

4.3

Méthode du point fixe "simple" ou du gradient

Pour annuler le gradient g(x), on construit une suite :
xi+1 = xi − λi g(xi )
.
Comme le but est de minimiser J, on choisit λi > 0 : on se déplace dans la direction
locale de la pente.

4.4. MÉTHODE DE NEWTON

21

Cette méthode est extrêmement robuste. Cependant, la convergence est très lente
dans les vallées du critère (la courbure est très faible dans une direction, et très forte dans
l’autre). Le Hessien est très mal conditionné à ces endroits.
Choix des λi Le choix d’un pas constant s’avère être une très mauvaise option. On
utilise la "steepest descent" (descente plus profonde).
On choisit à chaque étape le λ qui minimise la quantité :
f (λ) = J(xi − λg(xi ))
On utilise parfois la méthode du gradient conjugué, qui présente un faible coût calculatoire et de stockage.

4.4

Méthode de Newton

Il s’agit de la méthode de Newton-Raphson appliquée à l’optimisation, on construit la
suite suivante :
xi+1 = xi − H −1 (xi )g(xi )
Cette méthode assure une convergence quadratique, mais est moins robuste que la
méthode du gradient. Parfois, elle arrive à un maximum local et non à un minimum, dans
le cas ou le Hessien n’est pas défini positif. De plus, elle nécessite de calculer et de stocker
O(N 2 ) grandeurs et de résoudre un système linéaire.
On utilise également la méthode de quasi-Newton : on construit une approximation du
Hessien ou de son inverse.

4.5

Problème des moindres carrés linéaires

R

R

On cherche à minimiser J(x) = ||y − f (x)||22 , avec x ∈ n et y, f (x) ∈ N avec N ≥ n.
Si f (x) = F x, avec F ∈ MN,n ( ), on parle de problème des moindres carrés linéaires.
Dans ce cas, on cherche x
ˆ minimisant la quantité :

R

||y − F x||2 = (y − F x)T (y − F x) = y T y − 2y T F x + xT F T F x
On a donc :
= 0 ⇐⇒ 2F T F x − 2F T y = 0 ⇐⇒ x = (F T F )−1 F T y est la solution d’un système
linéaire.
F T F est une matrice carrée n,n définie positive, et l’équation F T F x = F T y est appelée
équation normale.
Souvent, les équations normales sont mal conditionnées :
∂J
∂x

cond(F T F ) = (cond(F ))2

22

CHAPITRE 4. OPTIMISATION

On effectue une décomposition QR de F, et on régularise le problème.
Pour une matrice F ∈ MN,n ( ), ∃Q ∈ MN,n ( ) et R ∈ Mn ( ) / F = QR telles que
T
Q Q = In et R soit triangulaire supérieure.
On a alors :

R

R

R

||Ax − b||2 = ||QRx − b||2 = ||Q(Rx − QT b)||2 = ||Rx − QT b||2
d’où :
x = R−1 QT b
Au prix d’une décomposition QR, on ne dégrade plus le conditionnement du système.

4.5.1

Régularisation de Tikhonov

Cette méthode consiste, plutôt que de minimiser ||y − F x||2 , à poser Jλ (x) = ||y −
F x||2 + λ||x||2 , λ > 0. Le terme en λ correspond à des contraintes a priori sur x.
Le minimiseur de Jλ est alors x
ˆλ = (F T F + λIn )−1 F T y.
Par la décomposition QR, on a :
F T F = RT QT QR = RT R
Si σ1 > ... > σn sont les valeurs propres positives de F T F , alors σ1 + λ > ... > σn + λ
sont celles de F T F + In . On a alors :
cond(F T F + λIn ) =

σ1 + λ
σ1
<
σn + λ
σn

On a donc réduit le conditionnement du système.

4.6

Moindres carrés non linéaires

La solution se construit dans le cas linéaires de façon itérative (méthodes de GaussNewton/Levenberg-Marquardt). A chaque nouvelle itération, on sera amené à résoudre les
équations normales.
Le problème est toujours le même : on veut trouver x
ˆ = argmin(||y − f (x)||2 ) =
2
argmin(||e(x)|| ).
P
2
On a alors J(x) = N
i=1 (||ei (x)||) .
La méthode de Newton consiste à poser xi+1 = xi − H −1 (xi )g(xi ).

4.6.1

Méthode de Gauss-Newton

La méthode de Gauss-Newton construit une approximation du Hessien tenant compte
de la forme particulière de J.

4.7. PROBLÈMES GÉNÉRAUX D’OPTIMISATION

23

On a :
g(x) =

∂J
∂x

=

PN

∂ei
i=1 2ei ∂x

= 2GT e, où G est la jacobienne de e.

Plutôt qu’utiliser H, avec :
Hkl =

N
X

∂ei ∂ei
∂ 2 ei
2
+ 2ei
∂xk ∂xl
∂xk ∂xj

i=1

!

on le remplace par son approximation :
a
Hkl
=


N
X
∂ei ∂ei

2

i=1

∂xk ∂xl

⇐⇒ H a = 2GT G

On pose donc : xi+1 = xi − (H a )−1 (xp )g(xp ) = xp − (GT G)−1 GT e(xp )
En pratique, à chaque itération, on trouve δx tel que GT Gδx = GT e, et on fait xp+1 =
xp − δx.

4.6.2

Méthode de Levenberg-Marquardt

Cette méthode apporte une correction supplémentaire à la méthode de Gauss-Newton :
cette fois-ci, on approxime le hessien par :
a
a + λI = 2GT G + λI , avec λ > 0 bien choisi.
HLM
= HGN
n
n

A chaque itération, on régularise les équations normales. On construit la suite :
xp+1 = xp − λg(xp )
.

4.7
4.7.1

Problèmes généraux d’optimisation
Méthode de Broyden

Les méthodes de quasi-Newton (BFGS, DFP) dérivent de la méthode de Broyden.
Cette méthode est l’extension de la méthode de la sécante en N dimensions.
On construit la suite :
xp = xp−1 + δx
ˆ p δx = −g(xp−1 ).
avec δx vérifiant H
ˆ
Pour que Hp soit une approximation du Hessien, on impose :
ˆ p (xp−1 − xp−2 ) = g(xp−1 ) − g(xp−2 )
H

24

CHAPITRE 4. OPTIMISATION

On a N 2 inconnues pour seulement N équations : la méthode de Broyden consiste
à ajouter N(N-1) contraites de façon à ce que le problème soit soluble. BFGS et DFP
ˆ p à être défini positif. La méthode de Broyden classique impose :
contraignent H
ˆ pu = H
ˆ p−1 u.
∀u⊥d = (xp−1 − xp−2 ), on a H


ˆp = H
ˆ p−1 I −
Ainsi, H

4.7.2

ddT
dT d



+

(f (xp−1 )−f (xp−2 )
dT d

Formule de Sherman-Morrison

On utilisera souvent la formule suivante :
Propriété 4.3
(A + uv T )−1 = A−1 −

4.7.3

A−1 uv T A−1
1+v T A−1 u

Méthode du gradient conjugué (ou méthode des directions conjuguées)

La convergence de cette méthode est linéaire mais plus rapide que le gradient simple.
Elle ne nécéssite que le stockage et calcul de O(N) valeurs.
Cette méthode peut être utilisée pour résoudre des systèmes linéaires de la forme Ax=b
où A est symétrique définie positive en N itérations et N produits matrice-vecteurs.
Application à la résolution de systèmes
Etant donnée une matrice A symétrique définie positive, on a :
x
ˆ solution de Ax = b ⇐⇒ x
ˆ minimise la forme quadratique J(x) = 21 xT Ax − bT x.
Initialisation :
– x1 = 0
– g(x1 ) = g1 = −b
– d1 = g1
Itérations :
i+1 )
– On construit xi+1 = xi − αi di et ∂J(x
=⇒ αi =
∂αi
– On détermine gi+1 = Axi+1 − b = gi − αi Ai di

dT
i
dT
i Adi

– On construit di+1 = gi+1 + βi di avec βi choisi tel que dTi+1 Adi = 0 soit βi = −
Gradient conjugué pour l’optimisation en grandes dimensions
Cette fois-ci, les itérations se font sous la forme :
– On trouve αi = argmin(J(xi − αdi et on déduit xi+1 = xi − αi di
– On calcule g(xi+1 ) = gi+1

T Ad
gi+1
i
dT
i Adi

4.7. PROBLÈMES GÉNÉRAUX D’OPTIMISATION
– On calcule di+1 = gi+1 −

T (g −g
gi+1
i
i+1 )
T
di+1 (gi −gi+1 )

25

Chapitre 5

Interpolation et extrapolation Dérivation et intégration
numérique
5.1

Généralités et applications

Sachant la valeur d’une fonction f en N points M1 , ..., MN d’enveloppe convexe E
(l’enveloppe convexe est l’ensemble des barycentres de ces points à coefficients positifs),
on peut vouloir :
– Interpoler : estimer f en un autre point à l’intérieur de E sans évaluer f
– Extrapoler : estimer f en un autre point à l’extérieur de E sans évaluer f
On souhaite donc déterminer fˆ(M ), M ∈ E sans avoir à faire une nouvelle évaluation
de f en M.
Savoir interpoler, c’est savoir dériver et intégrer à une dimension.
ˆ ) = D(fˆ) et I(f
ˆ ) = I(fˆ).
Le principe général de ces deux opérations est d’écrire : D(f
Les méthodes dites classiques sont basées sur l’interpolation polynômiale.
Exemple 5.1
Applications de l’interpolation :
– Approximation de fonctions coûteuses à évaluer : prospection pétrolière (le coût
d’un forage est très élevé), optimisation d’un moteur...
Applications de l’extrapolation :
– Prédictions météorologiques

26

5.2. PROBLÈMES EN UNE DIMENSION - INTERPOLATION POLYNÔMIALE

5.2
5.2.1

27

Problèmes en une dimension - Interpolation polynômiale
Polynômes de Lagrange

L’interpolation par un polynôme de Lagrange a pour objectif de trouver l’unique polynôme de degré N-1 valant f1 , ..., fN aux points distincts x1 , ..., xN , que l’on note fˆ.
Pour cela, on utilise la base suivante de l’espace vectoriel polynomial N −1 [X] : la
famille des (Li )i∈{1,...,N } , avec

R

n
Y

x − xj
x − xj
j=1j6=i i

Li (x) =

Li est le polynôme valant 1 au point xi et 0 aux autres points.
On a alors l’expression de fˆ suivante :
fˆ =

n
X

fi Li (x)

i=0

Il apparait donc de façon évidente des procédés de dérivation et d’intégration très
simples, car :

fˆ(n) (x) =

N
X

(n)

fi Li (x)

i=1

Z b

fˆ(x)dx =

a

N
X
i=1

Z b

fi

Li (x)dx
a

Cependant, les Li ne sont jamais explicitement utilisés en pratique car trop coûteux
(en O(N 2 ) opérations).
Si les points sont régulièrement espacés, il existe des formules toutes faites que nous
allons voir dans la suite.
Pour déterminer les coefficients du polynôme interpolateur, il faut résoudre un système
de Van der Monde. Pour déterminer la valeur du polynôme interpolateur en un point, on
pourra utiliser la méthode de Neville.

5.2.2

Système de Van der Monde

On cherche à résoudre le système ∀i ∈ [|1, N |],
a1 + a2 x + ... + aN xN −1 .
On a donc :

fˆ(xi ) = f (xi ). On note fˆ(x) =

28CHAPITRE 5. INTERPOLATION ET EXTRAPOLATION - DÉRIVATION ET INTÉGRATION NU










x1
x2

···
···

−1
xN
1
−1  
xN

2


1 xN

···

−1
xN
N

1
1
..
.

a1
..
.

 .
  ..
aN





 
 
 
=
 
 

f1
..
.
..
.
fN









Les matrices de Van der Monde ne sont pas bien conditionnées, même pour des N
faibles. On ne les utilisera pas souvent. La résolution de ce système est en O(N 3 ) opérations, contre O(N 2 ) pour l’utilisation des formules de Lagrange.

5.2.3

Formule barycentrique de Lagrange

wi
On réécrit Li (x) = L(x) x−x
, avec L(x) =
i

QN

j=1 (x

− xj ) et wi = QN

1

j=1,j6=i

(xi −xj )

. On a

donc :
fˆ(x) = L(x)

N
X
i=1

fi

wi
x − xi

Pour des abscisses données, la complexité devient O(N) opérations. Rajouter une nouvelle abscisse se fait en O(N) : on doit remplacer L(x) par (x − xN +1 )L(x) et remplacer
wi
wi par xi −x
.
N +1
On remarque que :
1 = L(x)

N
X
i=1

5.2.4

wi
x − xi

Limites de l’interpolation polynômiale

Le phénomène de Runge apparait assez rapidement : en effet, il arrive que plus on
accroit le nombre de points de l’interpolation, plus l’interpolant oscille. Cela peut être
limité en "resserrant" les abscisses. De plus, l’interpolation est sensible aux erreurs : ce
problème est intrinsèquement mal conditionné.

5.2.5

Interpolation par les splines cubiques

On cherche à interpoler par des fonctions cubiques par morceaux sur des intervalles
Ik = [xk , xk+1 ].
On a alors, sur l’intervalle Ik :
fˆ|Ik (x) = ak + bk x + ck x2 + dk x3
On a donc 4 inconnues par intervalle. De plus, on impose la continuité de f d’ordre C0 ,
C1 et C2 . On impose deux relations supplémentaires sur les bords.
On constate donc que fˆ00 est continue et affine par morceaux.

5.3. INTERPOLATION À D DIMENSIONS

29

x−x
k −x
+ uk+1 xkx−x
où les uk sont les valeurs de fˆ00 (pour
Ainsi, fˆ00 |Ik (x) = uk xk −xk+1
k+1
k+1
l’instant inconnues).
En intégrant deux fois l’expression de fˆ00 , on obtient :

1
uk
1 uk+1
fˆ|Ik (x) =
(x − xk+1 )3 −
(xk − x)3 + αk (x − xk+1 ) + βk (xk − x)
6 xk − xk+1
6 x − xk+1
Les deux constantes d’intégration s’obtiennent en écrivant les deux relations suivantes :
fˆ|Ik (xk ) = fk
fˆ|I (xk+1 ) = fk+1
k

Les dérivées secondes sont solutions d’un système tridiagonal qui se résout en O(N)
opérations.

5.2.6

Méthodes à noyaux

On cherche une approximation de la forme :
fˆ(x) =

X

wi Kµ (xi , x)

i

telle que fˆ(xj ) = i wi Kµ (xi , xj ) = f (xj ).
Les poids sont solution d’un système linéaire (matrice de Gram) : Kµ w = f .
On parle de méthode à noyau si Kµ est symétrique définie positive ∀(x1 , ..., xN ).
On peut interpréter Kµ (x, y) comme une mesure de similarité ou d’influence entre x
et y.
La résolution du système se fait en O(N 3 ) opérations.
P

Exemple 5.2
Les Radial Basis Functions (RBF) sont souvent utilisées : le noyau s’exprime en fonction de la distance uniquement entre x et y. On pourra utiliser par exemple la gaus2
sienne e−(µr) .
De nombreuses méthodes d’interpolation classiques (polynômes et splines) peuvent
être vues comme des cas particuliers ou des limites de méthodes à noyaux.

5.3

Interpolation à D dimensions

Si les points sont disposés sur une grille (par exemple pour le traitement d’image), les
méthodes envisagées précédemment sont toujours valables.
Sinon, il faudra réaliser un maillage pour pouvoir utiliser l’interpolation par les splines,
ce qui est coûteux en temps de calcul. Sinon, on pourra utiliser la méthode des noyaux.

30CHAPITRE 5. INTERPOLATION ET EXTRAPOLATION - DÉRIVATION ET INTÉGRATION NU

5.4
5.4.1

Dérivation numérique
Formules des différences finies

On suppose une fonction f connue en x1 , ..., xn régulièrement espacés d’un pas h.
Formules du premier ordre
Définition 5.1
On appelle différence finie arrière d’ordre 1 la quantité :
f (xi ) − f (xi − h)
f (xi ) − f (xi−1 )
=
= f 0 (xi ) + O(h)
xi − xi−1
h
O(h) est appelée l’erreur de discrétisation.
Définition 5.2
On appelle différence finie avant d’ordre 1 la quantité :
f (xi+1 ) − f (xi )
f (xi + h) − f (xi )
=
= f 0 (xi ) + O(h)
xi+1 − xi
h
O(h) est appelée l’erreur de discrétisation.

Formules du second ordre
Définition 5.3
On appelle différence finie centrée d’ordre 2 la quantité :
f (xi + h) − f (xi − h)
= f 0 (xi ) + O(h2 )
2h
Définition 5.4
On appelle différence finie arrière d’ordre 2 la quantité :
3f (xi ) − 4f (xi − h) + f (xi − 2h)
= f 0 (xi ) + O(h2 )
2h
La différence finie d’ordre N quelconque nécessite la connaissance de plus de points.
Propriété 5.1
Le calcul de la dérivée seconde peut se faire par la dérivée seconde centrée d’ordre 2 :
f (xi + h) − 2f (xi ) + f (xi − h)
= f 00 (xi ) + O(h2 )
h2

5.5. INTÉGRATION NUMÉRIQUE

31

Propriété 5.2
Le calcul de la dérivée seconde peut se faire par la dérivée seconde arrière d’ordre 2 :
2f (xi ) − 5f (xi − h) + 4f (xi − 2h) − f (xi − 3h)
= f 00 (xi ) + O(h2 )
h2

5.4.2

Choix d’un pas optimal

Il ne sert à rien de diminuer le pas h "à l’infini". L’erreur d’arrondi finit par submerger
l’erreur de discrétisation. On peut toujours trouver un compromis entre erreur de discrétisation et erreur d’arrondi.
Si possible, on utilisera l’extrapolation de Richardson pour raffiner la valeur de la
dérivée au point désiré.
Par exemple, si l’échantillonnage est uniforme, on calculera f (xi )−f4h(xi −4h) , f (xi )−f2h(xi −2h) ,
f (xi )−f (xi −h)
puis on passe à la limite.
h

5.5
5.5.1

Intégration numérique
Intégration à une dimension

On cherche à calculer une intégrale du type ab w(x)f (x)dx. Comment choisit-on w(x)
et comment gradue-t-on l’axe des abscisses ? La méthode de Gauss permet de répondre
au problème.
Lorsque l’intégrande varie fortement, et que la subdivision de l’intervalle est à choisir,
on utilisera les méthodes adaptatives.
Dans le cas d’abscisses régulièrement espacées, on utilise les méthodes de Newton-Cotes
composites ou la méthode de Romberg.
R

Méthode de Gauss en théorie
Pour tout w(x) > 0 tel que ab w(x)xk dx existe, on montre qu’il existe des abscisses
x1 , ..., xN et des poids p1 , ..., pN tels que :
R

Z b

w(x)f (x)dx =
a

N
X

pi fi

i=1

soit exact pour tout polynôme de degré ≤ 2N − 1.
Les xi sont les racines d’une famille de polynômes orthogonaux construite à partir de
w(x).
Exemple 5.3
– Polynômes de Legendre : w(x) = 1 sur [−1; 1]

32CHAPITRE 5. INTERPOLATION ET EXTRAPOLATION - DÉRIVATION ET INTÉGRATION NU

– Polynômes de Chebyshev : w(x) = 1 − x2 sur [−1; 1]
– Polynômes de Laguerre : w(x) = e−x sur [0; ∞[
2
– Polynômes d’Hermite : w(x) = e−x sur

R

Méthode de Gauss en pratique
Le cas le plus courant est le cas de Gauss-Legendre (w(x) = 1 sur [−1; 1]). Si l’intervalle
ne convient pas, on effectue un changement de variable pour se ramener à ce cas.
L’erreur commise vaut :
(b − a)2N +1 (N !)4 (2N )
f
(ξ)
(2N + 1)[(2N )!]3
Plus l’ordre est important, plus les xi s’accumulent sur les bords de l’intervalle : on a
atténuation du phénomène de Runge. Ce phénomène peut se résumer comme suit : l’augmentation du nombre N des points d’interpolation ne constitue pas une bonne stratégie
d’approximation.
Méthodes adaptatives
Il arrive parfois qu’une subdivision systématique de tous les sous-intervalles ne soit pas
optimale (phénomène de Runge).
L’idée des méthodes adaptatives consiste à ne raffiner que localement, en fonction
d’une erreur "a posteriori".
Par exemple, on utilise une quadrature d’ordre faible sur [a; b], ainsi que sur [a; (a+b)/2]
et [(a+b)/2; b]. Si le premier résultat est comparable à la somme des deux autres, on stoppe.
Sinon, on réitère en remplaçant [a; b] par [a; (a + b)/2] et [(a + b)/2; b].
Méthodes de Newton-Cotes
La méthode s’obtient via le polynôme interpolateur de Lagrange lorsque [a ;b] est
subdivisé en N-1 intervalles de longueur égale.
La méthode de Newton-Cotes est exacte pour les polynômes de degré N-1.
Exemple 5.4
Méthode des trapèzes (degré 1)
b−a
(f (a) + f (b)) '
2

Z b
a

f (x)dx −

(b − a)3 00
f (ξ)
12

Les méthodes de degré supérieur sont la méthode de Simpson (3 points, degré 2), la
méthode de Milne (4 points), la méthode de Boole (5 points) etc...
Cette méthode n’est jamais utilisée aux ordres elevés à cause du phénomène de Runge,
et est inutile aux ordres faibles.

5.5. INTÉGRATION NUMÉRIQUE

33

Méthodes de Newton-Cotes composites
On suppose [a ;b] subdivisé en P intervalles égaux de taille h, sur lesquels on applique
une méthode de Newton-Cotes d’ordre faible.
Exemple 5.5
Méthode des trapèzes composite
h
(f (a) + 4f (a + h) + ... + f (b)) '
2

Z b

f (x)dx −

a

b − a 2 00
h f (ξ)
12

Méthode de Romberg
La méthode de Romberg est une méthode récursive.
Pour la méthode des trapèzes composite, on a :
ˆ
I(h)
=

Z b

f (x)dx +
a

i
X

nf tyCk h2k

k=1

La méthode de Romberg consiste à appliquer la méthode de Richardson à notre méthode des trapèzes :
b
4ˆ h

I( ) − I(h)
'
f (x)dx + O(h4 )
3 2
3
a
A chaque itération de Richardson, on gagne 2 degrés sur l’erreur de discrétisation. Au
niveau 1, on a la méthode de Simpson composite. Au niveau N, on est plus précis que la
méthode de Newton-Cotes correspondante.

Z

5.5.2

Intégration en D dimensions

Fléau de la dimension
Une quadrature ayant une erreur en hn en 1 dimension a une erreur en hn/D en D
dimensions. Pour une précision donnée, P évaluations de f en une dimension deviennent P D
évaluations de f en D dimensions. Ainsi, accroitre la dimension, c’est accroitre l’imprécision
et le coût opératoire.
En D dimensions, il n’existe plus de notion d’intervalle. Les problèmes de la vraie vie ne
sont pas définis analytiquement. Si le domaine peut être approximé aisément (maillage ou
maillage adaptatif), des quadratures composites de Newton-Cotes sont possibles. Sinon,
on utilisera la méthode de Monte-Carlo exposée dans le paragraphe suivant.
Méthode de Monte-Carlo
Théorème 5.1
Loi des grands nombres
Si x est une variable aléatoire uniforme sur un domaine Ω, alors

34CHAPITRE 5. INTERPOLATION ET EXTRAPOLATION - DÉRIVATION ET INTÉGRATION NU

Z
N
1 X
1
lim
g(xi ) =
g(x)dx
N →+∞ N
Ω Ω
i=1
R

On souhaite calculer l’intégrale V f (x)dx.
La méthode de Monte-Carlo consiste à choisir :
– Ω un hypervolume contenant V, facile à calculer (un hypercube par exemple)
– x uniforme sur Ω
– g(x)=f(x) si x ∈ V
– g(x)=0 sinon
R
σf
Ω PN
On a N
i=1 g(xi ) = V f (x)dx ± Ω N 1/2 où σf dépend de la "lisseté" de la fonction à
intégrer et où N ne dépend pas de la dimension du problème.
De nombreux raffinements peuvent être apportés à cette méthode (pénalisation, utilisation de lois de probabilité non uniformes, etc...).

Chapitre 6

Problèmes aux valeurs initiales et
aux limites
6.1
6.1.1

Problèmes aux valeurs initiales
Contexte de l’étude

L’objectif est de s’intéresser à des problèmes régis par des équations du type :


 dy = f (y, t)
dt



y(0) = y0

Connaissant f, on cherche à déterminer une estimation de y(t) pour t>0 avec une bonne
précision et en un minimum d’opérations.

6.1.2

Notions élémentaires

L’objectif est d’approximer y(t) aux instants (t1 , ..., tn+1 ) par une suite :
yˆn+1 = yˆn + (tn+1 − tn )Φ(ˆ
yn+1 , yˆn , ..., yˆn−s , tn+1 , tn , ..., tn−s ; f )
Le choix d’une méthode de résolution correspond à un choix de l’application Φ. On
parlera de méthode explicite dans le cas particulier où yˆn+1 n’apparait pas dans le membre
de droite. Autrement, on parle de méthode implicite.
Dans le cas linéaire, la méthode peut être explicitée. Sinon, on utilise la méthode de
Newton ou bien une prédiction-correction.
Dans le cas où s=0, on parle de méthodes à un pas. A noter que les méthodes multi-pas
nécessitent généralement une initialisation à un pas.
On parle de méthode à pas constant si ∀i, ti+1 = ti + h.
A chaque étape, on redéfinit les conditions initiales ; on résout alors :
35

36

CHAPITRE 6. PROBLÈMES AUX VALEURS INITIALES ET AUX LIMITES



 dz = f (z, t)
dt



z(tn ) = yˆn

Définition 6.1
On appelle erreur locale de discrétisation la quantité :
n+1 = yˆn+1 − z(tn+1 )
Définition 6.2
On appelle erreur globale de discrétisation la quantité :
en+1 = yˆn+1 − y(tn+1 )
Définition 6.3
Une méthode est dite convergente si l’erreur globale tend vers 0 lorsque les pas tendent
vers 0.
Définition 6.4
Une méthode est dite d’ordre p si l’erreur locale est en O(hp+1 )
Définition 6.5
Une méthode est dite stable si pour tout pas h donné la suite yˆn est bornée.
Exemple 6.1
Méthodes d’Euler
On cherche à approximer la dérivée en t = tn . On a :
dy
yˆn+1 − yˆn
=
= f (ˆ
yn , tn )
dt
tn+1 − tn
On obtient donc la relation :
yˆn+1 = yˆn + (tn+1 − tn )f (ˆ
yn , tn )
La méthode d’Euler en dérivant vers l’avant est donc explicite. Si l’on avait utilisé
une dérivation vers l’arrière, la méthode serait implicite car on aurait :
yˆn+1 = yˆn + (tn+1 − tn )f (ˆ
yn+1 , tn+1 )

6.1. PROBLÈMES AUX VALEURS INITIALES

6.1.3

37

Méthodes de résolution

Méthodes de Runge-Kutta
L’objectif de ces méthodes est d’approximer f(y,t) entre tn et tn+1 via une quadrature.
On remplace les valeurs inconnues de y par des approximations.
Exemple 6.2
Méthode du point milieu de Runge (en estimant yˆn+ 1 par la méthode d’Euler) :
2

yˆn+1 = yˆn + hf (ˆ
yn+ 1 , tn +
2

yˆn+ 1 = yˆn +
2

h
)
2

h
f (ˆ
yn , tn )
2

La méthode de Runge est d’ordre 2 (mais ce n’est pas la seule méthode de RungeKutta d’ordre 2).

Exemples de méthodes multi-pas
Exemple 6.3
Méthode d’Adams-Bashforth d’ordre 2
Z tn+1
dy
tn

dt

Z tn+1

dt =

f (y, t)dt '

tn

h
(3f (ˆ
yn , tn ) − f (ˆ
yn−1 , tn−1 ))
2

Exemple 6.4
Méthode d’Adams-Moulton d’ordre 2 (trapèzes)
Z tn+1
dy
tn

dt

Z tn+1

dt =
tn

f (y, t)dt '

h
(f (ˆ
yn+1 , tn+1 ) + f (ˆ
yn , tn ))
2

Exemple 6.5
Méthode de Gear d’ordre 2
dy
1
'
(3ˆ
yn+1 − 4ˆ
yn + yˆn−1 )
dt |tn+1
2h
Le domaine de stabilité de ces méthodes diminue lorsque l’ordre augmente. Les méthodes d’Adams-Moulton et de Gear d’ordre inférieur à 3 sont A-stables. Il n’existe aucune
méthode multi-pas A-stable d’ordre strictement supérieur à 2.
Ces méthodes ne nécessitent qu’une nouvelle évaluation de f par pas de calcul. Elles
sont couramment utilisées dans des schémas prédiction-correction à pas adaptatif.

38

CHAPITRE 6. PROBLÈMES AUX VALEURS INITIALES ET AUX LIMITES

Méthodes adaptatives
Les méthodes à pas constant sont peu utilisées en pratique. Il est délicat d’utiliser un
schéma explicite sans crainte pour la stabilité.
Les méthodes adaptatives utilisent la différence entre une étape de prédiction et une
étape de correction pour déterminer si la taille courante du pas convient.
Exemple 6.6
Méthode Adams-Bashforth-Moulton d’ordre 2.
P
yˆn+1
' yˆn +

h
(3f (ˆ
yn , tn ) − f (yn−1 , tn−1 ))
2

C
yˆn+1
' yˆn +

h
P
(3f (ˆ
yn+1
, tn+1 ) − f (yn , tn ))
2

On montre que :
P
y(tn+1 ) − yˆn+1
= GP y (3) (τn+1 )h3

(6.1)

C
0
y(tn+1 ) − yˆn+1
= GC y (3) (τn+1
)h3

(6.2)

D’où :
C
P
0
(6.1) − (6.2) −→ yˆn+1
− yˆn+1
' (GC − GP )y (3) (τn+1
)h3

C
C
P
(2) −→ n+1 = y(tn+1 ) − yˆn+1
' (ˆ
yn+1
− yˆn+1
)

GC
= ˆn+1
GC − G P

Si |ˆ
n+1 | > emax , on divise le pas par deux jusqu’à satisfaction (ou plantage). Pour une
méthode d’ordre p, on estime un nouveau pas optimal à partir de :
ˆ = hn+1
h



emax

n+1 |



1
p+1

On impose évidemment une limite raisonnable à la valeur du nouveau pas.

6.2
6.2.1

Problèmes aux limites
Notations
2

∂ u
On notera ∂t u la dérivée partielle ∂u
∂t , et ∂xy u la dérivée partielle ∂x∂y .
0
On notera u˙ la dérivée droite temporelle du
dt , et u la dérivée droite spatiale
On notera Γ = ∂Ω le contour (ou la frontière) d’un domaine d’étude Ω.

Définition 6.6
Le produit scalaire de deux fonctions u et v sur l’espace Ω est défini par :

du
dx .

6.2. PROBLÈMES AUX LIMITES

39

Z

< u|v >=

u(x)v(x)dΩ


6.2.2

Cadre de l’étude

On suppose que l’on cherche à résoudre le problème suivant : il s’agit de trouver une
fonction u(t,x,y,z) définie sur [0; +∞[×Ω qui vérifie :
– Une équation aux dérivées partielles (EDP) :
Lt,x,y,z (u, ∂t u, ∂x u, ∂y u, ∂z u, ∂tt u, ∂tx u, ...) + f (t, x, y, z) = 0
(x, y, z) ∈ Ω, t > 0
– Des conditions aux limites (CL) :
Bt,x,y,z (u, ∂t u, ∂x u, ∂y u, ∂z u, ∂tt u, ∂tx u, ...) + q(t, x, y, z) = 0
(x, y, z) ∈ Γ = ∂Ω, t > 0
– Des conditions initiales :
u(0, x, y, z) = u0
∂t u(0, x, y, z) = v0
Définition 6.7
On dit que le problème est linéaire si et seulement si Lt,x,y,z et Bt,x,y,z sont linéaires
en u.
Définition 6.8
On dit que le problème est stationnaire si et seulement si il n’apparaît pas de dérivée
par rapport au temps dans la définition du problème.
Définition 6.9
On dit qu’il y a homogénéité :
– de l’EDP si f est nulle
– des CL si q est nulle
Définition 6.10
On dit que l’équation est scalaire si et seulement si u est à valeurs dans
on dit que l’équation est vectorielle

R. Autrement,

Définition 6.11
On dit que le problème est bien posé si et seulement si la solution existe, est unique,
et dépend continûment des sources f et q.

40

CHAPITRE 6. PROBLÈMES AUX VALEURS INITIALES ET AUX LIMITES

Exemple 6.7
Equations stationnaires linéaires
– Equation de Poisson (homogène, 2ème ordre en espace)
– Equation de Laplace (non-homogène, 2ème ordre en espace)
Exemple 6.8
Equations transitoires linéaires
– Equation de la chaleur / Black-Scholes (1er ordre en temps, 2ème ordre en espace)
– Equation des cordes / membranes / ondes (2ème ordre en temps, 2ème ordre en
espace)
– Equation des poutres d’Euler-Bernoulli (2ème ordre en temps, 4ème ordre en
espace)
Exemple 6.9
Equations transitoires non-linéaires
– Equation de Navier-Stokes (2ème ordre en temps, 2ème ordre en espace)
Exemple 6.10
CL pour EDP du 2ème ordre en espace
– Conditions aux limites de Dirichlet :
u = hD (x, y, z)

(x, y, z) ∈ ΓD

– Conditions aux limites de Neuman :
∂u
~ n = hN (x, y, z)
= ∇u.~
∂~n
– Conditions aux limites de Robin :
∂u
= hR (x, y, z)
∂~n
réunion disjointe.

u+α
Avec Γ = ΓD ∪ ΓR ∪ ΓN

6.2.3

(x, y, z) ∈ ΓN

(x, y, z) ∈ ΓR

Démarche générale

Cas des EDP stationnaires
On commence toujours par discrétiser le problème : on transforme u(x,y,z) en u
ˆ.
Dans le cas d’un système linéaire K u
ˆ + f = 0, on utilisera les méthodes classiques
associées (LU, QR, Gauss-Seidel...).
Dans le cas d’un système non linéaire R(ˆ
u) = 0, on utilisera les méthodes classiques
associées (NR, Quasi-Newton...).

6.2. PROBLÈMES AUX LIMITES

41

Cas des EDP transitoires
On commence par une discrétisation spatiale, puis une discrétisation temporelle : on
ˆ
transforme u(t,x,y,z) en u
ˆ(t), puis en u
ˆ(tn+1 ).
Ensuite, nous serons amenés à résoudre un système d’EDO linéaire de la forme K u
ˆ+
˙
¨
˙
¨
Cu
ˆ + Mu
ˆ + f = 0, ou bien non linéaire, de la forme R u
ˆ, u
ˆ, u
ˆ

6.2.4

Problèmes stationnaires

La méthode de Fourier : une méthode historique
En 1822, la théorie analytique de la chaleur est mise en place. L’idée de la méthode
de résolution de Fourier est la suivante : on part du constat que, sur un segment donné,
presque toutes les fonctions peuvent s’écrire comme une somme infinie de fonctions trigonométriques.
On souhaite alors approximer la solution d’une EDP du second ordre avec des CL
homogènes
u00 + f = 0
sous la forme d’une série tronquée
u
ˆ(x) =

P
X

up Np (x)

p=1

avec lim u
ˆ = u où les Np sont des fonctions trigonométriques vérifiant également les
P →∞

CL : Np (x) = sin(λp x) = sin( 2p−1
2 πx).
Comme une fonction continue est nulle si et seulement si tous ses coefficients sont
P
nuls, on obtient la résolution de R(ˆ
u) = − up λ2p Np + f = 0 en identifiant coefficient par
coefficient.
L’obtention d’un coefficient de Fourier up se fait par projection sur la base des Np :
up =< f |Np >=

Z 1

Np (x)f (x)dx
0

On aboutit alors à la résolution triviale d’un système linéaire diagonal :
1
− uk λ2k + fk = 0
2
Il y a cependant nécessité de calculer des intégrales.
L’extension au cas de CL non homogènes est triviale :
∀k,

u
ˆ(x) = u
¯(x) +

p
X
p=1

où u
¯(x) vérifie les CL non homogènes.

up Np (x)

42

CHAPITRE 6. PROBLÈMES AUX VALEURS INITIALES ET AUX LIMITES

Généralisation de la méthode de Fourier
Pourquoi s’astreindre à utiliser les fonctions trigonométriques comme base de l’espace
des fonctions ? D’autres bases existent : polynômes, fonctions polynômiales par morceaux,
ondelettes...
D’autre part, pourquoi s’astreindre à projeter l’EDP approchée sur cette même base
de fonctions ? Par exemple, on peut essayer de vérifier l’EDP :
– exactement en P points bien choisis : méthode de collocation
– en moyenne sur P sous-domaines bien choisis : méthode de sous-domaine
La méthode des résidus pondérés sera à la base de 99,9% des méthodes de résolution
d’EDP existantes.
Méthodes des résidus pondérés
Pour une équation avec des CL homogènes : R(u) = L(u) + f = 0, il s’agit de décomposer la solution sur une base de fonctions bien choisies :
u
ˆ(x, y, z) =

P
X

up Np (x, y, z)

p=1

R(ˆ
u) ne peut pas etre nul partout : on impose la nullité de la projection du résidu R(ˆ
u)
sur une autre base de fonctions de test φq bien choisie :
∀q ∈ {1, ..., P }

< R(ˆ
u)|φq >= 0

.
On arrive donc à un système linéaire de taille P × P si L est linéaire :
Ku
ˆ+f =0
Cette méthode s’étend aisément aux cas des CL non homogènes, ou bien au cas non
linéaire. On distingue les méthodes des types suivants :
– Galerkin (fonctions de base = fonctions test) : on obtient une matrice symétrique
(souvent définie positive), voire diagonale par les méthodes spectrales. Le calcul
d’intégrales est requis.
– Collocation (fonctions test = Diracs) : La matrice est quelconque, et aucun calcul
complexe hors résolution est requis.
– Sous-domaine (fonctions test = fonctions caractéristiques d’un sous-domaine) : La
matrice est quelconque. Le calcul d’intégrales est requis.
On peut utiliser des fonctions de base globales comme c’est le cas dans la méthode de
Fourier, ou bien des fonctions locales (méthodes des éléments finis et différences finies).
Méthode des différences finies - approche intuitive
Le principe de la méthode est le suivant : on commence par faire un quadrillage du
domaine. Ensuite, on remplace les dérivées par rapport à l’espace par des différences finies.

6.2. PROBLÈMES AUX LIMITES

43

Cette méthode aboutit à des matrices creuses (tridiagonales si le problème est à une
dimension). Elle peut être obtenue comme une méthode des résidus pondérés de type
collocation.
Exemple 6.11
2

On s’intéresse à résoudre le problème suivant : ddxu2 + f = 0 sur Ω = [0, 1]. On note
x0 , ..., xN une subdivision régulière de Ω de pas h. On a donc N+1 points et N intervalles
de taille h.
On a donc :
– N-1 équations données par l’équation différentielle aux N-1 points xq (q = 1, ..., N −
u(xq+1 )−2u(xq )+u(xq−1 )
1) :
+ f (xq ) = 0
h2
– Une équation donnée par la condition initiale u(x0 ) = 1
u(xN )−u(xN −1 )
– Une équation par la condition de la dérivée
=0
h
On a donc bien N+1 équations, et N+1 inconnues (les valeurs de la fonctions aux
points considérés.
Le système à résoudre est donc le suivant :











− h22
1
h2

0

1
h2

0

..

.

..

.

..

.

..

.

..

.

..

.

..

.

1
h2

1
h2
− h22





 
 
 
 
×
 
 
 
 

u(x1 )
..
.
..
.
..
.





 
 
 
 
+
 
 
 
 

u(xN −1 )

f (x1 )
..
.
..
.
..
.
f (xN −1 )




 
 
 
 
+
 
 
 


1
h2



0
..
.
..
.
0





=0




On constate que la matrice du système est très creuse.
Un problème de cette méthode est qu’elle ne s’applique qu’à un domaine quadrillable.
De plus, si P est élevé, on a de nombreux dangers qui apparaissent (espace mémoire, erreur
d’arrondi, conditionnement, etc...)
Méthode des éléments finis
A une dimension, cette méthode est identique à la précédente.
On commence par faire un maillage du domaine, puis on passe d’une formulation forte
du problème à une formulation faible :
– On effectue une projection et une intégration par parties (formule de Green)
– On choisit des fonctions de base/test polynômiales (linéaires si possible) par morceaux)
– u
ˆ est imposé à priori sur ΓD , les fonctions tests sont nulles sur ΓD .
Par exemple, on prendra :



 Np (x) linéaire par morceaux

N (x ) = 1

p p


 N (x ) = 0 si q 6= p
p q

44

CHAPITRE 6. PROBLÈMES AUX VALEURS INITIALES ET AUX LIMITES
P

On a alors u
ˆ = up Np (x) où up est la valeur de u
ˆ au noeud xp .
Cette méthode aboutit à une matrice creuse, plus simple et plus souple qu’avec la
méthode des différences finies : le domaine est plus complexe, le raffinement local est plus
simple.
Exemple 6.12
Formulation forte :
∆u + f = 0 sur Ω ⇐⇒ ∀v,

Z

(∆u + f )vdΩ = 0


Formulation faible :
∀v,

Z

f vdΩ −



Z

~
~
gradu.
gradvdΩ
+



Z

~u.~nvdΓ = 0
Γ

On décompose alors u sur une base de fonctions linéaires par morceaux, et v est
choisi dans une base de fonctions tests φq linéaires mar morceaux.

6.2.5

Problèmes transitoires

Les problèmes transitoires se résolvent de la même façon que les problèmes stationnaires, à la différence près que les coefficients dépendent du temps :
u
ˆ(t, x, y, z) =

P
X

up (t)Np (x, y, z)

p=1

La résolution aboutit à des systèmes d’EDO.
On utilise des schémas classiques (Runge-Kutta, Adams, Gear) dans les deux cas suivants :
– Equations d’ordre 1 en temps
– Equations d’ordre supérieur à 1 mais P petit
Sinon, lorsque l’ordre temporel est supérieur à 1 et P est grand, on pourra utiliser les
méthodes des différences finies en temps, ou le schéma de Newmark.




Télécharger le fichier (PDF)

Méthodes Numériques et Optimisation.pdf (PDF, 406 Ko)

Télécharger
Formats alternatifs: ZIP







Documents similaires


professeur benzine rachid optimisation sans contraintes tome1
professeur benzine rachid cours optimisation sans contraintes tome1
14 problemes ouverts sur le gradient conjugue sujets de theses de doctorat
gradient conjugue cas quadratique et non quadratique
methodes quasi newton
convergence et vitesse convergence methode du gradient

Sur le même sujet..