Transcript module doctoral

Calcul Scientifique pour les
Sciences de l’Ingénieur
Bruno Koobus
Franck Nicoud
MOTIVATION
MD "Calcul Scientifique"
2
Objectifs
• Techniques et outils de base
• Culture de base sur le calcul scientifique Exemples
• Résolution d’un problème académique
MD "Calcul Scientifique"
3
Motivation
• Les équations des sciences de l’ingénieur sont connues
depuis longtemps …
– Eq. de la chaleur (Fourier, 1807)
– Mécanique des fluides (Navier-Stokes, 1822)
– Elecromagnétisme (Maxwell, 1873)
• Les inconnues sont des fonctions scalaires ou
vectorielles de plusieurs variables: x,y,z,t
• Ces fonctions sont solutions d’équations aux dérivées
partielles (EDP) en général non-linéaires
MD "Calcul Scientifique"
4
Motivation
• Les cas où une solution analytique peut
être trouvée sont très rares et simples:
– géométrie simple
– équation linéaire
• Les solutions obtenues sont souvent très
lourdes même pour ces cas simplistes …
MD "Calcul Scientifique"
5
Exemple de résolution analytique
• Equation de convection diffusion 1D dans
un domaine infini
C( x,0)  C0 ( x)
lim C ( x, t )  0
x  
U0
C ( x, t ) ?
C
C
 2C
 U0
D 2
t
x
y
MD "Calcul Scientifique"
6
Exemple de résolution analytique
• Si la concentration est initialement de la
forme C(x,0)  C ( x)  C exp x / 4a 
on peut
obtenir la solution analytique …
2
0
C ( x, t )  C0
2
0
2


x  U 0t 
a
exp 
2
2
a  Dt
 4 a  Dt
2

MD "Calcul Scientifique"





7
Exemple de résolution analytique
• Effet de la diffusion
Re  
Re  20
U0  a
Re 
D
Re  200
Re  2
MD "Calcul Scientifique"
8
Exemple de résolution analytique
• Equation de la chaleur dans une cavité
rectangulaire
y

H
T
( x, H )  h  T ( x, H )  T0   0
y
T ( x, y) ?
T (0, y)  T0
0
 2T  2T
K
 2 
2
x
y

0
T ( x,0)  T0
MD "Calcul Scientifique"
T ( L, y)  T0
L
x
9
Exemple de résolution analytique
• Une méthode de séparation des variables
permet d’obtenir la solution analytique de
ce petit problème d’école …



 



k

k





 k
k

T ( x, y )  T0    Ak sinh
y 

cosh
y

1


sin

 

2 


L
L

  k  

 
 L
k 1




L




 k  hL 
 k  
H
H   1
 cosh
k L
 L  k 
 L  
Ak 
k
k
 k 
 k 

cosh
H   h sinh
H
L
L
L




 sinh
MD "Calcul Scientifique"


x


2K0
 1k  1
k 
k
10
Exemple de résolution analytique
Effet des fuites par convection
h0
h  0.05
h  0.01
MD "Calcul Scientifique"
h 1
11
Passage à des configurations plus
complexes …
• Les cas académiques sont très utiles pour
comprendre les phénomènes de bases, étudier
les effets des différents termes des équations
• Les problèmes posés en pratique sont
inaccessibles à la résolution analytique
• Il est alors naturel de faire des expériences …
MD "Calcul Scientifique"
12
Expériences …
• On imagine que l’on veut fabriquer un avion
• Il est exclu de réaliser toutes les mises au point
de la forme extérieure à partir de réalisations à
échelle 1
MD "Calcul Scientifique"
13
Expériences …
• On utilise donc des maquettes de taille réduites que l’on
met dans une soufflerie
Vitesse de l’avion
ou de l’air dans la
soufflerie
Mach ≈ 0.9 < 1
Re 
U0  L

Taille de la maquette:
entre 1% et 10% du vrai
avion
Viscosité cinématique
de l’air.
En gros 1.5 105 m2 / s
• Le nombre de Reynolds ne peut que très difficilement
être respecté (10 à 100 fois trop petit)
• Extension des données acquises au cas du vrai avion ?
MD "Calcul Scientifique"
14
Autre exemple …
Propergol solide
MD "Calcul Scientifique"
15
Expériences: de moins en moins !
• La tendance actuelle est à la diminution des expériences
à échelle 1 pour réduire les coûts de fabrication
• Les expériences sont réservées aux dernières mises au
point avant la certification et la mise en production
• Les phases de mise au point préliminaires font de plus
en plus appel à des « souffleries numériques »
• L’idée est de résoudre les EDP régissant le
fonctionnement du système numériquement puisque
l’approche analytique est impossible
MD "Calcul Scientifique"
16
Simulation numérique
• Lorsque l’on résout les EDP analytiquement, on cherche les
solutions dans des espaces fonctionnels de dimension infinie
(e.g. L2)
• Si on y arrive, on a accès à la solution (e.g. la vitesse de l’air autour
de l’avion) pour tout point de l’espace, et à chaque instant

