Analyse Numérique Appliquée

Julien Réveillon
François-xavier Demoulin
demoulin@[Link]

1
 Plan du cours
• Recherche des zéros
• Dérivée d’une fonction.
• Intégration numérique.
• Résolution d’équations temporelles.
• Résolution de systèmes d’équ. linéaires,
inversion de matrices.
• Interpolation et lissage de courbes

2
 0 d ’une fonction

Objectif : trouver les solutions d ’une fonction f ( x)  0

3
 0 d ’une fonction

Méthode de la bissection :

Méthodologie (pour 1 racine):

1. La racine recherchée est tout d’abord localisée
grossièrement dans un intervalle [x0,x1] grâce au fait
que le produit f(x0)*f(x1) est négatif.
2. La position x2 du centre de l’intervalle [x0,x1] est
alors déterminé x2 = (x0+x1)/2.
3. Si f(x0)*f(x2) est négatif alors x1 prend la valeur de
x2 sinon, x0 prend la valeur de x2.
4. Les étapes 2 et 3 sont effectuées à nouveau jusqu ’à
ce que la différence |x0-x1| soit inférieure à l’erreur
maximale que l ’on s’était fixée.

4
 0 d ’une fonction
Méthode de Newton :
Méthodologie :
Prenons le développement limité à l’ordre de 1 de la fonction considérée :

Or, notre désir le plus cher est d ’avoir f ( x1 )  0

f ( x0 )
Ce qui nous donne : x1  x0 
f ( x0 )
En suivant le même raisonnement, on obtient :
f ( x1 )
x2  x1 
f ( x1 )
f ( xn )
Soit, de manière générale : xn 1  xn 
f ( xn )
A nouveau, le calcul s’arrête lorsque xn 1  xn  
Afin d ’avoir une bonne chance de converger, il faut que x0 ne soit pas trop loin de la racine
5
 0 d ’une fonction
Méthode de la sécante :
Cette méthode permet d ’éviter le calcul de la dérivé de la fonction f(x) tout étant plus
efficace que la méthode de la bissection.

Le but du jeu est de déduire xn+1 en fonction de xn et xn-1 en appliquant le théorème de

Thales.

6
 0 d ’une fonction
Méthode de la sécante :

xn1  xn xn  xn 1

f ( xn 1 )  f ( xn ) f ( xn )  f ( xn 1 )
A nouveau, dans nos rêves les plus fous, nous souhaitons : f ( xn 1 )  0

f ( xn )
Soit : xn 1  xn  ( xn  xn 1 )
f ( xn )  f ( xn 1 )

A nouveau, le calcul s’arrête lorsque xn 1  xn  

7
 Calcul Dérivées

Deux cas de figure majeurs peuvent se présenter :

•devoir calculer la dérivée d ’une fonction analytique dont on connaît la valeur en
tout point: f (x )
•devoir calculer la dérivée d ’une fonction dont seules les valeurs en certains points
sont connues: f ( xi )
Dans le premier cas une résolution numérique est la bienvenue particulièrement si
l ’on ne sait pas (ou ne peut pas) effectuer une résolution algébrique.
Dans le second cas, la résolution numérique est nécessaire soit en utilisant
directement les valeurs discrètes de la fonction soit avec la mise en place du
mappage d’une fonction analytique sur les valeurs discrètes de la fonction que l ’on
souhaite dériver.

8
 Expansion en série de Taylor
Calcul de la dérivée 1ère de la fonction : f (x)
Expansion en séries de Taylor :
La fonction continue peut être développée en série de Taylor au voisinage de xi :
2 2 n n
 f  ( x  xi )   f  ( x  xi )   f
f ( x)  f ( xi )  ( x  xi )    2    n  H
 x i 2!  x i n!  x i
H représente les termes d’ordre élevé.
Il est possible d ’approcher la fonction en x  xi  1
2 2
 f  ( xi  1  xi )   f 
f ( xi  1)  f ( xi )  ( xi  1  xi )    2   H
 x i 2!  x i
Ou en x  xi  1
2 2
 f  ( xi  1  xi )   f 
f ( xi  1)  f ( xi )  ( xi  1  xi )    2   H
 x i 2!  x i 9
 Expansion en série de Taylor
Approximation de la dérivée en i (possibilité 1):
2 n 1 n
f
   f i  1  f i ( xi  1  xi )   f  ( xi  1  xi )   f
    2     n  H
 x i ( xi  1  xi ) 2!  x i n!  x  i
Approximation de la dérivée en i (possibilité 2):
2 n 1 n

 f f i  f i  1 ( xi  x i  1)   f  ( xi  1  xi )   f
     2    n  H
 x i ( xi  xi  1) 2!  x i n!  x i

Approximation de la dérivée en i (possibilité 3):

 f  fi  1  fi  1 ( xi  1  xi ) 2  ( xi  xi  1) 2   2 f 
     2  H
 x i xi  1  xi  1 2 xi  1  xi  1  x i

10
 Expansion en série de Taylor
Toutes les expressions vues précédemment sont exactes si tous les termes de la partie
droite sont retenus. Cependant des approximations sont nécessaires pour un calcul
pratique et il s’agit alors de tronquer les équations précédentes.

Après troncation, les dérivées s’écrivent :

 f   fi  1  fi Dérivée avant
 
 x i xi  1  xi Forward derivative
 f  fi  f i  1 Dérivée arrière
  
 x i xi  xi  1 Backward derivative

 f  fi  1  fi  1 Dérivée centrée
  
 x i xi  1  xi  1 Centered derivative