V x, y, z, t , x, y, z, t    0,
• Puisque l’on y arrive presque jamais, on décide de simplifier le
problème en ne cherchant la solution qu’en un nombre fini de points
de l’espace et pour un nombre fini d’instants

V xi , yi , zi , tn , i, n1,2,...,N1,2,...,M 
• On cherche alors les solutions de ce problème discrétisé dans un
espace vectoriel de dimension finie.
MD "Calcul Scientifique"
17
SIMULATION
NUMERIQUE
MD "Calcul Scientifique"
18
Maillages: profil d’aile d’avion
hybride
non structuré
structuré
MD "Calcul Scientifique"
19
Maillages: pot d’échappement
MD "Calcul Scientifique"
20
Maillages: injecteur + chambre
MD "Calcul Scientifique"
21
Maillages: injecteur + chambre
structuré
Non-structuré
MD "Calcul Scientifique"
22
Exemple de flamme turbulente
stable
Simulation aux Grandes Echelles – CERFACS – A. Sengissen
MD "Calcul Scientifique"
23
Exemple en biomédical
Iliac bifurcation – CHU Toulouse
3D model from
CT images
artheriography
MD "Calcul Scientifique"
24
FUNCTIONAL IMAGING
Iliac Bifurcation
Computational Grid
Wall Shear Stress
MD "Calcul Scientifique"
Tracers
25
Trois familles de méthodes
• Éléments finis
(B. Koobus)
• Volume finis
(non abordés dans ce module)
• Méthodes aux différences finies
(F. Nicoud)
MD "Calcul Scientifique"
26
Eléments finis en quelques mots
• A chaque instant, on cherche la solution de
l’EDP sous la forme
N


f h ( x )   f ii ( x )
i 1

fonctions i (x ) forment
• Les
une base de l’espace
de dimension N dans lequel on cherche à
approximer la ‘vraie’ solution f par fh.
• Les coefficients fi sont déterminés en imposant à
fh d’être la meilleure approximation de f dans
l’espace de dimension N choisi.
MD "Calcul Scientifique"
27
Éléments finis en quelques mots

• Un choix classique est de prendre i (x ) linéaire
par morceaux et égale à 1 au nœud i du
maillage, nulle partout ailleurs
0
1
Nulle
Linéaire
MD "Calcul Scientifique"
28
Éléments finis en quelques mots
K


ex
:

f


0





• On cherche à résoudre E(f)=0
• Avec l’approximation
commet une erreur E(fh)
N


f h ( x )   f ii ( x )
i 1
on
• La méthode de Galerkin consiste à dire que
cette erreur et orthogonale aux fonctions de
forme i (x)
  N
  
k  1,2,...,N ,  k ( x )E  f ii ( x ) dx  0
 i 1


• N équations, N inconnues …
MD "Calcul Scientifique"
29
Volumes finis en quelques mots
• Bien adaptés à des problèmes conservatifs du
type



f
 div F ( f )  0
t
• On intègre l’équation sur chaque cellule du
maillage et on utilise le théorème de la
divergence
n 1
n

Vi
fi

 fi
   Fk  nk dSk
t
faces k
de Vi
• Les inconnues sont les valeurs moyennes de f
sur chaque cellules Vi à l’instant n+1, soit fi n1
MD "Calcul Scientifique"
30
Volumes finis en quelques mots
• Les flux sur les faces de Vi sont calculés à partir
des valeur de f dans les cellules voisines

n3

n1
Vi
 
f i n 1  f i n
Vi
   Fk  nk dSk
t
faces k
de Vi

n2
MD "Calcul Scientifique"
31
Différences finies
• Contrairement aux éléments et volumes finis,
cette technique n’est pas adaptée aux maillages
non cartésiens
• Mais elle est très intuitive
• En 1D, les trois méthodes sont équivalentes
• Permet d’appréhender beaucoup de concepts
ou problèmes numériques communs aux
différentes méthodes
MD "Calcul Scientifique"
32
Différences finies
• L’idée est de remplacer les dérivées
partielles aux points de maillage par des
développement de Taylor
xi 3 xi  2
xi 1
i  3 i  2 i 1
xi
xi 1
xi  2
x
i
i 1 i  2
x
• Plutôt que de chercher f(x), on cherche les
valeurs de f aux nœuds du maillage, soit
fi=f(xi)
MD "Calcul Scientifique"
33
DERIVEES
PREMIERES
MD "Calcul Scientifique"
34
Dérivées premières
• Développement de Taylor au nœud i:
i  3 i  2 i 1
i
x
i 1 i  2

xi 1  xi  d 2 f
df
f i 1  f i  xi 1  xi 

dx xi
2
dx2
2

xi 1  xi  d 2 f
df
f i 1  f i  xi 1  xi 