11
 Expansion en série de Taylor
Calcul de la dérivée 2nde de la fonction : f (x)
Expansion en séries de Taylor, maillage régulier :
La fonction continue peut être développée en série de Taylor au voisinage de xi :
2
 f  (x)  2 f  ( x ) n   n f 
f ( xi  l )  f ( xi )  x         n   H
 x i 2!  x i
2
n!  x i
H représente les termes d’ordre élevé. Avec: x  xi  xi 1
Il est possible d ’approcher la fonction en x  xi  1
2 2
 f  ( x)   f 
f ( xi  1)  f ( xi )  x    2   H
 x i 2!  x i
et en x  xi  2
2 2
 f  4 (  x )   f
f ( xi  2)  f ( xi )  2x    2   H
 x i 2!  x i 12
 Expansion en série de Taylor

Approximation de la dérivée seconde en i :

2(x) 2   2 f 
f ( xi  2)  2 f ( xi  1)   f ( xi )   2
2!  x i

 2 f  f ( xi  2)  2 f ( xi  1)  f ( xi )
 2  2
 x i (  x )

13
 Expansion en série de Taylor

Les dérivées secondes s’écrivent :

 2 f  fi  2  2 fi  1  fi Dérivée avant (ordre 1)

 2   2
 x i (  x ) Forward derivative
 2 f  f i  2 fi  1  fi  2 Dérivée arrière (ordre 1)
 2   2
 x i (  x ) Backward derivative
 2 f  fi  1  2 fi  fi  1
 2   2
Dérivée centrée (ordre 2)
 x i (  x ) Centered derivative

14
 Expansion en série de Taylor
Calcul de la dérivée croisée de la fonction : f ( x, y )

Expansion en séries de Taylor, maillage régulier : xi , yj

 f   f 
f ( x i 1 , y j 1 )  f ( x i , y j )   x     y  
 x  i  y  j
(x) 2  2 f  (y ) 2  2 f  2(  y x )   2 f 
        
2!   x 2
i 2!   y 2
j 2!  xy  j

15
 Expansion en série de Taylor

2 2
x y
f i 1, j 1  f i , j  xf x'  yf y'  xyf xy''  f x''  f y''
2 2
' ' '' x 2 '' y 2 ''
f i 1, j 1  f i , j  xf x  yf y  xyf xy  fx  fy
2 2
' ' x 2 '' y 2 ''
''
f i 1, j 1  f i , j  xf  yf  xyf 
x y fx  fy
xy
2 2
' ' x 2 '' y 2 ''
''
f i 1, j 1  f i , j  xf  yf  xyf 
x y fx  fy
xy
2 2

2 f f i 1, j 1  f i 1, j 1  f i 1, j 1  f i 1, j 1 Dérivée croisée


xy 4xy Cross derivative
16
 Évaluation Polynomiale
Cette méthode consiste à approcher la fonction discrète
fi à dériver par un polynôme p(x) d ’ordre n.
Par exemple, définissons p ( x)  Ax 2  Bx  C

Tel que : p ( xi )  f i , p( xi 1 )  f i 1 , p ( xi  2 )  f i  2
Il nous faut déterminer A,B,C, 3 points sont nécessaires.
Prenons xi comme point de référence.  xi  0
p ( xi )  f i  Axi2  Bxi  C  C
p ( xi 1 )  f i 1  Axi21  Bxi 1  C  A(x) 2  B (x)  C
p ( xi  2 )  f i  2  Axi2 2  Bxi  2  C  A(2x) 2  B (2x)  C
17
 Évaluation Polynomiale
3 équations et 3 inconnues (A,B,C).

 f i  2  4 f i 1  3 f i f i  2  2 f i 1  f i
C  fi B A
2x 2(x) 2

Les dérivées se déduisent alors de la dérivée du polynôme

p ( x)
 2 Ax  B  2 p( x)
2
 2A
x x
Soit, en xi  0
f   f i  2  4 f i 1  3 f i 2 f  f i  2  2 f i 1  f i
 B 2
  2A 
x i 2x x i (x) 2
18
 Évaluation Polynomiale
Cas d ’un maillage irrégulier.
Le maillage est connu. Écrivons pour le point xi ,   cte
xi  0
xi 1  x
xi  2  (1   )x
xi 3  (1     2 ) x
Comme précédemment :
p ( xi )  f i  Axi2  Bxi  C  C
p ( xi 1 )  f i 1  Axi21  Bxi 1  C  A(x) 2  B (x)  C
p ( xi  2 )  f i  2  Axi2 2  Bxi  2  C  A((1   )x) 2  B ((1   )x)  C
19
 Évaluation Polynomiale
Cas d ’un maillage irrégulier.
Et finalement
C  fi

 f i  2  (1   ) 2 f i 1   (  2) f i
B
 (1   )x

f i  2  (1   ) f i 1  f i
A
 (1   )(x) 2

Le taux d ’accroissement des mailles doit rester faible.

20
 Calcul Dérivées
Différences à droite - Ordre 1
fi f i 1 fi2 f i 3 fi4
f
x 1 1
x
2
2  f
x 1 2 1
x2
3
3  f
x 1 3 3 1
x3
4
4  f
x 1 4 6 4 1
x4

21
 Calcul Dérivées
Différences à gauche - Ordre 1
f i 4 f i 3 f i 2 f i 1 fi
f
x 1 1
x
2
2  f
x 1 2 1
x2
3
3  f
x 1 3 3 1
x3
4
4  f
x 1 4 6 4 1
x4

22
 Calcul Dérivées
Différences centrée - Ordre 2
fi2 f i 1 fi f i 1 fi2
f
2x 1 1
x
2
2  f
x 1 2 1
x2
3
3  f
2x 1 2 0 2 1
x3
4
4  f
x 1 4 6 4 1
x4

23
 Calcul Dérivées
Différences centrée - Ordre 4
f i 3 f i 2 f i 1 fi f i 1 fi2 f i 3
f
12x 1 8 0 8 1
x
2
2  f
12x 1 16  30 16 1
x2
3
3  f
8x 1 8 13 0  13 8 1
x3
4
4  f
6x 1 12  39 56  39 12 1
x4

24
 Intégration de fonction
Calcul de l ’intégrale de la fonction : f (x)

Méthode des trapèzes.

Le principe de la méthode consiste à diviser la fonction en une
succession de segments. L’intégrale de la fonction étant
représentée par l ’aire entre les segments successifs et l’axe des
abscisses.

25
 Intégration de fonction
Méthode des trapèzes.
La surface d’un trapèze s’écrit donc

 xi f i 1  f i
  f ( x) dx  xi
 xi 1 2

xi  xi  xi 1

La somme des aires de l ’ensemble des trapèzes donnera une

approximation de l ’aire de la fonction entre les bornes [a, b] séparées
en n intervalles. b
1 n
a f ( x)dx  2 
i 1
( f i 1  f i )xi

26
 Intégration de fonction

Pour simplifier les expressions, les n intervalles sont considérés

comme égaux x  x  b  a
i
n

La surface s’écrit donc:

b n 1
x x  
 f ( x)dx  ( f 0  f1 )   ( f n 1  f n )   f 0  f n  2 f i )
a
2 2  i 1 

27
 Intégration de fonction
Méthode de Simpson (pas Omer, l ’autre)
La méthode consiste à remplacer la fonction entre les points xi-1, xi+1
par un arc de parabole passant par fi-1, fi, fi+1.
Elle se démontre grâce à un développement de Taylor. Définissons
I(x) comme l ’intégrale partielle de f(x) entre a et x.
x
dI ( x)
I ( x)   f ( x' )dx' soit  f ( x)
a
dx
Il est possible d ’écrire :
I ( xi 1 )  I ( xi  x), I ( xi 1 )  I ( xi  x)

I ( xi 1 )  I ( x )  xf ( x ) 
x 
2
f ( x ) 
x 
3
f ( x ) 
x 
4
f (3) ( xi )   (5)
i i i i
2 3! 4!

I ( xi 1 )  I ( xi )  xf ( xi ) 
x 
2
f ( xi ) 
x 
3
f ( xi ) 
x  ( 3 )
4
 (5)
f ( xi ) 28
2 3! 4!
 Intégration de fonction
Méthode de Simpson
En soustrayant les deux dernières relations, on obtient l’aire Ai entre
les bornes xi-1, xi+1. Ai  I ( xi 1 )  I ( xi 1 )
x 3 
Ai  2xf ( xi )  f ( xi )   (5)
3
La dérivée seconde est remplacée par son expression utilisant les
différences centrées :
f ( xi 1 )  2 f ( xi )  f ( xi 1 )
f ( xi )  2
  ( 2)
x 
Et finalement :
Ai  2xf i 
x
 f i 1  2 f i  f i 1   x 
3
 (2)   (5)
3 3
x
Ai   f i 1  4 f i  f i 1    (5) 29
3
 Intégration de fonction
Méthode de Simpson
L ’aire totale peut être déterminée:
b

 f ( x)dx  A  A
a
1 3  ....  An1
Il faut donc que le nombre d’intervalles n soit pair.
b
x   f x   f n
 f ( x)dx  0  4 f1  f 2     n  2  4 f n 1  f n    5
a
3 3 2
n ba
 5   5   4 
2 2x
 
b
x    f (a)  f (b)  4 f ( x )  2 f ( x )   4
n 1 n2

 f ( x)dx 
3 

i 1
i 
i2
i

a
 imp pair 
30
Expression exacte pour l’intégration d ’un polynôme jusqu ’à ordre 3
 Intégration de fonction
Méthode de Simpson
Calculer l ’intégrale de e(x) entre 0 et 4.
4
I   e x dx
0

n4
1
I  e 0  e 4  4e1  e3   2e 2   53,862846
3
Erreur : 5.10 3

31
 Intégration de fonction
Formule généralisée de Newton-Cotes
xn

 f ( x)dx   0 f ( x0 )  1 f ( x1 )     n f ( xn )
x0
xn
x n
 f ( x)dx  
A k 0
ak f ( x0  kx)
x0

Nom n A a0 a1 a2 a3 a4 a5 a6
Trapèze 1 2 1 1
Simpson 2 3 1 4 1
Villarceau 4 45 14 64 24 64 14
Hardy 6 140 41 216 27 272 27 216 41
32
 Intégration de fonction
Intégration sur un pas quelconque - Méthode générale
Procédure en 2 étapes :
1) Approximation de f(x) par un polynôme
2) Substitution de l’intégrale par une combinaison des f(xi)
b

 f ( x)dx  a
0 f ( x0 )  a1 f ( x1 )    an f ( xn )
a

avec x0  a et xn  b

Les coefficients ak sont déterminés tels que l’on ait une égalité stricte
lorsque l’on remplace f(x) par un polynôme quelconque de degré
inférieur ou égal à n.

33
 Intégration de fonction
Intégration sur un pas quelconque - Méthode générale