dx xi
2
dx2

2


2

 o xi 1  xi 
xi
2
 o xi 1  xi 
xi
• Ces développements font apparaître les
dérivées de f au nœud i uniquement
MD "Calcul Scientifique"
35
Dérivées premières
• Si les nœuds sont régulièrement espacés
x
i  3 i  2 i 1
i
df
x 2 d 2 f
f i 1  f i  x

dx xi
2 dx2
df
x 2 d 2 f
f i 1  f i  x

dx xi
2 dx2
f i 1  f i 1  0  2x
x
i 1 i  2
 
 o x 2
xi
 
 o x 2
xi
 
df
 0  o x 2
dx xi
MD "Calcul Scientifique"
36
Dérivées premières
• Si les nœuds sont régulièrement espacés,
la dérivée de f au nœud i est approximée
par
f i 1  f i 1
df
0
 D1 f i 
dx xi
2x
• Erreur d’approximation est o(x)
• Schéma centré d’ordre 2
MD "Calcul Scientifique"
37
Dérivées premières
• On peut manipuler les développements
limités pour obtenir d’autres
approximations de la dérivée première
f i 1  f i
df

•
 D1 f i 
 o(1)
dx xi
x
f i  f i 1
df

 D1 f i 
 o(1)
•
dx xi
x
MD "Calcul Scientifique"
38
Maillage non uniforme
• Développement de Taylor au nœud i:
i 3
i  2 i 1 i i  1
x
i2

xi 1  xi  d 2 f
df
f i 1  f i  xi 1  xi 

dx xi
2
dx2
2

xi 1  xi  d 2 f
df
f i 1  f i  xi 1  xi 

dx xi
2
dx2
2


2

xi
2
i 1  xi  xi 1

 o xi 1  xi 
 o xi 1  xi 
xi
i  xi 1  xi
MD "Calcul Scientifique"
39
Maillage non uniforme
• Développement de Taylor au nœud i:
i 3
i  2 i 1 i i  1
x
i2
i d 2 f
df
f i 1  f i   i

dx xi
2 dx2
 
2

2
i 1

 o i
xi
 i 1 d 2 f
df
f i 1  f i   i 1

dx xi
2 dx2

2
i
2

 f
2
i 1 i 1
 f
2
i i 1

 fi 
2
i 1
 
 o  i 1

2
i 1
2
i

2



df
 o 2i 12i
dx xi
2i 1 f i 1  f i 2i 1  2i  2i f i 1
df

 o 
MD
dx xi
 i "Calcul
1 i Scientifique"
i 1   i


xi
    i   i 1
2
i
2

40
Ordres plus élevés
• Maillage régulier
• En conservant plus de termes dans les
développements on obtient les schémas à
l’ordre 4 et 6 suivants
•
 f i 2  8 f i 1  8 f i 1  f i 2
df

 o(x3 )
dx xi
12x
•
fi 3  9 fi 2  45 f i 1  45 fi 1  9 fi 2  fi 3
df

 o(x5 )
dx xi
60x
MD "Calcul Scientifique"
41
Formules décentrées
•
f i 1  f i
df

 o(1)
dx xi
x
ordre 1 aval
•
f i  fi 1
df

 o(1)
dx xi
x
ordre 1 amont
•
 f i 2  4 f i 1  3 f i
df

 o(x)
dx xi
2x
ordre 2 aval
•
f i 2  4 f i 1  3 f i
df

 o(x)
dx xi
2x
ordre 2 amont
MD "Calcul Scientifique"
42
Comparaison des schémas
• Problème modèle 1D: Eq. de convection
f
f
U0
 0,  2 m  x  8 m, U 0  1 m/s
t
x
• Conditions limites et initiale:
f (2, t )  f (8, t )  0
f ( x,0)  exp x / 4a , a  0.2 m
2
2
t  5s
t  0s
MD "Calcul Scientifique"
43
Test numérique
df i
 U 0 Df i  0, i