b

 f ( x)  1   1dx  b  a  a0  a1   an
 a
 b
b2  a 2
 f ( x)  x   xdx   a0 x0  a1 x1   an xn
 a
2


b
 n n b n1  a n1 n n n
 f ( x )  x   x dx   a x
0 0  a x
1 1   a x
n n
 a
n  1

34
 Intégration de fonction
Intégration sur un pas quelconque - Méthode générale
Soit un système de n+1 équations à n+1 inconnues

1 1  1   a0   ba 
 x0 x  xn   a1   (b 2  a 2 ) 2 
1
     
         
 xn x n  xnn  an  (b n1  a n1 ) n  1
 1

Système linéaire à résoudre.

Si les points sont équidistant, on retrouve les coefficients de
Newton-Côtes.

35
 Intégration de fonction
Limites infinies

Il s ’agit de calculer : I ( x)   f ( x)dx
a

Il faut tout d ’abord s ’assurer que l ’intégrale

converge, ensuite, 2 méthodes sont envisageables :
1) Changement de variable permettant de
transformer des bornes infinies en bornes finies
(que l ’on sait résoudre).
2) Séparation de l ’intégrale en 2.
 b 

 f ( x)dx  f ( x)dx   f ( x)dx

a a b

36
 Intégration de fonction
Limites infinies
La première des 2 intégrales s’évalue de manière
classique. Pour la seconde, différentes conditions
peuvent intervenir.
•Si b est suffisamment grand, la seconde intégrale
peut être négligeable.
Pour savoir si b est suffisamment grand, on calcule
b 2b
I   f ( x)dx et I   f ( x)dx
a b

Il faut I  I

37
 Intégration de fonction
Limites infinies

•Si on connaît la forme asymptotique g(x) de f(x) alors, pour x

suffisamment grand, f(x)=g(x) et

 b 

 f ( x)dx   f ( x)dx   g ( x)dx

a a b

Souvent la solution analytique de cette deuxième intégrale est connue.

•Une dernière possibilité est de faire un changement de variable

sur la seconde intégrale 
 f ( x)dx
b
38
 Intégration de fonction
Singularités de l’intégrande
Les bornes sont finies mais la fonction n ’est pas
définie pour l ’une des deux bornes.
b

 f ( x)dx
a

Les bornes sont finies mais la fonction n’est pas

définie pour l’une des deux bornes.

1
1
dx  2 x 0  2
1
Exemple I  
0
x

39
 Equation d’évolution
Condition initiale
y0  y (t  0)

Equation différentielle ordinaire: EDO

dy (t )
 F ( y, t ) ordre 1
dt
ou
d 2 y (t ) dy
2
 F ( , y, t ) ordre 2
dt dt

40
 Equation d’évolution
Réduction de l’ordre de l’EDO:
On transforme l’EDO d’ordre élevée en un système d’EDO d’ordre 1
Exemple:
d2y
2
 ay  cos(bt )
dt
dy (t )
On pose: u (t ) 
dt
On obtient alors le système suivant:

 du (t )
 dt   ay (t )  cos(bt )

 dy (t )  u (t )
 dt 41
 Résolution méthode d’Euler
dy (t )
Soit à résoudre:  F ( y, t )
dt

Méthode d’Euler explicite

On connaît à ti=t0+iΔt: yi=y(t=ti) et on cherche yi+1
Dévellopement de Taylor:
dy
yi 1  yi  t  O ( t )
dt ti
Soit h  t , on passe de i à i  1suivant :
ti 1  ti  h

 yi 1  yi  hF ( yi , ti )  O (1)
42
 Méthode implicite
dy (t )
Soit à résoudre:  F ( y, t )
dt

Méthode d’Euler implicite

Développement de Taylor à partir de ti+1:
dy
yi  yi 1  t  O (t )
dt ti1

Soit h  t , on passe de i à i  1suivant :

ti 1  ti  h

 yi 1  yi  hF ( yi 1 , ti 1 )  O (1)
43
 Convergence

Un méthode est dit convergente si en diminuant le pas (ici le pas de

temps) on tends bien vers la solution exacte:

dy dy
yi 1  yi  t  O ( t ) yi 1  yi  t  O ( t )
dt ti dt ti1
dy dy
si t  0 les dérivées 
dt t i1 dt t i
donc les deux méthodes tendent vers le même résultat

de plus O( t )  0
donc les deux méthodes tendent vers la solution exacte

44
 Stabilité

Un méthode est dit stable si l’erreur n’augmente pas démesurément au

cour du calcul, c’est à dire jusqu’à donner des valeurs impossible.
Exemple:
t
dy y 
  et à t  0 , y  y0 solution  y (t )  y0 e 
dt 

Si y0>0 initialement tout les yi seront positifs, une méthode

conduisant à des y<0 serait instable.
Appliquer Euler explicite et implicite et discuter la stabilité en
fonction de Δt/τ

45
 Crank-Nicholson
dy (t )
Soit à résoudre:  F ( y, t )
dt
dy
yi 1  yi  t  O ( t )
dt ti
dy
Développement de Taylor: yi 1  yi  t  O (t )
dt ti1

t  dy dy 
yi 1  yi    O(t 2 )
2  dt ti dt ti1 
Soit h  t , on passe de i à i  1suivant :
ti 1  ti  h

 h
y
 i 1  y i  F ( yi , ti )  F ( yi 1 , ti 1 )   O(2) 46
2
 Prédicteur-Correcteur
dy (t )
Soit à résoudre:  F ( y, t ) et F non linéaire en Y
dt
On veut utiliser un méthode implicite, ou contenant une partie
implicite.
Il faut connaître: F ( yi 1 , ti 1 )
1 - On évalue yi*1à partir d' une méthode explicite : (prédicteur)
yi*1  yi  hF ( yi , ti )  F ( yi 1 , ti 1 )  F ( yi*1 , ti 1 )
2 - On utilise la méthode implicite : (correcteur)
h

yi 1  yi  F ( yi , ti )  F ( yi*1 , ti 1 )
2

On peut inclure un processus itératif pour faire converger les étapes
47
1et 2
 Méthode Runge-Kutta
dy (t )
Soit à résoudre:  F ( y, t )
dt
On par de (ti,yi):

 1  hF ( yi , ti )
 h
 2  hF ( yi  1 , ti  )
2 2
2 h
 3  hF ( yi  , ti  )
2 2
 4  hF ( yi   3 , ti  h)
1
yi 1  yi   1  2 2  2 3   4 )  O (4) 
6
48
 Résolution d’un système d ’équations
La résolution de systèmes d’équations est utilisée dans de très
nombreux domaines et les champs d’application sont très vastes.

• mathématique (par définition)

• physique (méca des fluides, des structures,
échelle atomique jusqu’à astronomique, …)
• chimie, biologie, … sciences en général
• informatique (réseaux, bases de données, …)
• économique et statistique (modèles financiers,
prévisionnels, …)
49
 Système à résoudre
Le système à résoudre peut s’écrire sous la forme suivante :

a11 X 1  a12 X 2    a1nXn  Y 1 qui se résume en

a 21 X 1  a 22 X 2    a 2 nXn  Y 2 AX  Y
 avec
an1 X 1  an 2 X 2    annXn  Yn X  X 1, X 2,  , Xn, Y  Y 1, Y 2,  , Yn
et

 a11 a12  a1n 

a 21 a 22  a 2n
A 
   
 an1 an 2  ann 

50
 Méthode de Cramer

La méthode consiste à utiliser le calcul de déterminants.

xi  det( Ai ) / det( A)

det( A) est le déterminant de la matrice A.

det( Ai ) est le déterminant de la matrice A mais la i ème colonne de la matrice
a été remplacée par le vecteur colonne Y.

A éviter :
(Dujardin)
Cramer Gauss
n=5 0.001 s 0.003 s
n=6 0.006 s 0.005 s
n=10 30 mn 0.025 s
n=12 1 heure 0.043 s
n=200 600 siècles…0.29 s

51
 Méthode de Gauss-Jordan

Il s’agit de diagonaliser la matrice du système par

combinaison linéaire des lignes ou des colonnes:
 a11 a12  a1n   x1   y1 
a 21 a 22  a 2n  x2  y 2
     
      
 an1 an 2  ann   xn   yn 


 1 0  0   x'1   y '1 
 0 1  0   x'2   y '2 
     
       
 0 0  1   x' n   y ' n 
    
52
 Méthode de Gauss-Jordan
Appliquons la méthode pour un système 4x4:
 a11 a12 a13 a14   x1   y1 
a 21 a 22 a 23 a 24   x 2   y 2  Les lignes sont prises les unes après
       les autres, le pivot est le premier
 a 31 a 32 a33 a 34   x3   y 3  élément non-nul de la ligne.
a 41 a 42 a 43 a 44   x 4   y 4 


1 a'12 a'13 a'14   x1   y '1

a 21 La première ligne est donc divisée
a 22 a 23 a 24   x 2   y 2 
       par le pivot a11
 a 31 a32 a33 a34   x3   y 3 
a 41 a 42 a 43 a 44   x 4   y 4 


53
 Méthode de Gauss-Jordan
1 a'12 a'13 a'14   x1   y '1  La première colonne des autres
0 a '22 a'23 a '24   x 2  y '2  lignes est annulée en soustrayant ces
     
0 a '32 a'33 a'34   x 3   y '3  lignes avec la première ligne
0 multipliée par les termes à éliminer.
 a '42 a'43 a '44   x 4  y '4 
Ln  Ln  L1 * an1
1 a'12 a'13 a '14   x1   y '1 
0 1 a' '23 a' '24   x 2  y ' '2 La seconde ligne est à présent
     
0 a '32 a'33 a'34   x 3   y '3  divisée par son pivot.
0 a '42 a'43 a'44   x 4  y '4 

1 0 a' '13 a' '14   x1   y ' '1  La seconde colonne des autres
0 1 a' '23 a' '24   x 2   y ' '2 
      lignes est annulée en soustrayant ces
0 0 a' '33 a' '34   x3   y ' '3  lignes avec la seconde ligne
0 0 a' '43 a' '44   x 4   y ' '4  multipliée par les termes à éliminer.

54
 Méthode de Gauss-Jordan

iv
1 0 0 0  x1   y 1 
0     iv  On applique le pivot de la même
1 0 0 x2  y 2
    manière sur les 3ème et 4ème
0 0 1 0  x3   y iv 3 
0      lignes.
 0 0 1  x 4   y 4 
iv

iv
 x1   y 1
Le système est résolu.
 x 2   y iv 2 
  
 x3   y iv 3 
 x 4   iv 
   y 4 

55
 Méthode de Gauss-Jordan
Problèmes éventuels de résolution :

•Le pivot peut être nul ce qui rend la division impossible.

•Le pivot peut être très petit ce qui peut rendre la
solution instable numériquement.
•Le pivot se vectorise ou se parallélise difficilement

Solution optimale :