dt
• Équation semi-discrète
f i 3
fi2
f i 1
fi
i  3 i  2 i 1
i
fi2
x
i 1 i  2
x
f i 1
• On calcule les fi entre t=0 et t=5 s à partir
des deux schémas
df i
f i  f i 1
U0
0
dt
x
amont ordre 1
df i
f i 1  f i 1
U0
0
dt
2x
MD "Calcul Scientifique"
centré ordre 2
44
Test numérique
400 noeuds
200 noeuds
100 noeuds
amont
ordre 1
centré
ordre 2
MD "Calcul Scientifique"
45
Ordre 2 centré / Ordre 1 amont
• Ordre 1 introduit de la diffusion … (cf
solution analytique avec Re=2)
• Ordre 2 centré « exact » avec 400 points
• Ordre 2 centré déforme le signal si le
nombre de points est plus petit
• Ordre 2 meilleur que ordre 1
MD "Calcul Scientifique"
46
Test numérique
400 noeuds
200 noeuds
100 noeuds
centré
ordre 4
centré
ordre 2
MD "Calcul Scientifique"
47
Ordre 2 centré / Ordre 4 centré
• Ordre 4 « exact » dans tous les cas
considérés ici
• Ordre 2 centré « exact » avec 400 points
• Ordre 4 meilleur que ordre 2
MD "Calcul Scientifique"
48
Test numérique
400 noeuds
200 noeuds
100 noeuds
amont
ordre 2
centré
ordre 2
MD "Calcul Scientifique"
49
Ordre 2 centré / Ordre 2 amont
• Ordre 2 centré et amont « exacts » avec 400
points
• Ordre 2 centré et amont déforment le signal si le
nombre de points est plus petit, mais pas de la
même manière
• Ordre 2 amont amortit plus le signal
• L’ordre ne dit pas tout sur un schéma …
MD "Calcul Scientifique"
50
Analyse spectrale
• Cas d’une fonction harmonique
f ( x)  Reexp( jkx ) 
df
 Re jk exp( jkx )
dx
• Schéma centré d’ordre 2
fi 1  f i 1
df
f i  Reexp(jkix),

dx xi
2x
df
 sin(kx)

 Re  jk
exp( jkix)
dx xi
kx


• L’erreur commise est
sin( kx )
k x
MD "Calcul Scientifique"
51
Signification de kx
• Sinusoïde de période L décrite avec N
points
• x = L / N,
kx  0
(exact)
k = 2/L donc kx = 2 / N
kx 

4
kx 
MD "Calcul Scientifique"

2
kx  
52
Analyse spectrale
• Tout se passe comme si on résolvait l’équation
f
sin( kx) f
U0
0
t
kx x
• Les différentes longueurs d’onde ne se
déplacent pas à la même vitesse
centré
ordre 2
MD "Calcul Scientifique"
exact
53
Analyse spectrale
• Équation effective
SCHEMA
Centré
ordre 2
Amont
ordre 1
Amont
ordre 2
Centré
ordre 4
f
f
 U 0 E (kx)
0
t
x
ReE(kx)
ImE(kx)
sin( kx )
0
k x
sin( kx )
cos( kx )  1
k x
k x
 cos( 2kx)  4 cos( kx)  3
sin( kx)
2  cos( kx) 
kx
kx
sin( kx)
0
4 MDcos(
kScientifique"
x) 
"Calcul
54
3kx
Analyse spectrale
MD "Calcul Scientifique"
55
Lien avec l’ordre du schéma
• Dans la limite kx → 0, la vitesse de propagation tend
vers U0
• La vitesse avec laquelle l’erreur tend vers zéro dépend
de l’ordre du schéma
• Au voisinage de 0, Re(E(kx)) = 1+O((kx)n), avec n
l’ordre du schéma
• La partie imaginaire de l’erreur n’est pas directement
reliée à l’ordre du schéma
• Les schémas centrés sont non dissipatifs: Im(E(kx)) = 0
MD "Calcul Scientifique"
56
Dispersion
• La vitesse de propagation effective n’est égale à la vitesse
théorique que dans la limite kx → 0
• Une perturbation peut donc être propagée trop lentement ou
trop vite
jkx
jk ' x
e
, k  k ' ne sont pas propagées à la
et
• Les fonctions e
même vitesse en général
• Que se passe-t-il lorsque l’on convecte f (x) ?
MD "Calcul Scientifique"
57
Déformation du signal
• On peut décomposer cette fonction comme une somme
de fonctions harmoniques (en rendant f périodique
éventuellement)
jkx
ˆ
f ( x)   f k e
• La solution théorique après t s de simulation est
f ( x  U 0t )   fˆk e jk ( x U 0t )
jk ( x  E ( kx )U 0t )
jkx
e
e
• Numériquement le mode
devient
• La solution numérique est donc
g ( x  U 0t )  
gˆ k

e jk ( x U 0t )  f ( x  U 0t )
e jk (1 E ( kx )) U 0t
MD "Calcul Scientifique"
58
DERIVEES
SECONDES
MD "Calcul Scientifique"
59
Dérivées secondes
• Maillage régulier
• On utilise le fait que
d2 f
dx2

xi
d  df 
dx  dx  xi
• En appliquant l’opérateur D à D f
d2 f
2
dx
d2 f
2
dx
0
1
0
1 i 1
xi
0
D
0 0
1 f i 1  D f
 D1 D1 f i 
2x
D
0
1 i
 o(x)
f i  2  2 f i  f i 2
f 
 o(x)
2
4x
0, 2
2
i
xi
MD "Calcul Scientifique"
60
Problème de localité
• Si les nœuds sont régulièrement espacés
x
i  3 i  2 i 1
d2 f
dx2
xi
i
i 1 i  2
x
f i  2  2 f i  f i 2