Intervertir les lignes de manière à ce que le pivot soit

maximal (en valeur absolue) sur la ligne considérée.

56
 Méthode de Gauss-Jordan
Exercice : décrire la procédure mathématique de la méthode

57
 Méthode de Gauss-Jordan
Procédure, résolution de : AX  Y
ai(,0j)  A(0), i et j  1, n
 (0)
ai ,n 1  yi
pour k  0 à n  1 :

i  k 1  ( k 1) ( k 1) (k )
  a k 1, j  a k 1, j a k 1, k 1
j  k  1 à n  1
i  1 à n, i  k  1 ( k 1) (k ) (k ) ( k 1)
  ai, j  ai, j  a .a
i , k 1 k 1, j
j  k 1 à n 1 
xi  ai(,nn)1

58
 Inversion de matrice

Soit à résoudre le système: AX  Y

1 1
Par définition : AA  A AI

Donc: A
1
 AX  Y   IX  A 1
YXA Y 1

Si on réécrit le système : AX  IY
En transformant: AX  IX  X
On aura aussi: IY  A 1Y

Ainsi on obtiendra la matrice inverse en appliquant à la matrice

unité les opérations que l’on applique à A pour la rendre unité
59
 Inversion de matrice
2 1 1
 
Soit A une matrice à inverser : A  1 2 1
1 1 2
 
On associe à A la matrice unité :
 2 1 1  1 0 0 
  
AI   1 2 1  0 1 0 
 1 1 2  0 0 1 
  
La méthode de Gauss-Jordan au deux matrices simultanément:

1 1 2 1 2 1 2 0 0   1 1 2 1 2  1 2 0 0 
     
1 2 1  0 1 0   0 3 2 1 2   1 2 1 0 
1 1 2  0 0 1   0 1 2 3 2   1 2 0 1 
     
60
 Inversion de matrice

Finalement, on obtient :
 1 0 0  3 4  1 4  1 4 
  
 0 1 0   1 4 3 4  1 4 
 0 0 1   1 4  1 4 3 4 
  
soit
 3 4 1 4 1 4 
1
 
A   1 4 3 4 1 4 
 1 4 1 4 3 4 
 

Vérifier que A-1 est bien la matrice inverse ….

61
 Méthode de Gauss
Il s’agit de trianguler la matrice du système afin de réduire le
temps calcul lié à l ’annulation des termes au dessus du pivot.

 a11 a12  a1n   x1   y1 

a 21 a 22  a 2n  x2  y 2
     
      
 an1 an 2  ann   xn   yn 


 1 a'12  a '1n   x'1   y '1 

0 1  a'2 n   x'2   y '2 
     
        
0 0  1   x' n   y ' n 

62
 Méthode de Gauss
1 a'12 a'13 a '14   x1   y '1 
0 1 a' '23 a' '24   x 2  y ' '2 Les lignes sont divisées par leurs
     
0 a '32 a'33 a'34   x 3   y '3  pivots.
0 a '42 a'43 a'44   x 4  y '4 

1 a ' '12 a' '13 a ' '14   x1   y ' '1  Toutes les valeurs sous les pivots
0 1 a' '23 a' '24   x 2   y ' '2  sont mises à 0
     
0 0 1 a' '34   x3   y ' '3 
0 0 0 1   x 4   y ' '4 

 x 4  Y ' '4
 x3  y' '3  a' '34 x 4 La solution est trouvée par
 inversion du système

 x 2  y ' '2  a' '24 x 4  a' '23 x3 triangulaire.
 x1  y ' '1  a' '14 x 4  a' '13 x3  a' '12 x 2
63
 Décomposition LU
La méthode de Gauss réduit une matrice complète en une matrice triangulaire. La
procédure peut être effectuée de manière plus formelle en exprimant la matrice A
originale comme le produit d’une matrice triangulaire inférieure et d’une matrice
triangulaire supérieure. A peut donc être décomposée, dans quasiment tous les cas,
selon l ’expression suivante :

A  LU
Où L est une matrice triangulaire inférieure (lower) et U une matrice triangulaire
supérieure (upper). Afin que cette factorisation soit unique, il est possible d ’imposer
le fait que les éléments diagonaux de L ou U soient unitaires.
U est la matrice obtenue par une décomposition de Gauss et L est la matrice composée
des facteurs multiplicatifs utilisés dans la décomposition de Gauss. De plus, une
matrice unique (A) peut être utilisée pour stocker L et U.

64
 Décomposition LU
Le système à résoudre étant : AX  Y et sachant que : A  LU

Il s ’agit de déterminer au préalable Z tel que : LZ  Y

Avant d’en déduire X de l ’expression : UX  Z

Pourquoi utiliser une décomposition LU et non pas la méthode de Gauss ?

Parce que ...
La décomposition LU s ’effectue sans intervention au niveau du vecteur Y qui
peut être quelconque. Cela signifie donc que plusieurs systèmes ayant la même
matrice A mais des parties droites Y différentes peuvent être résolus.

65
 Matrices tridiagonales
Système couramment rencontré

 a11 a12 0 0   x1   y1   b1 c1 0 0   x1   y1 
a a22  0   x2   y 2  a b  0   x   y 
 21       2 2  2    2 
0   an 1n      0   cn 1    
0
 0 ann 1 ann   xn   yn  0 0 a
 n b x   y 
n  n   n

Économie de mémoire

a1 , a2 ,, an , b1 , b2 ,, bn , c1 , c2 ,, cn 

avec a1  cn  0

66
 Matrices tridiagonales - Choleski

Il est possible d’utiliser la méthode de Gauss.

Double balayage de Choleski : une descente et une remontée

 b1 x1  c1 x2  y1
a x  b2 x2  c2 x3  y2
 2 1
  
a x  bi xi  ci xi 1  yi
 i i 1
an xn 1  bn xn  yn
67
 Matrices tridiagonales - Choleski
Chaque élément xi peut être exprimé en fonction de xi 1
c1 y1
x1   x2 
b1 b1
 c1a2  y1
pour la ligne suivante :   b2  x2  c2 x3  y2  a2
 b1  b1

De manière générale, la descente s’écrit pour i = 1 à n:


x  A x  B
i i 1
 i i

 ci  A0  0
 Ai   de plus i  1  
 ai Ai 1  bi  B0  0
 yi  ai Bi 1
B
 i 
 ai Ai 1  bi 68
 Matrices tridiagonales - Choleski
Les binômes suivants sont donc connus :
( A0 , B0 ), ( A1 , B1 ), ( A2 , B2 ), , ( An , Bn )
La remontée se fait de n à 1 .
an xn 1  bn xn  yn
an  An 1 xn  Bn 1   bn xn  yn
yn  an Bn1
xn   Bn
an An 1  bn
 xn1  An 1 xn  Bn 1
x  
 n 2


 x1  
69
 Méthode itérative

Les méthodes décrites précédemment (ainsi que d’autres qui sont similaires) sont
utilisables dans tous les cas de figure. Cependant, il y a des situations où il peut être
intéressant d’utiliser une méthode itérative. Notamment dans le cas des matrices
creuses (sparse matrices) qui sont de très grandes dimensions et qui ont par définition
de nombreux éléments nuls (cela signifie que l ’on a un système avec beaucoup
d ’équations qui contiennent elles-même très peu d ’inconnues). Dans ce cas, il peut
être plus économique d ’utiliser une méthode itérative comme celle de Jacobi par
exemple. Il s ’agit alors d ’initialiser les calculs avec une approximation de la solution
puis d ’utiliser chaque itération pour se rapprocher de la solution exacte. Si le nombre
d’itérations reste raisonnable, cette méthode est la plus économique.

70
 Méthode itérative - Jacobi

Reprenons le système AX  Y (A matrice, X et Y vecteurs)

a11 X 1  a12 X 2    a1nXn  Y 1
a 21 X 1  a 22 X 2    a 2 nXn  Y 2


an1 X 1  an 2 X 2    annXn  Yn

Isolons les inconnues de la manière suivante :

X 1  (Y 1  a12 X 2    a1nXn ) / a11
X 2  (Y 2  a 21 X 1    a 2 nXn ) / a 22

Xn  (Yn  an1 X 1    ann  1 Xn  1) / ann 71
 Méthode itérative - Jacobi
En choisissant de manière arbitraire un premier vecteur
solution : X 0  ( X 10 , X 20 ,, X n0 )

Il est alors possible d ’obtenir le vecteur solution du niveau

supérieur : X 1  ( X 11 , X 21 ,, X n1 )
X 1t 1  (Y1  a12 X 2t    a1n X nt ) / a11
X 2t 1  (Y2  a21 X 1t  a2 n X nt ) / a22

X nt 1  (Yn  an1 X 1t  ann 1 X nt 1 ) / ann
Et ainsi de suite, nous obtenons X 2 , X 3 , , X 

72
 Méthode itérative - convergence
Le calcul se termine lorsque le résidu (différence de deux
solutions successives) atteint une valeur minimale R  
Le résidu peut se calculer de manières différentes :

n t 1 t
 i X  X i
 i 1
 n t 1 t 2
  Xi  Xi
 i 1
R n
X it 1  X it

 i 1 X it 1

n t 1 t 2
 Xi  Xi
  X t 1
 i 1 i
73
 Méthode itérative - Gauss-Seidel

La méthode de Gauss Seidel est quasiment identique à celle

de Jacobi à une différence près : dès qu ’un élément de X t 1
est obtenu, il est introduit dans la résolution sans attendre le
pas suivant :

X1t 1  (Y1  a12 X 2t    a1n X nt ) / a11 X1t 1  (Y1  a12 X 2t    a1n X nt ) / a11
X 2t 1  (Y2  a21 X1t 1    a2n X nt ) / a22 X 2t 1  (Y2  a21X1t    a2n X nt ) / a22
 
t 1 t t
X nt 1  (Yn  an1 X1t 1    ann1 X nt 11 ) / ann X n  (Yn  an1 X1    ann1 X n1 ) / ann
 
Gauss - Siedel Jacobi

74
 Convergence

Pour la méthode de Gauss Seidel la convergence est assurée si

le système est à matrice diagonale dominante:

n
aii   aij , i
j 1
j i

75
 Facteur de sous relaxation

On utilise généralement un facteur de sous-relaxation soit pour

améliorer la vitesse de convergence soit pour aider la
convergence:
  0;2 : facteur de sous relaxation
si   0;1  favorise la convergence
si   1;2  augmente la vitesse de convergence
X 
t 1
i final
t
 
 X  X
i
t 1
i approché
X t
i 
76
 Systèmes linéaires - Conditionnement
Reprenons deux systèmes : AX  Y
 x1  x2  1  x1  0
 
1.002 x1  x2  1  x2  1
 x1  x2  1  x1  0.5
 