4 x 2
x
La dérivée seconde approximée de cette fonction
est nulle !!
MD "Calcul Scientifique"
61
Dérivées secondes
• Déduire la dérivée seconde des
développements de Taylor
i  3 i  2 i 1
i
df
x 2 d 2 f
f i 1  f i  x

dx xi
2 dx2
df
x 2 d 2 f
f i 1  f i  x

dx xi
2 dx2
d2 f
dx2
x
i 1 i  2
xi
xi
x 3 d 3 f

6 dx3
x 3 d 3 f

6 dx3
 
 o x 3
xi
 
 o x 3
xi
f i 1  2 f i  f i 1
D f 
 o(x)
2
x
0 ,1
2
i
xi
MD "Calcul Scientifique"
62
Problème de localité
• Si les nœuds sont régulièrement espacés
x
i  3 i  2 i 1
d2 f
dx2
xi
i
i 1 i  2
x
f i 1  f i  f i 1

x 2
x
La dérivée seconde approximée de cette fonction
est non nulle, mais pas infinie …
MD "Calcul Scientifique"
63
Analyse spectrale
• Cas d’une fonction harmonique
f ( x)  Reexp( jkx) 

d2 f
2

Re

k
exp( jkx)
2
dx

• Schéma centré d’ordre 2 à 2
f i  Reexp( jkix),
d2 f
dx2
xi
d2 f
dx2

xi
f i 1  2 f i  f i 1
x 2
 cos(kx)  1

 Re 2
exp(
jki

x
)

x 2


• L’erreur commise est
1  cos( kx )
2
k 2 x 2
MD "Calcul Scientifique"
64
Analyse spectrale
• Schéma centré d’ordre 2 à 4
f i  Reexp( jkix),
d2 f
dx2
xi
d2 f
dx2

xi
f i  2  2 f i  f i 2
4x 2
 1 cos(2kx)  1

 Re 
exp(
jki

x
)

x 2
2

• L’erreur commise est
1  cos( 2kx)
2 k 2 x 2
MD "Calcul Scientifique"
65
Analyse spectrale
• Les erreurs sont réelles uniquement, donc
pas de convection numérique
2
4
MD "Calcul Scientifique"
66
Interprétation volumes finis
• Si les nœuds sont régulièrement espacés
i  3 i  2 i 1
d2 f
dx2

xi
i
i 1 i  2
df
df

dx xi1 / 2 dx xi1 / 2
x
x
f i 1  2 f i  f i 1

x 2
• Utile pour les coefficients de diffusivité variables
MD "Calcul Scientifique"
67
Retour sur le schéma amont ordre 1
• Rappel: ce schéma introduit beaucoup de
dissipation par comparaison avec le centré
d’ordre 2
f i  f i 1 f i 1  f i 1 x f i 1  2 f i  f i 1


• En effet:
x
2x
2
x 2
• Utiliser ce schéma revient donc à résoudre
f
f U 0 x  2 f
 U0

t
x 
2 x 2
D
avec un schéma centré d’ordre 2
MD "Calcul Scientifique"
68
Laplacien
x
j2
y
j 1
j
j 1
j2
i  3 i  2 i 1
f
xi

f i 1, j  2 f i , j  f i 1, j
x
2
i

x
i 1 i  2
f i , j 1  2 f i , j  f i , j 1
MD "Calcul Scientifique"
y
2
 o(x, y)
69
INTEGRATION
TEMPORELLE
MD "Calcul Scientifique"
70
Schéma en temps
• La solution est évaluée en une succession de
temps discrets
fi n  f ( xi , tn )
• Le pas de temps t est le plus souvent constant
xi  i  x
tn  n  t
• Itération n: passer de tn à tn+1
MD "Calcul Scientifique"
71
Schémas en temps classiques
• Même démarche que pour les dérivées en
espace: développements de Taylor en temps
• Euler explicite:
f
t
n
f f

t
n
n 1
 o(1)
• Euler implicite:
f
t
n
f n1  f n

 o(1)
t
MD "Calcul Scientifique"
72
Intégration en temps
• Problème aux valeurs initiales:
df
 K ( f , t ), f (t0 )  f 0
dt
• Algorithme exact
f
n 1
t n1
 f 
n

tn
df
dt  f n 
dt
t n1
 K ( f , t )dt
tn
• Comment estimer l’intégrale de K ?
MD "Calcul Scientifique"
73
Intégration en temps
• Méthode à un pas:
f n1  f n  t  ( f n , t n , t )
• Méthode d’ordre p si:
f
n 1
f
 ( f n , t n , t )  O t p
t
 
n
MD "Calcul Scientifique"
74
Augmentation de l’ordre
• Il suffit de prendre un développement de Taylor
plus grand …
2
2
p
p
df
d
f

t
d
f

t
f n1  f n  t  2
 ...  p
 o(t p )
dt
dt 2
dt p!
• Et de définir ( f n , t n , t ) comme


t
 ( f , t , t )  K  K  K K