1.002 x1  x2  1.001  x2  0.5
Un tel écart dans les solutions pour une différence d ’1/1000
du système signifie que l ’on a un mauvais conditionnement.
Les données peuvent provenir de mesures expérimentales,
un mauvais conditionnement peut donner de nombreux
résultats différents dans le champ d’incertitudes des
mesures + problème arrondit machine. 77
 Systèmes linéaires - Conditionnement
Critère de conditionnement (Ralston) :
Un système est probablement mal conditionné si les
éléments de la matrice inverse A1 sont différents de
plusieurs ordres de grandeur de la matrice A .
En effet :
 1  1  500 500
A  et 1
A  
1.002  1   501 500
Ce critère ne peut être utilisé si A possède une dimension
physique car A1 n ’aurait pas la même dimension et ne
pourrait lui être comparée.
78
 Systèmes linéaires - Conditionnement
Critère général de conditionnement :
Si la jème colonne de la matrice A est notée a j et a pour
composantes ( A1 j , A2 j , , Anj ) . Le facteur de conditionnement
est défini par :

det A

a1 a2  an

Ce facteur est compris entre 0 et 1.

Si   1 alors le système est mal conditionné.

79
 Equation de la chaleur
2
Soit à résoudre: d y ( x)
2
 F ( x, y )
dx
Le long d’une barre de longueur L divisé en n intervalle Δx .

80
 Racines d ’un polynôme
Méthode de Bairstow
Le but du jeu est de calculer les racines d’un polynôme d ’ordre n>2.
(1) Pn ( x)  a0 x n  a1 x n 1    an 1 x  an
Quelle que soit la valeur de n, il est possible d’écrire:
(2) Pn ( x)  x 2  px  q b0 x n  2  b1 x n 3    bn 3 x  bn 2   Rx  S
A chaque couple (p,q) correspond un ensemble de valeur b,R et S soit :
b0  b0 ( p, q), , R  R( p, q), S  S ( p, q)
Il s ’agit maintenant de trouver le couple (p,q) tel que R et S soient nuls, dans ce cas les
racines du polynôme d ’ordre 2 x 2  px  q  0 sont aussi racines du polynôme P.

Lorsque cette étape est franchie, il nous faut trouver les racines du polynôme d’ordre n-2
ce qui est fait en utilisant la même méthode.

81
 Intégration de fonction
Singularités de l’intégrande
•La première chose à faire est d ’essayer d ’éliminer la singularité
de manière mathématique : intégration par partie, changement de
variable,…
•Si on calcule l ’intégrale par une méthode de Simpson, l ’intervalle
doit exclure la borne singulière.

1
1
I  dx avec  très petit

x

82
 Racines d ’un polynôme
Méthode de Bairstow

Déterminons les facteurs b, R et S en fonction de p et q. Il s ’agit tout simplement de

développer l’expression de la relation (2) et de l’identifier à (1).
Posons tout d ’abord  R  bn 1
(3) 
S  pbn1  bn
Par identification, nous obtenons alors

(4)
bk  ak  pbk 1  qbk  2

b 2  b1  0
Les relations (4) permettent d ’obtenir par récurrence bk en fonction de p et q ainsi que
R et S. Il s ’agit à présent de trouver P et Q tels que R=S=0. Pour cela, on utilise une
méthode récursive (Newton-Raphson).

83
 Racines d ’un polynôme
Méthode de Bairstow
 R ( p0 , q0 )  R0
Pour un p0 et un q0 arbitraires, il est possible de calculer 
 S ( p0 , q0 )  S 0
 R ( p 0  p , q 0  q )  0
Il faut alors déterminer p et q tels que 
 S ( p 0  p , q 0  q )  0
0  R0  p R p 0  q R q 0
Avec un D.L. d ’ordre 1 : 
0  S 0  pS p 0  q S q 0

0  R0  p bn 1 p 0  q bn 1 q 0

Soit, d ’après la relation (3)
0  S 0  p bn 1  p0 bn 1 p 0  bn p 0  
 q  p0 bn 1 q 0  bn q 0 

84
 Racines d ’un polynôme
Méthode de Bairstow

Calculons les différentes dérivées partielles grâce au système (4)

bk p  bk 1  p(bk 1 p )  q bk 2 p

Posons ck 1   bk p pour k  1, n  1 alors

ck  bk  pck 1  qck  2

(5) k  0, n  1 
c1  c 2  0
bk q  bk 2  p(bk 1 q)  q bk  2 q
Posons ck 2   bk q pour k  2, n alors

ck  bk  pck 1  qck  2

(6) k  0, n  2

c1  c 2  0
85
 Racines d ’un polynôme
Méthode de Bairstow

Il s ’avère que ck  ck (résultat de la comparaison entre (5) et (6)). De ce fait,

p  P D
q  Q D

2
 D  c n  2  cn 1cn  3
P  b c  b c
n 1 n  2 n n 3

Q  bn cn 2  bn 1cn 1

86
 Racines d ’un polynôme
Méthode de Bairstow, procédure.

1) Nous avons pn (x), polynôme d ’ordre n. Choisissons p0 et q0 arbitrairement.

2) Détermination de la famille des bn. (Par identification; relation 4)
3) Détermination de la famille des cn. (Relation 5)
4) Calcul de p1  p0  p et q1  q0  q . On boucle sur l’étape 1 avec
cette fois ci p0=p1 et q0 = q1. On stoppe la boucle lorsque R   et S  
5) On résout le polynôme x 2  p j x  q j  0 qui nous donne deux racines de
pn(x).
6) On recommence l ’ensemble du problème avec pn-2(x) et ainsi de suite jusqu ’à
obtenir un polynôme d ’ordre 2 ou 1.

87