2
n
n
n
n
t
n
f
n

 K ttn  K nft K n  ( K n K nf  K tn ) K f  K n K tfn  K n K nff
2
MD "Calcul Scientifique"

t 2
 ...
6
75
Passage pratique à l’ordre 2
• Plutôt que d’estimer des dérivées d’ordre
élevées, il est préférable de mieux choisir les
endroits où on évalue K
• On part nden
( f , t , t )  A1 K ( f n , t n )  A2 K ( f  t , t  t )


Kn
• Par identification jusqu’à l’ordre 1 inclus


t
K  K K K
 ( A1  A2 ) K n  t  K tn  t  K nf
2
1
Kn


A1  A2  1
2
2
n
n
t
n
f
n
MD "Calcul Scientifique"
76
Runge-Kutta d’ordre 2
• Schéma à deux « étapes »:
k1  K ( f n , t n )
 n k1t n t 
k2  K  f 
,t  
2
2

f n 1  f n  t  k 2
MD "Calcul Scientifique"
77
Passage pratique à l’ordre p
• On part de
k1  K ( f , t )
n
n
k 2  K  f  k1 21t , t   2 t 
n
n
...
 n i 1

n
ki  K  f   k j  ij t , t   i t 
j 1


p
( f n , t n , t )   Aj k j
j 1
• Par identification on obtient les coefficients
qui permettent d’obtenir l’ordre p
MD "Calcul Scientifique"
78
Passage pratique à l’ordre p
• Nombre de coefficients à fixer:
i
 ij
Ai
p  1 1  0
p( p  1)

p( p  1) / 2   2 p 
1
2

p

• Système non-linéaire sous déterminé
MD "Calcul Scientifique"
79
Runge-Kutta ordre 4
• Problème aux valeurs initiales:
k1  K ( f n , t n )
t 
 n t
n
k 2  K  f  k1 , t  
2
2

t 
 n t
n
k3  K  f  k 2 , t  
2
2


k 4  K f n  tk3 , t n  t
f
n 1

t
 f   k1  2k 2  2k3  k 4 
6
n
MD "Calcul Scientifique"
80
Autres Schémas en temps classiques
• Crank-Nicolson:
f
n 1
 
 
t
 f   K f n 1  K f n
2
n

• Adams-Bashforth
t
n 1
n
f  f   3K f n  K f n 1
2
   

• Runge-Kutta d’ordre p
f
n 1

t  K 


f 
k!
p
n k
n
k 0
MD "Calcul Scientifique"
81
ANALYSE
SPECTRALE
MD "Calcul Scientifique"
82
Facteur d’amplification
• Forme de la solution
f ( x, t )  fˆ (t )  e jkx
f ( xi , tn )  fi n  fˆ n  e jkxi
• Comment l’amplitude d’une perturbation de
nombre d’onde k évolue-t-elle en temps ?
ˆf n1  Aˆ  fˆ n
MD "Calcul Scientifique"
83
Stabilité Von Neumann
• On s’intéresse désormais à l’équation
complètement discrétisée, y compris pour le
terme temporel
•
Aˆ  1 : signal amorti
•
Aˆ  1 : exact
•
Aˆ  1 : instabilité
MD "Calcul Scientifique"
84
Exemple 1
• Convection pure:
f
f
U0
0
t
x
• Schéma centré ordre 2 en espace
• Euler explicite en temps
n
n
f

f
f i n 1  f i n  tU 0 i 1 i 1
2x
MD "Calcul Scientifique"
85
Exemple 1 - suite
• Pour une perturbation du type
jkxi
n
ˆ
fi  f  e
n
tU 0
ˆ
A  1 j
sin( kx)
x
• Ce schéma est (inconditionnellement) instable
k : Aˆ  1
MD "Calcul Scientifique"
86
Exemple 2
• Décentré amont d’ordre 1 en espace
• Euler explicite en temps
n
n
f

f
i 1
f i n 1  f i n  tU 0 i
x
tU 0
ˆ
1  cos( kx)  j sin( kx)
A  1
x
• Ce schéma est (conditionnellement) stable
Aˆ  1, k
tU 0
 CFL 
1
x
MD "Calcul Scientifique"
87
Exemple 3
• Centré ordre 2 en espace
• Runge-Kutta d’ordre 2 en temps
n
n
2
f

f
U

t
n 1
n
f i  f i  U 0 t i 1 i 1  0
2x
2
2
f i n 2  f i n f i n  f i n2

2x
2x
2x
2


CFL
ˆA  1  jCFL sin(kx) 
cos(2kx)  1

4


• Ce schéma est inconditionnellement instable
MD "Calcul Scientifique"
88
Exemple 3 - suite
• Centré ordre 2 en espace
• Runge-Kutta d’ordre 2 en temps
f
K  f   U 0
x


ˆ n  e jkxi
K
f
sin(kx)
ˆ
K
  jU 0
ˆf n  e jkxi
x
k ˆk

t
K
ˆ
A
k!
k 0
3
• Ce schéma est inconditionnellement instable
MD "Calcul Scientifique"
89
Exemple 4
• Centré ordre 2 en espace
• Runge-Kutta d’ordre 3 en temps
CFL  1.8
ˆk


t
K
Aˆ  

 k  0 k! 
3
k
sin( kx)
ˆ
K   jU 0
x
• Ce schéma est conditionnellement stable
Aˆ  1, k
 CFL  1.73
MD "Calcul Scientifique"
90
Exemple 5
• Centré ordre 2 en espace
• Runge-Kutta d’ordre 4 en temps
CFL  1.8
ˆk


t
K
Aˆ  

 k  0 k! 
4
k
sin( kx)
ˆ
K   jU 0
x
• Ce schéma est conditionnellement stable
Aˆ  1, k
 CFL  2.82
MD "Calcul Scientifique"
91
Exemple 6
• Diffusion pure:
f
2 f
D 2
t
x
• Schéma centré ordre 2 en espace
• Euler explicite en temps
n
n
n
f

2
f

f
i
i 1
f i n1  f i n  tD i 1
x 2
MD "Calcul Scientifique"
92
Exemple 6 - suite
• Pour une perturbation du type
jkxi
n
ˆ
fi  f  e
n
2tD
ˆ
cos( kx)  1
A  1
2
x
• Ce schéma est conditionnellement stable
Aˆ  1, k

tD
Fo 
 0.5
2
x
MD "Calcul Scientifique"
93
Zone de stabilité
• Dans le cas général d’une équation linéaire
df
 K ( f , t)
dt
• On définit l’opérateur associé à K

ˆ n  e jkxi
K
f
Kˆ 
fˆ n  e jkxi

• Le schéma induit une équation du type
k ˆk
p


t
K 
ˆ
ˆ
ˆ
A  P tK ,
ex : A  

 k  0 k! 
 
MD "Calcul Scientifique"
94
Zone de stabilité
• On trace la courbe Aˆ  1 dans le plan complexe
tKˆ , tKˆ 
RK1
RK2
RK3
RK4
MD "Calcul Scientifique"
95
Zone de stabilité
• Le schéma complet est stable ssi la trajectoire
de tKˆ est contenue dans la zone de stabilité
Convection pure
RK4
tKˆ   jCFL sin(kx)
RK3
CFL  2.82
RK2
Diffusion pure
tKˆ  2Fo cos(kx) 1
1
Fo 
2
MD "Calcul Scientifique"
RK1
96
Effet de l’ordre
Convection pure ordre 2
tKˆ   jCFL sin(kx)
 
 tKˆ  2.82
CFL  2.82
RK4
Convection pure ordre 4
CFL  8 sin(kx) 
ˆ


tK   j
6   sin(2kx) 
CFL  2.05
Convection pure ordre 6
CFL
ˆ
45 sin( kx)  9 sin( 2kx)  sin(3kx) 
tK   j
30
MD "Calcul Scientifique"
CFL  1.77
RK3
RK2
RK1
97
RESOLUTION
DE SYSTEMES
LINEAIRES
MD "Calcul Scientifique"
98
Positionnement du problème
• Toutes les méthodes numériques
conduisent à l’établissement de systèmes
linéaires
• On doit résoudre AX  B
• A: carrée, B: vecteur colonne, X: inconnu
• Taille N: O(10)-O(106)
• Pour les méthodes DF, EF et VF, les
matrices sont souvent très creuses
MD "Calcul Scientifique"
99
Méthodes directes
• On connaît le compte exact d’opérations
• On obtient la solution exacte si les
ordinateurs étaient infiniment précis
• En pratique, c’est moins clair …
MD "Calcul Scientifique"
100
Méthodes directes
• Gauss-Jordan
– plusieurs variantes (pivot)
– nécessite environ N3 opérations (1 x et 1 +)
• Elimination de Gauss et substitution
– nécessite environ N3/3 opérations
• Décomposition LU
– Nécessite du « pivoting »
– nécessite environ N3/3 opérations
– 1 résolution de système triangulaire
demande environ N2/2 opérations
MD "Calcul Scientifique"
101
Méthodes directes
• Des méthodes directes spécifiques à
certains types de matrices existent
• Décomposition Cholesky
– A symétrique définie positive
– A=LLt
– nécessite environ N3/6 opérations (2 fois
moins que LU)
MD "Calcul Scientifique"
102
Méthodes itératives
• On ne connaît pas le compte exact
d’opérations
• On obtient la solution exacte si la
convergence complète est atteinte
• Permet un gain de temps énorme pour les
gros systèmes …
MD "Calcul Scientifique"
103
Méthodes itératives
• Jacobi
AX  B; A  L  D  U
DX k 1  ( L  U ) X k  B
• Rayon spectral: plus grande valeur propre
1

D
(L  U )
de la matrice d’itération
• Nombre d’itérations pour diviser le résidu
par 10 de l’ordre  ln10 / ln  Jacobi
• Laplacien en DF d’ordre 2: N2/2
MD "Calcul Scientifique"
104
Méthodes itératives
• Gauss-Seidel
AX  B; A  L  D  U
( L  D) X k 1  UX k  B
• Nombre d’itérations pour diviser le résidu
par 10 de l’ordre de N2/4 (Laplacien en DF
d’ordre 2)
MD "Calcul Scientifique"
105
Méthodes itératives
• Successive Overrelaxation (SOR)
AX  B;
A  L  D U

X k 1  X k   ( L  D) 1 ( L  D  U ) X k  B
0  2
optimal 

2
2
1  1   Jacobi
• Nombre d’itérations pour diviser le résidu
par 10 de l’ordre de N/3 (Laplacien en DF
d’ordre 2).
MD "Calcul Scientifique"
106
Méthodes itératives
• A définie positive:
résoudre AX  B revientà
1
minimiser f ( X )   AX , X    B, X 
2
alors f  AX-B  0
MD "Calcul Scientifique"
107
Minimisation
• Plus grande pente:
X k 1  X k  f
 fixé " petit" ou optimal
MD "Calcul Scientifique"
108
Exemple
• On cherche à résoudre
AX  B
1 1  x  1
1 2  y   1

   
• C’est-à-dire à minimiser
1 t
f ( x, y )  X AX  X t B
2
1 1  x 
1
f ( x, y )  x y 
 x



2
1 2  y 
MD "Calcul Scientifique"
1
y  
1
109
Exemple: Plus grande pente
• C’est une forme quadratique
1
0 solution
 
 2
X0   
1 
2
r0   
1 
MD "Calcul Scientifique"
110
Exemple: Plus grande pente
• C’est une forme quadratique
 1 
X1  

3
/
2


 0. 5
r1   
  1
MD "Calcul Scientifique"
111
Exemple: Plus grande pente
• C’est une forme quadratique
  1 / 2
X2  

1
/
2


 1 
r2   
1 / 2
MD "Calcul Scientifique"
112
Exemple: Plus grande pente
• C’est une forme quadratique
 0 
X3  

3
/
4


 1/ 4 
r3  

 1 / 2
MD "Calcul Scientifique"
113
Exemple: Plus grande pente
• C’est une forme quadratique
1 / 4
X4   
1 / 4
1 / 5
r4   
1 / 4
PAS FINI !
MD "Calcul Scientifique"
114
Gradient conjugué
• A symétrique définie positive
• Minimiser sur l’ensemble des directions
explorées0
r0  B  AX
d 0  r0
P our k  0,1,...:
t
k k
rr
k  t
d k Adk
rkt1rk 1
 k 1  t
rk rk
xk 1  xk   k d k
rk 1  rk   k Adk
d k 1  rk 1   k 1d k
MD "Calcul Scientifique"
115
Exemple: Gradient conjugué
• C’est une forme quadratique
1
0 solution
 
 2
X0   
1 
2
r0   
1 
MD "Calcul Scientifique"
116
Exemple: Gradient conjugué
• C’est une forme quadratique
 1 
X1  

3
/
2


 0. 5
r1   
  1
 1 
d1  


3
/
4


MD "Calcul Scientifique"
117
Exemple: Gradient conjugué
• C’est une forme quadratique
1
X2   
0 
FINI !
MD "Calcul Scientifique"
118
Gradient conjugué
• On montre que:
 k  m, r d m  0 
t
k
Ad  d
t
k
m
 0  rkt rm  0
• Puisqu’il existe au plus N vecteurs
orthogonaux dans un EV de dimension N,
l’algorithme converge au plus en N
itérations
• Ne nécessite que des produits matricevecteur ou produits scalaires
MD "Calcul Scientifique"
119
Minimisation
• GC plus efficace lorsque les vp de A sont
peu distribuées
• GC préconditionné: C 1 AX  C 1B
r0  B  AX 0
z0  C 1r0
d 0  z0
CA
P our k  0,1,...:
t
k k
zr
k  t
d k Adk
xk 1  xk   k d k
zk 1  C 1rk 1
z kt 1rk 1
 k 1  t
z k rk
MD "Calcul Scientifique"
rk 1  rk   k Adk
d k 1  z k 1   k d k
120
Cas non-symétrique
• AtA est toujours symétrique définie positive … Mais le
conditionnement de AtA X= AtB est (A)2
• Generalized Minimum Residual: GMRES (www.cerfacs.fr)
• Trouver une approximation de la solution de AX=B sous la
forme Xm=X0+VmY
• Vm: base orthogonale de l’espace de Krylov
Vect{r0,Ar0, A2r0 ,…,Am-1r0}, r0=B-AX0
MD "Calcul Scientifique"
121