Transcript Powerpoint - Institut de Planétologie et d`Astrophysique de Grenoble

Méthodes numériques pour
l’astrophysique
M2 « Astrophysique et Milieux Dilués »
Hervé Beust
Laboratoire d’Astrophysique de Grenoble
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
1
Méthodes numériques pour
l’astrophysique
• Techniques de base
• Estimateurs et statistique
• Modélisation de données
• Résolution numérique d’équations différentielles
• Equations aux dérivées partielles
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
2
Méthodes numériques pour
l’astrophysique
• Techniques de base
• Résolution numérique d’équations
• Intégration numérique
• Minimisation de fonctions
• Estimateurs et statistique
• Modélisation de données
• Résolution numérique d’équations différentielles
• Equations aux dérivées partielles
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
3
Résolution numérique d’équations
• Problème : résoudre numériquement une équation de type
f(x) = 0, en général par une méthode itérative.
• Le problème est très différent suivant que f est une
fonction scalaire ou vectorielle.
• En général, les méthodes classiques fonctionnent assez
bien pourvu que
– On soit sûr qu’il y ait une racine;
– On parte du voisinage de la racine.
– Ca marche mieux si on encadre la racine entre deux réels a et b et
si la fonction est monotone dans l’intervalle.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
4
Différents cas de figure...
Cas standard
septembre 2011
 Cas plus délicats 
Master 2 AMD - Méthodes
numériques pour l'astrophysique
5
Méthode de dichotomie
• On prend c=(a+b)/2 et on calcule f(c). Si f(c) est du même
signe que f(a), la racine est entre c et b; sinon entre a et c.
On recommence avec le nouvel intervalle.
Do while (abs(a-b)>…)
c = (a+b)/2
if f(c)*f(a)>0 then
a = c
else
b = c
end if
end do
• Cette méthode a l’avantage de marcher tout le temps
pourvu qu’il y ait une racine. Elle ne nécessite pas le
calcul de la dérivée.
• C’est une méthode d’ordre 1 a convergence
moyennement rapide.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
6
Autres méthodes d’ordre 1
• Méthodes de la sécante et de la fausse position
• Ces méthodes convergent en
général plus vite que la
dichotomie, mais
pas toujours…
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
7
Méthode de Newton-Raphson
• On part d’un point, on approxime la
fonction par sa tangente, on
recommence avec la racine trouvée,
etc…
f(x+h) ≈ f(x)+h f ’(x)+(h2/2) f ’’(x)+..
A xn+1 = xn - f(xn)/f ’(xn)
• Cette méthode ne nécessite pas d’encadrer la racine au préalable, mais juste
de partir d’un point. Elle fonctionne dans le complexe.
• C’est une méthode d’ordre 2 très efficace.
• MAIS, il faut être sûr d’être au voisinage de la racine. La méthode peut
échouer sinon.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
8
Méthode de Newton-Raphson : Problèmes
• Il faut être sûr d’être au voisinage d’une racine…
• Il faut aussi pouvoir calculer la dérivée. Sinon:
f’(x) ≈ (f(x+h)-f(x-h)) / 2h
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
9
Amélioration : Newton-Raphson quartique
•
On peut améliorer la convergence
si on sait calculer des dérivées
d’ordre supérieur
f ( xn )
f ' ( xn )  12 1 f ' ' ( xn )
f ( xn )
3  
f ' ( xn )  12  2 f ' ' ( xn )  16  22 f ' ' ' ( xn )
xn 1  xn   3
•
Exemple : équation de Képler
M = u - e sin u
septembre 2011
IMPLICIT NONE
REAL*8
M,
! Mean anomaly
&
E,
! Eccentricity
&
U,
! Eccentric anomaly
&
F,F1,F2,F3, ! Intermediaires
&
DELA
! Correction
REAL*8, PARAMETER :: EPSI=1d-12
LOGICAL
OK
INTEGER
I
f ( xn )
1  
f ' ( xn )
2  
SUBROUTINE KEPLER0(M,E,U)
U=M
OK = .FALSE.
I=0
DO WHILE (.NOT.OK)
F = U-E*SIN(U)-M
F1 = 1.d0-E*COS(U)
F2 = U-M-F
F3 = 1.d0-F1
DELA = -F/F1
DELA = -F/(F1+0.5d0*DELA*F2)
DELA = -F/(F1+0.5d0*DELA*F2+DELA*DELA*F3/6.d0)
U = U+DELA
I = I+1
OK = ((ABS(DELA).LT.EPSI).OR.(I.GE.500))
END DO
END
Master 2 AMD - Méthodes
numériques pour l'astrophysique
10
En multidimensions …
• En dimensions plus que 1, c’est
beaucoup plus difficile.
 f1 ( x1 ,, xn )  0

.......................... 
 f ( x , , x )  0
n
 n n
 
f ( x)  0
• La seule méthode générale, c’est
Newton-Raphson, sans garantie de
succès…




1
xn1  xn  J xn ( f ( xn ))
• Jxn = Matrice jacobienne en x
f i
( J xn ) i , j 
( xn )
x j
septembre 2011
• A chaque étape il faut
résoudre un système linéaire
• La convergence est
quadratique pouvu qu’on soit
au voisinage de la racine…
Master 2 AMD - Méthodes
numériques pour l'astrophysique
11
Intégration numérique
• Problème : calculer numériquement une intégrale de la forme

b
a
f ( x) dx ou

b
a
w( x) f ( x) dx avec
w( x)  0
• w(x) est une fonction de poids
• La technique consiste toujours à faire une combinaison linéaire
d’évaluations de la fonction f en un certain nombre de points xi , i=1..n

b
a
n
w( x ) f ( x ) dx  ai f ( xi )
i 1
• Tout repose dans le choix des coefficients ai et des points xi.
• Les formules convergent vers l’intégrale lorsque n ∞
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
12
Formules à points régulièrement espacés
x1  a, xn  b,
w( x )  1, h  f ( xi 1 )  f ( xi )
• Formule des trapèzes :

b
a
f ( x ) dx  h
1
2
1
f ( x1 )  f ( x2 )  f ( x3 )    f ( xn 1 )  f ( xn )  O  2 
n 
1
2
• Formule de Simpson :

b
a
 13 f ( x1 )  43 f ( x2 )  23 f ( x3 )  43 f ( x4 ) 
1
f ( x ) dx  h 

O
 4

2
4
1


f
(
x
)

f
(
x
)

f
(
x
)
n 
n 2
n 1
n
3
3
3


septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
13
Quadratures de Gauss – polynômes orthogonaux
•
•
•
•
Idée : ne plus imposer que les xi soient régulièrement espacés.
Choisir les xi et les ai de manière à ce que la formule d’approximation soit
exacte pour des polynômes de degré le plus élevé possible.
C’est possible jusqu’à des polynômes de degré 2n-1 : Pour un choix donné de a,
b, w(x), et n il existe une unique jeu de xi et les ai telle que l’approximation soit
exacte pour toute fonction f polynôme de degré ≤ 2n-1.
La théorie est liée à celle des polynômes orthogonaux. On définit
n
Pn   ( X  xi )
i 1
•
La suite des polynômes Pn lorsque n varie est orthogonale au sens du produit
scalaire suivant
b
f g   w( x) f ( x) g ( x) dx
a
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
14
Quadratures de Gauss – polynômes orthogonaux
• Les Pn obéissent à la récurrence suivante :
P0 ( x )  1
P1 ( x )  x 
xP0 P0
P0 P0

xPn Pn 
Pn Pn
Pn 1 ( x )   x 
Pn 1 ( x )
 Pn ( x ) 
P
P
P
P
n
n 
n 1
n 1

• Les xi se calculent comme les racines de Pn et les ai vérifient :
Pn 1 Pn 1
ai 
Pn 1 ( xi ) Pn( xi )
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
15
Quadratures de Gauss – polynômes orthogonaux
Cas particuliers :
(a,b)
w(x)
Nom des polynômes
Récurrence
(-1, 1)
1
Polynômes de Legendre
(n+1) Pn+1(x) = (2n+1) xPn(x) - nPn-1(x)
(-1,1)
(1-x²)-1/2
Polynômes de Tchébychev
Pn+1(x) = 2xPn(x) - Pn-1(x)
(0,+∞)
xc e-x
Polynômes de Laguerre
(n+1) Pn+1(x) = (-x+2n+c+1) xPn(x) (n+c) Pn-1(x)
(-∞,+∞)
e-x²
Polynômes de Hermite
Pn+1(x) = 2xPn(x) - 2nPn-1(x)
Dans le cas de Gauss-Legendre, on a en outre
Et pour tout intervalle (a,b) on a

b
a
septembre 2011
ai 
2
2
(1  xi2 )Pn( xi )
ba 1 ab ba 
f ( x ) dx 
f

y  dy
2 1  2
2

Master 2 AMD - Méthodes
numériques pour l'astrophysique
16
Exemple résolu : Gauss-Legendre
SUBROUTINE GAULEG(X1,X2,X,W,N)
IMPLICIT NONE
INTEGER*4
N,M,
! Ordres
&
I,J
! Indices
REAL*8
X(N),
! Tableau de racines
&
W(N),
! Tableau de poids
&
Z,
! Racine
&
X1,X2,XM,XL, ! Bornes
&
Z1,P1,P2,P3, ! Intermediaires
&
PP
! Derivee
LOGICAL
OK
! Test d'arret
REAL*8, PARAMETER :: EPS=3.d-14,
&
PI=3.14159265358979323846d0
M = (N+1)/2
XM = 0.5d0*(X2+X1)
XL = 0.5d0*(X2-X1)
! Conversion d’intervalle
DO I = 1,M
Z = COS(PI*(I-0.25d0)/(N+0.5d0)) ! Guess initial
OK = .FALSE.
DO WHILE(.NOT.OK)
P1 = 1.0d0
P2 = 0.0d0
DO J = 1,N
septembre 2011
P3 = P2
P2 = P1
P1 = ((2.0d0*J-1.0d0)*Z*P2-(J-1.0d0)*P3)/J
END DO
! Calcul de Pn(Z) par récurrence
PP = N*(Z*P1-P2)/(Z*Z-1.0d0)
Z1 = Z
Z = Z1-P1/PP ! Newton-Raphson
OK = (ABS(Z-Z1).LT.EPS)
END DO
X(I) = XM-XL*Z
X(N+1-I) = XM+XL*Z
W(I) = 2.d0*XL/((1.0d0-Z*Z)*PP*PP)
W(N+1-I) = W(I)
END DO
END
Master 2 AMD - Méthodes
numériques pour l'astrophysique
17
Intégration multidimensionnelle
• C’est en général beaucoup
plus difficile.
• Technique de base :
intégrales emboîtées
• 2 difficultés :
– Nombre de points nd
– Ca ne marche que si le
domaine à intégrer n’est pas
trop compliqué
b
 d ( x ) f ( x, y ) dy  dx
a c( x )

• Dans les autres cas :
Intégration Monte-Carlo
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
18
Intégration Monte-Carlo
• On veut intégrer une fonction f sur un domaine W. La moyenne de f
sur W vaut par définition
1
f W

W
f dV
• On tire au hasard un nombre n de points dans l’ensemble et on
approxime la moyenne par
1 n
f 
f (x )

n
i
i 1
• L’erreur moyenne sur l’intégrale est
a convergence en n-1/2, lente
W
f
2
 f
2
n
• Si W n’est pas calculable directement, on inclut le volume W dans un
volume V plus grand, facilement calculable, et on prolonge f par 0 en
dehors de W. On intègre ensuite sur V. Mais il faut choisir un volume
V proche de W sinon, la convergence est encore plus lente.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
19
Minimisation de fonctions
•
•
•
•
Problème : Trouver numériquement un point minimum
(ou maximum) d’une fonction f donnée
Un minimum peut être local ou global : les méthodes
itératives garantissent en général la convergence vers
un minimum local.
Le problème est très différent suivant que la fonction
est mono- ou multidimensionnelle
Il y a deux familles de méthodes : les méthodes avec
gradient et celles sans gradient.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
20
Méthodes en dimension 1
• On veut minimiser une fonction f d’une variable réelle.
• Avant de converger, il faut encadrer un minimum : trouver 3 points
a<b<c, tels que f(b)<f(a) et f(b)<f(c).
• Pour encadrer, on déplace le triplet dans le sens de la descente jusqu’à
trouver un encadrement.
• Première méthode simple une
fois qu’on a un encadrement:
On tire un réel d entre a et c,
on évalue f(d), et on cherche un
nouvel intervalle d’encadrement…
puis on recommence…
• Ca marche toujours, mais ce
n’est pas très rapide.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
21
Méthode : interpolation par paraboles inverses
• Idée : près d’un minimum, la fonction s’approxime bien
par une parabole.
• Quand on a un encadrement (a,b,c), on calcule la fonction
de degré 2 (parabole) qui passe par ces trois points et on
en cherche le minimum d. On a
1 ( b  a ) 2  f ( b)  f ( c )   ( b  c ) 2  f ( b)  f ( a ) 
d b
2 (b  a ) f (b)  f (c)  (b  c) f (b)  f (a )
• Le nouvel encadrement est (a,d,b) ou (b,d,c).
• La méthode de Brent combine les deux méthodes
précédentes.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
22
Méthodes multidimensionnelles
• On cherche à minimiser une fonction f(x1,…,xn) = f(x)
• La plupart des méthodes procèdent à des minimisations
successives le long de plusieurs directions.
• Minimisation le long d’une direction: étant donné deux
vecteurs x et n, trouver le réel l qui minimise f(x+ln).
a On utilise une méthode de dimension 1
• A partir de x+ln, on recommence à minimiser dans une
autre direction n’.
• Problème : En minimisant le long de n’, il ne faut pas
détruire la minimisation le long de n.
a Les directions n et n’ doivent être conjuguées.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
23
Directions conguguées
• On se place en un point P
 

f
1
2 f
f ( P  x )  f ( P)  
xi  
xi x j  

2 i , j xi x j 
i xi P
P
avec
   1

 f ( P)  b  x  x  A  x  
2
 
2 f
b  f  et A ij 
( matrice Hessienne )
P
xi x j 
P
• Le changement de gradient dans un déplacement infinitésimal vaut
 


 f  A  x 
• Une fois qu’on a minimisé le long de u, si on recommence le long de v,
il faut que le gradient reste perpendiculaire à u.
 
 


u   f  0  u  A  v  0
• Les directions u et v sont alors dites conjuguées.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
24
Méthode de Powell
•
•
•
Idée : faire des minimisations successives le long de directions conjuguées.
Problème : Comment trouver des directions conjuguées en chaque point ?
La méthode de Powell minimise le long de directions qui tendent à être
conjuguées au fil des itérations : On a à chaque instant n directions ui, i=1..n.
Au départ, on prend ui=ei (vecteurs de base). On répète ensuite la séquence
suivante:
– On part d’une position de départ x0.
– On fait n minimisation successives le long des n directions. On appelle x1,… xn les
points trouvés.
– On remplace u0 par u1, puis u1 par u2, …, un-1 par un.
– On remplace un par xn-x0
– On minimise dans la direction de un à partir de xn et on remplace x0 par le point
trouvé.
•
•
•
La méthode a un défaut : elle tend à produire des directions ui linéairement
dépendantes
Solution : Réinitialiser les directions (ui=ei) toutes les n itérations.
Remarque : Cette méthode est une méthode sans gradient.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
25
Méthodes avec gradient
• Quand on sait calculer le gradient f en tout point…
• Méthode de descente maximale : On part de x, on minimise dans la
direction de -f , et on recommence.
• Cette méthode ne donne pas de bons résultats car les directions
successives ne sont pas conjuguées.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
26
Méthode de Fletcher-Reeves-Polak-Ribiere
• On part d’un point
 x0, on fixe une direction u0, et on

calcule g0  f 
x0

 

• A chaque étape on dispose de xi , ui et gi  f 
xi
• On minimise le long de ui. On appelle xi+1 le nouveau

point trouvé. On calcule g  f
i 1
• On calcule ui+1 comme



ui 1  gi   i ui

xi 1

 
g  g   g
avec  i  i 1  i i 1
gi  gi
• On recommence.
• Cette méthode garantit que les directions ui sont
conjuguées. Elle est recommandée quand on n’a aucune
idée de la matrice Hessienne.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
27
Méthodes numériques pour
l’astrophysique
• Techniques de base
• Estimateurs et statistique
• Propriétés et exemples
• Convergence
• Recherche d’estimateurs
• Estimation par ajustement
• Modélisation de données
• Résolution numérique d’équations différentielles
• Equations aux dérivées partielles
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
28
Estimateurs et statistiques
•
•
•
•
•
Une série de mesures doit être considérée comme une séries de réalisations
d’une ou plusieurs variables aléatoires
Le phénomène aléatoire dépend d’un paramètre physique θ que l’on
cherche à estimer.
Un estimateur est une fonction T arbitraire d’un échantillon de données
(X1,…,Xn) censée représenter (estimer) le paramètre inconnu θ
En général, on a plutôt une suite d’estimateurs (Tn) de θ en fonction de la
taille n de l’échantillon
Exemple : On cherche à estimer la moyenne (espérance) d’une variable
aléatoire. Si on fait n tirages successifs, un estimateur de la moyenne est
Xn 
•
•
X1    X n
n
On dit que l’estimateur est convergent si pour tout ε>0, on a
lim PTn       0
n
La loi faible des grands nombre montre que Xn est un estimateur convergent
de l’espérance.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
29
Propriétés des estimateurs
• Si (Tn) est un estimateur convergent du paramètre θ, si φ est une
fonction de IR dans IR, continue en θ, alors φ(Tn) est un estimateur
convergent de φ(θ).
• Exemple : X suit la loi uniforme sur [0,θ]. Xn est un estimateur
convergent de la moyenne = θ/2  Tn=2Xn est un estimateur
convergent de θ.
n
• On pose Y=ln X. On a alors E(Y)=ln θ – 1 (calcul). Donc n1  ln X i
i 1
est un estimateur convergent de ln θ – 1 . Et ensuite
•
1 n

Tn  exp  n  ln X i  1
 i 1

est un autre estimateur convergent de θ
On peut en créer beaucoup d’autres. Comment sélectionner le
meilleur ?
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
30
Propriétés des estimateurs
• Variance de Tn : V[Tn] = E[(Tn-E(Tn))2]
• Erreur quadratique de Tn par rapport à θ:
EQ(Tn,θ) = E[(Tn-θ)2]
• Biais de Tn par rapport à θ : B(Tn,θ) = E[Tn-θ].
• Propriété 1 : si EQ(Tn,θ))0 quand n∞, alors l’estimateur
Tn est convergent.
• Un estimateur est meilleur qu’un autre si son erreur
quadratique est inférieure.
• Avec un estimateur biaisé, E[Tn] ≠ θ (décalage). On préfère
en général des estimateurs sans biais (B(Tn,θ) = 0).
• Un estimateur est asymptotiquement sans biais si B(Tn,θ)0
quand n∞.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
31
Exemples
• On reprend X, loi uniforme entre 0 et θ. On a Tn et Tn’,
estimateurs de θ.
• E(Tn) = θ  B(Tn,θ) = 0; EQ(Tn,θ) = θ2/(3n).
• B(Tn’) ~ θ/(2n) ; EQ(Tn’,θ) ~ θ2/n
a Tn est meilleur que Tn’ (non biaisé et plus petite erreur
quadratique).
• Autre estimateur de θ : Tn″ = max {X1,…,Xn}, convergent et
biaisé.
• E(Tn″) = nθ/(n+1)  B(Tn″,θ) = -θ/(n+1);
EQ(Tn,θ) = 2θ2/(n+1)(n+2) ~ 2θ2/n2.
• On construit alors Tn″′=(n+1)Tn″ /n, non biaisé…
EQ(Tn,θ) = θ2/n(n+2) ~ θ2/n2.
C’est le meilleur des 4 (non biaisé et convergence en n2)
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
32
Convergence des estimateurs
• EQ(T,θ) = E[(T-θ)2]= E[(T-E(T)+E(T)-θ)2]
= E[(T-E(T))2] + (E(T)-θ)2 + 2(E(T)-θ)×E[T-E(T)]
= Var(T) + B(T,θ)2
• Conséquence : Si un estimateur est sans biais (même
asymptotiquement), et si sa variance tend vers 0, il est
convergent.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
33
Estimateurs standards
• On a (X1,...,Xn), un échantillon de
loi inconnue.
• Estimateur de l’espérance
convergent et non biaisé
Var(Xn)=Var(X)/n
• Estimateur de la variance
Convergent, mais biaisé (variance
empirique)
E[Sn2]=(n-1)/n×Var(X)
• Estimateur non biaisé
(variance empirique non biaisée)
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique


n
Xn
i 1
Xi
n
2
X
i1 i
n
Sn2 
n
 X n2
n
Vn 
Sn2
n 1
34
Recherche d’estimateurs
• Je cherche à estimer un paramètre physique à partir de
statistiques : Comment trouver un bon estimateur ?
• Ceci nécessite de faire une hypothèse sur le processus
physique, par exemple une loi de probabilité
(modélisation de données)
• Exemple : méthode des moments. Pour k>0 E[Xk] et
E[(X-E(X))k] sont les moments de la variable aléatoire X.
– Je suppose une loi de probabilité d’une forme particulière,
dépendant de k paramètres a1,...,ak que je cherche à estimer.
– Les k premiers moments peuvent s’exprimer en fonction des
a1,...,ak  Les paramètres s’expriment en fonction des moments.
– Il suffit d’avoir des estimateurs des moments (p.ex la moyenne et
la variance) pour en déduire des estimateurs des paramètres.
– Inconvénient de la méthode : les estimateurs qu’on tire sont
souvent peu précis.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
35
Estimation par ajustement
• Idée : définir une métrique sur les lois de probabilités, et trouver la loi
de probabilité P qui minimise une « distance » par rapport à la loi
empirique Pe.
• Si (x1,…,xn) est un échantillon observé, (c1,…,ck) l’ensemble des
valeurs prises par les xi’s, on définit Pe(ci)=n(ci)/n, où n(ci) est le
nombre de fois où ci est réalisé (fréquence empirique).
• Distance du χ2 :
2
k
D 2 ( P, Pe )   2  
i 1
( P(ci )  Pe (ci ))
P(ci )
• Distance de Kolmogorov-Smirnov = comparaison entre les fonctions
de répartition théoriques (F) et empiriques (Fe)
DKS ( F , Fe )  max F (ci )  Fe (ci )
i 1,,n
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
36
Méthodes numériques pour
l’astrophysique
• Techniques de base
• Estimateurs et statistique
• Modélisation de données
• Régression linéaire
• Modèles linéaires de moindres carrés
• Modèles non-linéaires : méthode de Levenberg-Marquardt
• Recuit simulé et algorithmes génétiques
• Résolution numérique d’équations différentielles
• Equations aux dérivées partielles
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
37
Modélisation de données par régression
•
•
•
•
•
On mesure une série de données (x1,y1)… (xn,yn). On cherche à trouver une relation
physique y=F(x).
Dans la pratique on se donne un modèle dépendant d’un nombre k de paramètres
α1,…,αk:
y=f(α1,…,αk; x).
On va chercher les paramètres α1,…,αk qui collent au mieux avec les observations
(xi,yi) .
Chaque mesure yi vient avec son erreur standard σi2. Yi=f(α1,…,αk; xi)+Ei , où Ei est
une variable aléatoire de variance σi2.
On va chercher les paramètres qui minimisent
 y  f (a1 ,, a k ; xi ) 

 2    i
i
i 1 

n
•
2
Toute la difficulté consiste à
1.
2.
Trouver les paramètres réalisant le minimum
Estimer l’incertitude sur les paramètres.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
38
Régression linéaire
• But : modéliser des données
par une relation linéaire
y = ax+b
• Méthode : Minimiser le χ2
pour trouver a et b
 y  (axi  b) 

 2    i
i
i 1 

n
2
n
  2
xi ( yi  axi  b)
SS xy  S x S y



2
0

a

 a
2

aS xx  bS x  S xy
i
i 1



 2
n
S S  S x S xy
 aS x  bS  S y
   2 yi  axi  b  0
b  xx y
 b


 i2
i 1
n
n
n
n
1
xi
yi
xi2
xy
S   2 , S x   2 , S y   2 , S xx   2 , S xy   i 2 i ,   SS xx  S x2
i 1 σ i
i 1 σ i
i 1 σ i
i 1 σ i
i 1 σ i
n
avec
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
39
Régression linéaire (2)
• Erreur sur a et b :
– Var(yi)=σi2  Varaiyiai2 σi2
(Somme de variables aléatoires indépendantes)

SS  S x S y
 a  xy



S S  S x S xy
b  xx y


a 
S

septembre 2011
n

 Sxi  S x  yi
a


 2




 i
i 1
 
n
b    S xx  S x xi  yi2


 i
i 1 
; b 
S xx

2
n

S
 Sxi  S x  1
 2 
 Var (a )   
  i


i 1 

2
n
Var (b)   S xx  S x xi  1  S xx

2



 i
i 1 
; cov(a, b)  
Sx

Master 2 AMD - Méthodes
numériques pour l'astrophysique
;  ab  
Sx
SS xx
40
Modèles linéaires de moindres carrées
• On cherche un modèle qui dépend linéairement des paramètres
y = f(α1,…,αk; x) = α1X1(x)+∙∙∙+αkXk(x)
où les X1,…,Xk sont des fonctions données à l’avance
• On va chercher à minimiser
2
k

n  y 
a X j ( xi ) 

i
j 1 j

2
  

i
i 1 


en disant ∂χ2/∂αi=0 pour i=1…k
• On introduit la matrice A, n×k, et le vecteur b de longueur n tels que
Aij 
septembre 2011
X j ( xi )
i
et bi 
yi
i
Master 2 AMD - Méthodes
numériques pour l'astrophysique
41
Modèles linéaires de moindres carrés
• ∂χ2/∂αj=0 pour j=1…k revient à
k
c a
j 1
ij
j
 di
ou encore
n
avec
cij  
X i ( xl ) X j ( xl )
 l2
l 1
c  t A  A
et
n
et d i  
l 1
d  t A  b
yl X i ( xl )
 l2
…autrement dit résoudre un système linéaire ! (équations normales)
• On appelle fij = [c]-1ij . On a
 n yl X j ( xl ) 
ai   c d j   fij 

2

j 1
j 1
l
 l 1

k
 n X j ( xl ) X p ( xl ) 
 Var (ai )   f ij f ip 
  f ii
2

j 1 p 1
1
l
l 



k
k
c jp
septembre 2011
k
1
ij
i 
Master 2 AMD - Méthodes
numériques pour l'astrophysique
fii
et cov ai , a j   fij
42
Modèles non linéaires
• y = f(α1,…,αk; x) où la dépendance est quelconque.
• Le jeu de paramètres (α1,…,αk) minimisant le χ2 doit être
trouvé de manière itérative.
• C’est conceptuellement identique à un problème de
minimisation dans un espace de dimension k.
• La non-linéarité ne garantit pas l’existence d’un minimum
unique  problèmes de minima locaux  nécessité de
faire de nombreux essais !
• Quand on peut calculer le gradient de la fonction f par
rapport aux paramètres, la méthode recommandée est celle
de Levenberg-Marquardt.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
43
Méthode de Levenberg-Marquardt
•
•
On appelle α le vecteur (α1,…,αk) courant.
Si le modèle est linéaire, la fonction χ2(α) est une forme quadratique




 2 (a )    v  a  12 a  H  a
•
v est un vecteur et H une matrice k×k
Dans ce cas on sait où est le minimum : αm = H-1∙v. De manière équivalente, si
 2 
on part d’un vecteur α0 , on a


1

am  a0  H    (a0 )
•
•

Si la fonction n’est pas trop éloignée d’une forme quadratique, cette formule
peut-être une bonne formule d’itération.
Si la fonction est plus compliquée, on se contentera d’une simple descente de
 2 
gradient


am  a0  cte   (a0 )
•
•
Idée de la méthode : Alterner entre les deux formules
Problème : Il faut la matrice H = matrice Hessienne
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
44
Méthode de Levenberg-Marquardt (2)
2



y

f
(
a
;
x
)
i
 2 (a )    i

i
i 1 



n
 2
yi  f (a ; xi ) f (a ; xi )
 2
 2d j
a j


a
i 1
i
j



n

2  2
1  f (a ; xi ) f (a ; xi )
 2 f (a ; xi ) 
 2 2 
  yi  f (a ; xi ) 
  2c jl
a j a l


a

a

a

a
i 1
i 
j
l
j
l 



n

 2 
a m  H    (a0 )

1


k
 c a
l 1
jl
l
 d j j  1 k (1)
 2 

am  cte   (a0 )  a j  cte  d j
•
•
(2)
Problème : Les termes cjl dépendent des dérivées secondes de f qui ne sont pas
forcément accessibles
Souvent on laisse tomber les termes en question…


n
1  f (a ; xi ) f (a ; xi ) 
c jl   2 


a

a
i 1  i 

j
l

septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
45
Méthode de Levenberg-Marquardt (3)
•
Deux questions :
1.
2.
•
Comment décider de l’alternance entre les deux formules (1) et (2) ?
Comment fixer la constante dans la formule (2) ?
Levenberg-Marquardt :
1.
2.
cte ~ 1/cjj . En fait cte = 1/(λcjj)  (2) devient λcjjδαj = dj
(1) et (2) peuvent se combiner en une seule équation. On définit C’= C +
λ×diag(c11,…,ckk) et on remplace (1) et (2) par
k
 c a
l 1
3.
•
jl
l
 d j j  1 k (3)
(cjj  c jj (1  l ); cjl  c jl j  l )
Quand λ>>1, (3) revient à (2); quand λ<<1, (3) équivaut à (1)
Recette :
1.
2.
3.
4.
On part d’un choix initial α, et on calcule χ2(α). On prend un petit λ =0.001
On résout (3) et on calcule l’incrément a. On calcule χ2(αa.
Si χ2(αa  χ2(α), on est loin du minimum  on multiplie λ par 10 et on
revient à 2.
Si χ2(αa < χ2(α), on approche  On divise λ par 10, on remplace α par
αa, et on revient à 2.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
46
Intervalles de confiance des paramètres
• Une fois qu’on a obtenu les paramètres (α1,…,αk) qui
minimisent, on souhaite obtenir des intervalles de
dispersion (δα1,…,δαk) des paramètres à un niveau de
confiance donné.
• La méthode de Levenberg-Marquardt fournit la matrice des
covariances F = C-1. Mais celle-ci n’est utilisable que si
les erreurs de mesure sont distribuées de manière
gaussienne ! Si ce n’est pas le cas, on peut en déduire
n’importe quoi !!
• Si les erreurs ne sont pas gaussiennes (ou si on ne sait pas),
il faut avoir recours à d’autres méthodes statistiques
beaucoup plus robustes…
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
47
Avec la matrice des covariances…
• Si les erreurs sont gaussiennes, alors on a a j   f jj avec un degré
de confiance de 68% (1σ dans la loi normale).
• Attention : le jeu d’intervalles de confiances δαj, j=1..n ne constitue pas
une région de confiance pour les k paramètres conjointement…
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
48
Méthodes statistiques
ou quand les erreurs ne sont pas gaussiennes…
• Simulation de jeu de données Monte-Carlo
– Idée : tirer plusieurs (N) jeux de données fictives dans les intervalles de
confiance ±σi (de manière gaussienne…), et effectuer N minimisations.
– On observe la distribution des jeux de paramètres α obtenus, on en déduit
des intervalles de confiance.
• Méthode Bootstrap
– Au lieu de tirer des jeux de données fictives, on tire au hasard n données
(avec éventuelle répétition) parmi les n de base, et on refait une
minimisation.
– On répète l’opération N fois, on obtient la distribution des α .
• Avantage de ces méthodes : Elles sont simples, robustes et ne
présupposent rien sur la distribution des erreurs de mesure.
• Désavantage : Elles sont coûteuses en temps de calcul
(il faut refaire N minimisations).
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
49
Méthode MCMC
(Monte-Carlo + Chaînes de Markov)
• Idée : ne plus chercher d’intervalle de confiance pour les paramètres,
mais échantillonner leur distribution statistique.
• C’est une amélioration de la méthode de simulation de jeu de données.
• On a un vecteur de données y. Pour un jeu de paramètres a, on connaît


p y a , probabilité d’observer y sachant a. On cherche p a y
• Par le théorème de Bayes
 
 

 
 
  p y, a  p y a  pa 
pa y  

 

p y 
p y 

• On suppose connaître pa  (prior).
p y a  pa 
  
 py a pa da
• Une chaîne de Markov = une suite de modèles ai, où chaque
ai+1 se

déduit de ai via une probabilité de transition ptr a i 1 a i . A la fin, la

distribution des ai échantillonne p a y
 
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique


50
Méthode MCMC (2)
• Définition : La chaîne de Markov est réversible si

 
 
 
pa y  ptr a i a i 1   pa i 1 y  ptr a i 1 a i 
• Théorème : Si la chaîne est réversible, apériodique et irréductible, ell

converge vers p a y
• Problème : Il n’est pas facile de trouver une probabilité de transition
 
réversible. On peut en construire une à partir d’une autre qtr a i 1 a i
qui ne l’est pas, par l’algorithme de Metropolis-Hastings
 


 
 
 
ptr a i a i 1   qtr a i 1 a i  a i 1 a i 
 
 





p
a
y
q
a
 
i 1
tr
i 1 a i 
avec  a i 1 a i   min   
,1
 
 pa i y qtr a i a i 1  
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
51
Méthode MCMC (3)
• Méthode à chaque pas :





1. Connaissant ai, on tire au hasard un ai+1 en suivant qtr a i 1 a i
 
2. On calcule 2(ai+1) et 2(ai).
p a i 1 y
2 
2
1
3. On détermine le rapport

exp


a


 
i 1
2
 
pa y 
 
ai 1 
i
4. On tire un nombre u au hasard entre 0 et 1.
 a i 1 a i  . Si u<, on accepte le pas, sinon on revient à ai.
• Souvent, on prend qtr symétrique
    2

 

v  a i  v  a i 1  
1
qtr a i 1 a i  
exp 

2
2
2

2v
v


dans une direction v. On fait tourner les directions.
5. On calcule


• On peut prendre qtr non symétrique quand on considère le modèle
uniforme dans des combinaisons de paramètres. Dans ce cas, on
multiplie par le rapport des jacobiens.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
52
Méthode MCMC (4)
• Test d’arrêt : Statistique de Gelman-Rubin.
• On fait tourner 10—20 chaînes en parallèle. On teste la
convergence des variances des différentes chaînes.
• Quand la convergence est acquise (ça peut être long !), on
laisse tourner les chaînes plus longtemps en retenant tous
les modèles calculés. Ceux-ci échantillonnent la
distribution des paramètres.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
53
Problèmes de minima locaux…
• La méthode de Levenberg-Marquardt permet de converger vers jeu de
paramètres α qui minimise localement le χ2. Si le modèle est nonlinéaire, il peut y avoir d’autres minima  Comment savoir si on est
dans le bon ?
• Technique 1 : Je recommence N fois la minimisation en partant de
points de départ différents
– Efficace s’il n’y a pas trop de minima locaux et si la dimension de
l’espace n’est pas trop élevée ( ≤ 5-6)
• Technique 1 bis : Technique 1 + Recuit simulé
– Permet de relier des minima locaux voisins
• Technique 2 : Algorithmes génétiques
– La seule technique permettant d’explorer correctement un espace des
paramètres de dimension élevée.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
54
Recuit simulé
• A employer comme tel ou (mieux) à combiner
avec Levenberg-Marquardt;
• On introduit une fonction « température » T(n)
décroissante (souvent linéairement) en fonction du
nombre d’itération. La température initiale T0 est
élevée.
• A chaque itération, on teste un saut aléatoire de la
solution, qui induit une variation de χ2 : Δχ2.
– Si Δχ2<0, on applique le saut;
– Si Δχ2>0, on l’applique avec la probabilité e-Δχ2/T(n)
• Cette technique peut permettre de sortir d’un
minimum local.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
55
Algorithmes génétiques
•
•
Au lieu de partir d’un seul point dans l’espace des paramètres et de
minimiser, je prends N points de départ (~100) différents
(Génération 0 d’individus).
Je classe les individus par χ2, et je construis ma génération suivante
comme suit :
1.
2.
3.
4.
Je garde ceux qui ont les meilleurs χ2 ;
J’en construis d’autres en « croisant les caractères » d’individus «
parents » pris dans la génération précédente, en choisissant de
préférence des parents ayant de bons χ2
J’en construis quelques autres en introduisant des « mutations » sur des
individus de la génération précédente.
Eventuellement, j’applique quelques pas de minimisation (LevenbergMarquardt…) sur ma population et j’obtiens ma génération suivante.
Puis retour en 1…
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
56
Méthodes numériques pour
l’astrophysique
• Techniques de base
• Estimateurs et statistique
• Modélisation de données
• Résolution numérique d’équations différentielles
• Méthodes de type Runge-Kutta
• Contrôle du pas
• Méthode de Bulirsch et Stoer
• Equations mal conditionnées (Stiff)
• Codes N corps
• Conditions au limites en plusieurs points : relaxation
• Equations aux dérivées partielles
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
57
Résolution numérique d’équations
différentielles
• Un type de problème qui se pose couramment en
astrophysique (intégration de systèmes
dynamiques, calculs MHD, etc…)
• Deux grands types de problèmes … avec des
méthodes d’approche très différentes
– Les équations différentielles ordinaires (ODEs)
– Les équations aux dérivées partielles (PDEs)
• Dans tous les cas, la nature des conditions aux
limites conditionne la résolution.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
58
Equations différentielles ordinaires
• C’est résoudre le problème
• Remarques :
dy
 f ( x, y ) avec
dx
y ( x0 )  y0
– y n’est pas nécessairement scalaire. Ce peut être un vecteur  système
différentiel
– Les équations d’ordre supérieur se ramènent à ce schéma :
 dy1
 dx  y2
d y
dy
dy

f
(
x
,
y
,
)

(
y

y
,
y

) 

1
2
2
d
y
dx
dx
dx
 2  f ( x, y1 , y2 )
 dx
2


dY
 F ( x, Y )
dx
– avec une condition aux limites de la forme y(x0)=y0, dy/dx(x0)=y′0
(problème de type conditions initiales)
– Si les conditions aux limites sont de la forme y(xa)=ya et y(xb)=yb, le
problème est de nature différente (et la méthode d’approche aussi…)
(problème de conditions aux limites en deux points distincts)
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
59
Méthodes de type Runge-Kutta
• Pour un problème de type « conditions initiales ».On part de y(x0), on
en tire y(x0+h≡x1), puis y(x0+2h≡x2), etc…
h est le pas (de temps), supposé petit.
Les diverses méthodes sont d’autant plus précises que h est petit.
• La méthode la plus simple : Euler (explicite)
dy
( xn )  yn  hf ( xn , yn )
dx
Méthode d’ordre 1 : yn+1=y(xn+1)+O(h2)
peu précise (+asymétrique, peu stable)
Variante : Euler implicite, plus stable, mais non linéaire
yn 1  yn  h 
•
•
yn1  yn  hf ( xn1, yn1 )
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
60
Méthodes de type Runge-Kutta (2)
• Runge-Kutta d’ordre 2 (ou méthode du point
médian) yn+1=y(xn+1)+O(h3)
k1  hf ( xn , yn )

1
1
k2  hf ( xn  2 h, yn  2 k1 )
 yn 1  yn  k2
• Runge-Kutta d’ordre 4: yn+1=y(xn+1)+O(h5)
k1  hf ( xn , yn )
k  hf ( x  1 h, y  1 k )
n
n
2
2 1
 2
1
1
k3  hf ( xn  2 h, yn  2 k2 )
k  hf ( x  h, y  k )
n
n
3
 4
 yn 1  yn  16 k1  13 k2  13 k3  16 k4
• Remarque : Si le système est vectoriel les
ki’s sont des vecteurs…
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
61
Contrôle du pas
•
•
Comment s’assurer de la qualité de l’intégration ? Comment doit-on adapter le
pas d’intégration au fil du calcul ?
Première méthode : doublement de pas.
On calcule y(x+h) de deux façons.
– En prenant un pas h → y1
– En prenant deux pas h/2 → y2
•
•
En comparant y1 et y2, on teste si l’intégration est bonne. On pose
Δ=y2-y1  h5. Si l’erreur qu’on souhaiterait avoir est Δ0, alors il faudrait choisir
1/ 5
0
h  h 

Stratégie :
1/ 5
– Si Δ>Δ0, on recommence en diminuant le pas. On prend
– Si Δ<Δ0, on peut augmenter le pas en posant
septembre 2011
h  Sh 

h  Sh  0

Master 2 AMD - Méthodes
numériques pour l'astrophysique
0

1/ 4
( S  0.9)
62
Contrôle du pas (2)
• Que doit valoir Δ0 ? On doit avoir au moins Δ0  |y|, mais il faut se
méfier des situations où y passe par 0. Dans la pratique, on utilise
dy  
dy


0     y  h  ou min  yi  h i
dx  
dx





• Si un pas est bon, quel valeur choisir ? y1 ou y2 ? Une combinaison…
 y ( x  h )  y1  c  h 5  O (h 6 )

5
6
6
1
 y ( x  h )   151 y1  16

h
15 y2  O ( h )  y2  15   O ( h )
6
 y ( x  h )  y2  2c   2   O (h )

• Deuxième méthode : Fehlberg-Cash-Karp. Au lieu d’utiliser une seule
formule de Runge-Kutta et de doubler le pas, on utilise deux formules
différentes et on compare les résultats. Dans la pratique, on utilise une
formule d’ordre 5 et une d’ordre 4.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
63
Contrôle du pas (3) : Felhberg-Cash-carp
•
La formule d’ordre 5 est de la forme
 k1  hf ( xn , yn )
 k  hf ( x  a h, y  b k )
n
2
n
21 1
 2


 k  hf ( x  a h, y  b k  b k  b k  b k  b k )
n
6
61 1
62 2
63 3
64 4
65 5
 6
 yn 1  yn  c1k1  c2k2  c3k3  c4k4  c5k5  c6k6  y ( xn 1 )  O (h 6 )
•
La formule d’ordre 4 utilise les mêmes coefficients ai et bij, mais des coefficients c′i
yn 1  yn  c1k1  c2 k2  c3k3  c4 k4  c5k5  c6 k6  y ( xn 1 )  O(h5 )
différents
j
i
1
2
3
ai
4
5
bij
1
2
1/5
1/5
3
3/10
3/40
9/40
4
3/5
3/10
-9/10
6/5
5
1
-11/54
5/2
-70/27
35/27
6
7/8
1631/55296
175/512
575/13824
44275/110592
septembre 2011
253/4096
Master 2 AMD - Méthodes
numériques pour l'astrophysique
ci
c′i
37/378
2825/27648
0
0
250/621
18575/48384
125/594
13525/55296
0
277/14336
512/1771
1/4
64
Contrôle du pas (4) : Felhberg-Cash-carp
• L’estimation de l’erreur se fait au moyen de
6
  yn 1  yn 1   (ci  ci)ki  h 5
i 1
• Δh5, donc la correction du pas en fonction de l’erreur se
fait de la même façon que pour la méthode de doublement
de pas.
• Si le pas est accepté, on prend yn+1 comme valeur (erreur
O(h6)).
• Les méthodes de ce type donnent en général de meilleurs
résultats que les méthodes de doublement de pas.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
65
Méthode de Bulirsch et Stoer
•
•
Une méthode en général plus puissante
que Runge-Kutta, mais si la solution
est bien régulière…
Partant de x0 , on va calculer y(x0+H)
(H grand). On prend une méthode de
Runge-Kutta
– Si on prend n1 pas (h1=H/n1), on
obtient une estimation yest(x0+H,H/n1);
– Si on prend maintenant n2 pas (n2>n1),
on obtient une nouvelle (meilleure)
estimation yest(x0+H,H/n2)
•
Idée de Bulirsh-Stoer : Considérer
yest(x0+H,h) comme une fonction de h
(ou n=H/h) et extrapoler cette fonction
à h→0 (ou n→∞).
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
66
Méthode du point médian modifiée
(ou différences centrées)
• C’est la méthode de Runge-Kutta qu’on utilise pour BulirshStoer. On veut aller de x0 à x0+H en n pas de h=H/n.
 z0  y ( x0 )
 z  z  hf ( x , z )
 1
0
0 0

 zm1  zm1  2hf ( x0  mh, zm ) pour m  1,2,, n - 1
 y ( x0  h )  yn  12 zn  zn 1  hf ( x0  H , zn )
• C’est une méthode d’ordre 2 (moins précise que RK4). Mais
avantage pour Bulirsch-Stoer : L’erreur ne contient que des

puissances paires de h :
y ( x0  H )  yn   ci h 2i
i 1
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
67
Extrapoler ?
• On a des valeurs yest(x0+H,h)≡g(h) pour différentes valeurs de h :
h1,…,hk Comment extrapoler à h=0 ?
• Extrapolation polynômiale : Je calcule l’unique polynôme de degré k1 (interpolation de Lagrange) qui passe par les k points de mesure
k
(X  h )
k
P( X )   g (hi ) 
i 1
j
j 1
j i
k
 (h  h )
j 1
j i
i
j
que l’on calcule de la manière suivante, et on
évalue en X=0.
 P1  g (h1 ), P2  g (h2 ), , Pk  g (hk )
( P , P )  P , ( P , P )  P ,, ( P , P )  P
12
2
3
23
k 1
k
( k 1) k
 1 2
( P12 , P23 )  P123, ( P23 , P34 )  P234 ,, ( P( k 2 )( k 1) , P( k 1) k )  P( k 2 )( k 1) k
.................

( P123( k 1) , P234k )  P1k  P
avec la récurrence Pi ( i 1)( i  m ) 
septembre 2011
( X  hi  m ) Pi ( i 1)( i m1)  ( hi  X ) P( i 1)( i 2 )( i  m )
hi  hi m
Master 2 AMD - Méthodes
numériques pour l'astrophysique
68
Extrapoler (2) ?
• Extrapolation rationnelle : On utilise des fractions rationnelles au lieu
de polynômes. Les résultats sont en général meilleurs
 R  0; R1  g ( h1 ), R2  g ( h2 ), , Rk  g ( hk )
( R , R , R )  R , ( R , R , R )  R ,, ( R , R , R )  R
12
2
*
3
23
k 1
*
k
( k 1) k
 1 * 2
( R12 , R2 , R23 )  R123, ( R23 , R3 , R34 )  R234 ,, ( R( k 2 )( k 1) , Rk 1 , R( k 1) k )  R( k 2 )( k 1) k

( R123, R23 , R234 )  R1234, ( R234 , R34 , R345 )  R2345,, ( R( k 3)( k 2 )( k 1) , R( k 2 )( k 1) , R( k 2 )( k 1) k )  R( k 3)( k 2 )( k 1) k
...................

( R123( k 1) , R23( k 1) , R234k )  R1k  R
avec la récurrence Ri ( i 1)( i  m )  R( i 1)( i  m ) 
R( i 1)( i  m )  Ri( i  m 1)
R
 Ri( i  m 1) 
 X  hi  

 1  ( i 1)( i  m )
 1
X

h
R

R
i m  
( i 1)( i  m )
( i 1)( i  m 1) 

• On utilise en général la séquence n=2,4,6,8,…,16 (h=H/n) et on
s’arrête lorsque la correction est suffisamment petite.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
69
Equations mal conditionnées (stiff)
•
•
Des situations où les méthodes traditionnelles fonctionnent mal…
Exemple 1 : équation (y′=-cy, y(0)=1) avec c>0 (solution y(x)=e-cx), que je
choisis de résoudre par la méthode d’Euler
( y  cy )  yn1  yn  hcyn  (1  hc) yn
•
qui diverge clairement pour h>2/c (|1-hc|>1)  il faut choisir h<2/c sinon
on s’éloigne de la solution.
Exemple 2 : Système différentiel (Y′=-C.Y, Y(0)=Y0), où C est une matrice
symétrique définie positive (solution Y(x) = exp(-Cx).Y0)
( Y  C  Y)  Yn1  (I  hC)  Yn
qui diverge dès qu’une des valeurs propres de (I-hC) sort de [-1..1] ,
c’est-à-dire qu’on doit avoir h<2/λmax (valeurs propres de C)
•
Exemple 3 : (u’=998u+1998v, v’=-999u-1999v, u(0)=1, v(0)=0)
Solution (u(x)=2e-x-e-1000x, v(x)=-e-x+e-1000x)  il faut h<1/1000 , même si
dans la solution, le terme en exp(-1000x) est très vite négligeable
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
70
Equations mal conditionnées (stiff)
• Exemple 4 : équation (y″=y,y(0)=1,y′(0)=-1). Solution
y(x)=e-x. Les méthodes numériques finissent toutes par
diverger, car on introduit des éléments de la solution
parasite y(x)=ex…
• Que faire ?
– Eviter d’avoir des valeurs propres très différentes 
dédimensionnement
– Faire très attention à l’adaptation du pas
– Employer des méthodes implicites
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
71
Méthodes implicites
• Exemple 1
y
( y  cy )  yn 1  yn  hcyn 1  yn 1  n
toujours stable !
1  hc
• Exemple 2
( Y  C  Y)  Yn 1  (I  hC) 1  Yn
toujours stable aussi !
• De manière générale, les méthodes implicites sont plus
stables. Mais
– Les équations à résoudre à chaque pas sont souvent non-linéaires
yn1  yn  f ( xn , yn1 )
– Il est difficile de définir des schémas implicites d’ordre supérieur
 Méthodes semi-implicites de Rosenbrock
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
72
Méthodes semi-implicites

y n 1  y n  hf ( xn , y n 1 )  y n  h f ( y n )  J y n  ( y n 1  y n )
• On approxime…

 y n 1  y n  h I  hJ y n

1

 f (yn )
• J = matrice jacobienne. On doit résoudre un système linéaire à chaque
pas…
• Méthodes de Rosenbrock : généralisation à des schémas d’ordre
s
supérieur
y n 1  y n   ci k i
i 1
I  
ij
hJ y n

i 1
i 1


 k i  hf  y n  aijk j   hJ y n    ijk j
j 1
j 1


pour
i  1,, s
• Les γ, γij, αij, ci et l’entier s (ordre de la méthode) sont des
caractéristiques de la méthode
• A chaque pas, il faut inverser s matrices.
• On ajuste le pas en comparant avec une autre formule avec des c′i.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
73
Codes N corps
• Le problème (gravitationnel) des N corps est un
exemple de système différentiel d’ordre 6N qui
peut s’intégrer par les méthodes classiques.
• Il y a deux types de difficultés spécifiques à ce
problème
– Lorsque N est modéré (mécanique céleste) on a souvent
besoin d’intégrer pendant longtemps avec une grande
stabilité  méthodes symplectiques.
– Lorsque N est grand, la difficulté réside dans le calcul
des N(N-1)/2 termes de forces  codes autogravitants
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
74
Intégration symplectique
• La particularité des systèmes N corps est d’être des systèmes
hamiltoniens.
• On considère un système dynamique régi par un Hamiltonien
conservatif H(xi,pi)
:
dxi H

dt pi
,
dpi
H

dt
xi
• Pour toute quantité q(xi,pi)
 q H q H 
dq
  F ( q)
  

dt i  xi pi pi xi 
• La solution donnant q(t) à partir de q(t-t) est


t 2F 2
q(t )  exp(tF )q(t  t )  1  tF 
 q(t  t )
2


septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
75
Intégration symplectique (2)
• Une intégration numérique, c’est trouver q(t) connaissant q(t-t
• Les méthodes classiques ne garantissent pas H=cte
(Runge-Kutta, Bulirsh & Stoer …)
• Une méthode
symplectique
conserve exactement
H ou un autre
Hamiltonien voisin
de H. Du coup,
l’erreur faite sur H
est bornée
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
76
Intégration symplectique (3)
• Hypothèses de base de l’intégration symplectique :
H = HA+HB (F = A+B) où on sait intégrer HA et HB
(on sait calculer exp(τA) et exp(τB))
• Souvent, on suppose en plus HB<<HA (HB/HA ~ ε <<1)
• On a alors
2

exp(tA). exp(tB)  exp t  A  B   Ot

a En intégrant d’abord HA pendant t puis HB pendant t,
on réalise un intégrateur symplectique d’ordre 1 : On
résout exactement un Hamiltonien
H integ  H  O (t )
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
77
Intégration symplectique (4)
•
On a aussi
 
 
exp t B . exptA. exp t B  expt ABOt 3
2
2
a On obtient un intégrateur symplectique d’ordre 2 en
intégrant 1) HB pendant t/2,
2 HA pendant t, 3 HB pendant t/2.
Dans ce cas,
2
H integ  H  Ot
•
•

Il y a aussi des méthodes d’ordre 4, 6, 8…
En fait, compte tenu de HB/HA~ε on a même
H integ  H  Ot 2 
a On peut intégrer avec un grand pas de temps
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
78
Intégrateurs symplectiques pour codes N corps
• Exemple le plus simple : la méthode T+U
– H=Hamiltonien du problème à N corps = E. cinétique + E. potentielle
= T+U
– Variable conjuguées cartésiennes : xi=x,y,z.. pi=mvx,mvy,mvz…
– T = T(pi) et U=U(xi)  On sait intégrer séparément T et U
– Méthode symplectique avec HA=T, HB=U. Mais on n’a pas
HB<<HA
• Deuxième exemple: la méthode MVS (Mixed Variable
Symplectic)
– Hypothèse : (mi<<m0 pour tout i = 1,…,n-1) (système planétaire)
– HA = Somme des Hamiltoniens Képlériens (on sait intégrer…)
– HB = Le reste = Les perturbations mutuelles  dépend que des xi’s)
• La condition HB<<HA reste vérifiée si
1) L’objet 1 est beaucoup plus massif que les autres
2) Les distances mutuelles ne sont pas trop petites
(pas de rencontres proches)
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
79
Les problèmes de rencontres proches
• Lorsque deux corps passent proches l’un de l’autre, on a
des perturbations importantes et rapides.
• L’idéal est de diminuer le pas de temps dans ce type de
situation.
• Mais une intégration symplectique demande que le pas de
temps soit constant (Hinteg=f(τ))  difficulté ! Deux
variantes :
– On laisse tomber la symplecticité le temps de la rencontre proche
(code RMVS Levison & Duncan)
– On modifie le découpage HA+HB le temps de la rencontre proche.
On est toujours symplectique mais le temps de calcul est plus long
(code MERCURY Chambers)
– On découpe les potentiels en cercles concentriques avec un pas de
temps adapté à chaque tranche (Code SyMBA Duncan, Lee &
Levison).
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
80
Codes autogravitants
• On veut calculer un système N corps avec N grand
• La difficulté réside surtout dans le calcul des forces. Les différents types
de codes se différencient par la méthode choisie.
• Codes directs (PP): On calcule directement les N(N-1)/2~N2/2 termes de
forces. C’est lent et le temps de calcul est N2. On arrive à simuler
quelques 104 particules.
• Codes en arbre : tempe de calcul  N×ln N ~108 particules
– Les forces proches sont calculées directement
– Pour les termes lointains, les objets sont regroupés en groupes cubiques
appelée nœuds.
– La contribution des nœuds lointains est développée en harmoniques
sphériques et tronquée.
2
l 


 
 
ai    U ( rj  ri )    U ml ( rj  ri )
j
k l 0 m   l




 
termes directs
septembre 2011
termes lointains
Master 2 AMD - Méthodes
numériques pour l'astrophysique
81
Codes autogravitants (2)
•
Codes particule-maille (PM): On découpe l’espace en J3 cellules cubiques de
côté a, et on suppose que la masse dans chaque cellule est concentrée en son
centre. Le potentiel dans la cellule (i,j,k) se calcule comme
U i , j ,k  
•


 / a  terme d' adoucissem ent 
C’est une convolution discrète. En passant par l’espace de Fourier, il n’y a qu’à
multiplier les transformées de Fourier  temps de calcul en J×ln J
(~109 cellules)
Codes particule-particule/particule-malle (P3M) : Les codes PM manquent de
résolution pour traiter les forces proches. On divise alors le calcul en deux :
–
–
–
–
•
GM (i, j, k )
2
2
2
2
( i, j,k  ) ( i , j ,k ) a i   j   k    / a 

Les forces lointaines sont calculées par la méthode PM
Les termes proches sont calculées individuellement (PP)
On atteint ~107 particules
Variante : TMP (Tree/Particle Mesh), où la partie PP est calculée avec un code en
arbre. On atteint ~1010 particules !
Codes à grille adaptatives (AP3M) : On crée des sous-grilles plus fines dans la
partie P3M là où la densité de particules est plus grande (Codes MLAPM, ART,
RAMSES…)
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
82
Conditions aux limites en plusieurs points
• Cas général : on a un système différentiel de la
forme dY/dx = f(x,Y), où Y=(y1,…yn) est un
vecteur de dimension n, avec n1 conditions de la
forme gi(x1,Y)=0 (à satisfaire en x=x1) et n2
conditions hi(x2,Y)=0 (en x=x2≠x1) (n=n1+n2)
• Cas particulier classique : une équation scalaire du
second ordre y″=f(x,y,y′) avec y(x1)=y1 et
y(x2)=y2.
• Deux méthodes d’approche : méthodes de tir et
méthodes de relaxation.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
83
Méthodes de tir
• Idée : On va intégrer en partant de x1 à partir d’un point initial qui
vérifie les n1 premières conditions. On ajuste ensuite le point de départ
de manière à ce que les n2 autres en x=x2 soient satisfaites.
• Le point de départ est Y(x1), en vecteur à n composantes. Les n1
conditions en x=x1 laissent un « espace » de dimension n2 pour choisir
Y(x1). On dira Y(x1)=Y(x1,v1,…,vn2). On écrit V=(v1,…,vn2).
• Pour un choix de V, on intègre par une méthode classique jusqu’à
x=x2. On obtient Y(x2,V). On regarde ensuite fk=hk(x2,Y(x2,V)) pour
k=1,…,n2. (F=(f1,…,fn2))
• On va chercher à fixer V de manière à obtenir F=0  Résoudre un
système d’équations linéaires. Si on est capable de calculer les dérivées
partielles (∂Fi/∂Vj) (même numériquement), on peut utiliser NewtonRaphson. Sinon, on peut faire de la dichotomie.
• Inconvénient : Chaque essai nécessite une intégration, et même
plusieurs si on veut les dérivées partielles numériquement.
Numériquement, ça peut être long !
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
84
Méthodes de relaxation
• On remplace l’équation différentielle par un schéma de différences
finies :
– On divise l’intervalle en p segments de longueur h=(x2-x1)/p; on pose
ti=x1+i×h, i=0,…,p (t0=x1, tp=x2); on pose Yi=Y(ti).
– On approxime
dY  ti  ti 1  Y(ti )  Y(ti 1 ) Yi  Yi 1



dx  2 
ti  ti 1
h
– Et on remplace l’équation par
 t  t Y  Yi 1 
Yi  Yi 1  hf  i i 1 , i
 i  1,..., p
2
2


– Et les conditions aux limites s’écrivent
g j ( x1 , Y0 )  0 j  1,..., n1 et h j ( x2 , Yp )  0 j  1,..., n2
• Ceci donne un système non linéaire d’équations (dimension
p×n+n1+n2 = (p+1)×n) que l’on doit résoudre pour trouver Y0,…,Yp.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
85
Méthodes de relaxation (2)
• Cas particulier : une équation scalaire du second ordre
y″=f(x,y,y′), avec y(x1)=y0 et y(x2)=yp donnés.
– On ne transforme pas en système différentiel.
On se place sur les points ti
d2 y
y (ti 1 )  y (ti 1 )  2 y (ti ) yi 1  yi 1  2 yi
(
t
)


i
2
2
1
dx
h2
4 (ti 1  ti 1 )
– Et on remplace l’équation par
y y 

yi 1  yi 1  2 yi  h 2 f  ti , yi , i 1 i 1  i  1,..., p  1 avec y0 et y p donnés.
2h 

• Ca permet de diminuer la dimension du système
d’équations.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
86
Méthodes de relaxation : résolution
• On doit résoudre l’équation aux différences finies :
 t  t Y  Yi 1 
Yi  Yi 1  hf  i i 1 , i
 i  1,..., p
2
2


g j ( x1 , Y0 )  0 j  1,..., n1 et h j ( x2 , Yp )  0 j  1,..., n2
• On peut réécrire cela F(Y)=0 avec Y=(Y0,…,Yp) qui se
résout par Newton-Raphson (dimension n(p+1)). A chaque
pas, on cherche l’incrément ΔY à appliquer à Y :(JF =
Jacobienne de F)
• Yn(i-1)+j =Yi,j . On a
J F   Y  Y
Fi , j
 t  t Y  Yi 1 
Fi , j  Yi , j  Yi 1, j  hf  i i 1 , i
 0 si i  i ou i  1
 
2
Yi, j
 2

 La matrice JF est bloc-diagonale !
 Il faut tenir compte de cette forme pour résoudre le
système (élimination Gaussienne bloc par bloc) !
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
87
Méthodes de relaxation : plus loin
• Dans certains cas, il peut y avoir des conditions au limites
au milieu de l’intervalle considéré  Ca rajoute un bloc
spécial dans la milieu de la matrice.
• Dans d’autres situations, on cherche à ce que les p+1
points ti ne soient pas régulièrement répartis, mais puissent
être alloués de manière dynamique au cours du processus
de relaxation
– La solution = considérer le vecteur T=(t0,…,tp) comme un
ensemble de variables supplémentaires et donner des équations à
satisfaire.
– Les codes d’évolution stellaire fonctionnent de cette façon…
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
88
Méthodes numériques pour
l’astrophysique
• Techniques de base
• Estimateurs et statistique
• Modélisation de données
• Résolution numérique d’équations différentielles
• Equations aux dérivées partielles
• Types d’équations
• Equations elliptiques : Méthodes de relaxation
• Equations elliptiques : Surrelaxation
• Equations hyperboliques : schémas de Lax, leapfrog…
• Equations paraboliques : schémas FTCS, de Crank-Nicholson…
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
89
Equations aux dérivées partielles (PDEs)
• Une gamme de problèmes très vastes en astrophysique :
HD, MHD = codes eulériens.
• Méthode de base : différences finies, mais il y a d’autres
approches.
• Trois familles d’équations : hyperboliques, paraboliques,
elliptiques
• Deux types de problèmes :
– Hyperboliques, paraboliques = problèmes avec des conditions
initiales  propagation, diffusion = évolution temporelle.
 problème = stabilité du schéma de discrétisation
– Elliptiques = problèmes statiques avec des conditions aux limites
complexes
 problème = efficacité de la convergence vers la solution
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
90
Types d’équations aux dérivées partielles
 2u ( x, y )
 2u ( x, y )
 2u ( x, y )
A
B
C
  0
2
2
x
xy
y
• Equation hyperbolique B2-4AC>0 . Exemple l’équation des
2
ondes (c = vitesse de propagation)
 2u
2  u
t
2
c
x 2
• Equation parabolique B2-4AC=0 . Exemple l’équation de
diffusion (D = coefficient de diffusion) u   u 
t

D 
x  x 
• Equation elliptique B2-4AC<0 . Exemple l’équation de
Poisson (ρ = terme source)
 2u  2u
 2   ( x, y )
2
x
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
y
91
Equations elliptiques :
problèmes de conditions aux limites
• Dans ce type de problème, la forme des conditions aux limites compte au
moins autant que l’équation elle-même…
• Discrétisation : On introduit un réseau de points (xj,yl), pour j=0...J,
l=0…L, réalisant un maillage de l’espace considéré.
x j  x0  jh ( j  0,, J ),
• On approxime
u j 1,l  2u j ,l  u j 1,l
 2u
(
x
,
y
)

j
l
2 x2
h2
et
yl  y0  lh (l  0,, L)
u j ,l 1  2u j ,l  u j ,l 1
 2u
(
x
,
y
)

j
l
2 y 2
h2
• L’équation de Poisson se réduit à
u j 1,l  u j 1,l  u j ,l 1  u j ,l 1  4u j ,l  h2  j ,l
pour 1  j  J  1 et 1  l  L-1
plus des conditions aux limites pour j=0, j=L, l=0, l=L .
• Mis bout à bout, la résolution se ramène à un système linéaire
A  u  b avec u  (u0,0 , u0,1 ,, u0,L , u1,0 ,, u1,L , u2,1,,, uJ ,L )
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
92
Equations elliptiques :
Résolution du système
• La résolution de l’équation aux différences finies se ramène donc à la
résolution d’un système linéaire. Simple ?
• OUI MAIS : La taille du vecteur u est (J+1)(L+1)  La taille de la
matrice A est (J+1)(L+1)×(J+1)(L+1) . Si J=L=100, u est de taille
10000 et la matrice A contient 108 éléments !! Impossible par des
méthodes générales…
• En fait la matrice A est très creuse : tridiagonale avec bordures
• Toute méthode de résolution doit tenir compte explicitement de la
forme de la matrice
• Les méthodes directes (pivot, LU, gradient conjugué…) ne sont
efficaces que pour des tailles modérées (<300×300) de grilles, et si on
a l’espace mémoire suffisant pour stocker la matrice…
• Autrement, les bonnes méthodes sont les méthodes de relaxation
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
93
Equations elliptiques :
Méthodes de relaxation
• En général, une équation elliptique du second ordre se réduit toujours
à une équation de la forme
a j ,l u j 1,l  b j ,l u j 1,l  c j ,l u j ,l 1  d j ,l u j ,l 1  e j ,l u j ,l  f j ,l
équivalente au système A.u=b. Idée : On décompose A en E-F où E
est facilement inversible (diagonale ou tridiagonale)
(E  F)  u  b  E  u  F  u  b
• On part d’un choix initial u(0) et itère l’équation par une procédure de
point fixe:
E  u( k 1)  F  u( k )  b
• Dans la pratique on a A = L+D+U (L= diagonale inférieure, D =
diagonale, U = diagonale supérieure)
– Méthode de Jacobi : E = D, F=-(L+U);
– Méthode de Gauss-Seidel : E = L+D, F=U
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
94
Equations elliptiques :
Méthodes de relaxation
• Les deux méthodes sont convergentes mais lentement. Gauss-Seidel
est un peu plus efficace.
• L’erreur décroît comme ρs-k , où ρs = rayon spectral de la matrice
E-1.F. Plus la grille est grande, plus ρs est proche de 1 
convergence lente ! Typiquement,  s  1  aJ 2   sk ~ 1  kaJ 2
• Solution : méthode de surrelaxation (SOR). On part de Gauss-Seidel
(L  D)  u ( k 1)  U  u ( k )  b 

u ( k 1)  (L  D) 1.  U  u ( k )  b



 u ( k )  (L  D) 1.  (L  D  U)  u ( k )  b  u ( k )  (L  D) 1 w ( k )
avec w(k) = vecteur résidu, et on remplace par
avec 1<ω<2
septembre 2011
u( k 1)  u( k )  (L  D)1  w( k )
Master 2 AMD - Méthodes
numériques pour l'astrophysique
95
Equations elliptiques :
Méthode de surrelaxation (SOR)
•
•
Avec un bon choix de ω, la méthode converge plus vite. On montre que le
choix optimal est
2
op 
2
1  1   Jacobi
Dans ce cas le rayon spectral de la méthode SOR est
SOR

 Jacobi

1  1-  2
Jacobi

 Jacobi  1 
•
•
•
a




2
 op  2 
2 2a
J
 SOR  1 
2 2a
J
J2
La convergence est plus rapide !
Problème : le gain n’est réel que si ω ≈ ωop. Or il est très difficile de
connaître ρJacobi (dépend du problème et de ses conditions aux limites)
Quand on ne sait pas, on prend une valeur standard (Neumann-Dirichlet)
1
2
 Jacobi   cos
septembre 2011


 cos 
J
L
Master 2 AMD - Méthodes
numériques pour l'astrophysique
96
Equations elliptiques :
Méthode de surrelaxation (SOR) : Formules pratiques
•
Le but est de dégager des formules pratiques pour la SOR. On reprend l’équation
générale
a j ,l u j 1,l  b j ,l u j 1,l  c j ,l u j ,l 1  d j ,l u j ,l 1  e j ,l u j ,l  f j ,l
•
Le résidu wj,l s’écrit
w j ,l  a j ,l u j 1,l  b j ,l u j 1,l  c j ,l u j ,l 1  d j ,l u j ,l 1  e j ,l u j ,l  f j ,l
•
La formule d’itération
u
( k 1)
u
(k )
  ( L  D)  w
1
(k )
 u
( k 1)
j ,l
u
(k )
j ,l

w j ,l
e j ,l
à condition de prendre les variables dans la bon ordre pour faire l’inversion de L+D
DO J = 2,….
DO L = 2,…
W = A(J,L)*U(J+1,L)+B(J,L)*U(J-1,L)+C(J,L)*U(J,L+1)+D(J,L)*U(J,L-1)+E(J,L)*U(J,L)-F(J,L)
U(J,L) = U(J,L)-OMEGA*W/E(J,L)
END DO
END DO
•
Amélioration : accélération de Tchébychev. On prend ω=1 au départ, et on le change à
chaque itération pour le faire tendre vers ω=ωop
 ( 0)  1,  (1) 
septembre 2011
1
1   Jacobi / 2
,  ( k 1) 
1
1   ( k )  Jacobi / 4
Master 2 AMD - Méthodes
numériques pour l'astrophysique
pour k  1
97
Equations de propagation (hyperboliques)
•
•
On cherche une fonction u(x,t) satisfaisant une équation hyperbolique; on
connaît u(x,0) (condition initiale). On cherche à connaître l’évolution dans le
temps.
On va discrétiser spatialement et temporellement. Par exemple, si u(x,t), on
écrit
n
x j  x0  jx ( j  0,, J ), tn  t0  nt (n  0,1,), u j  u( x j , tn )
•
•
Un schéma de discrétisation est un formule permettant de calculer ujn+1
(j=1,…,J) en fonction des ujn (j=1,…,J) et éventuellement ujn-1 (j=1,…,J)
Important : La stabilité d’un schéma (von Neumann). Un schéma doit être
stable pour être praticable. Un mode propre du schéma c’est
unj   neikjx
•
•
(k  nombre d' onde;    (k )  facteur d' amplificat ion)
ujn+1/ujn=ξ . Stabilité  |ξ(k)| ≤ 1 pour tout k (sinon on pourrait trouver un
mode croissant exponentiellement)
On injecte la forme d’un mode propre dans l’équation du schéma pour trouver
ξ(k)…
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
98
Equations de conservation de flux
• La plupart des équations de propagation peuvent se réécrire comme
des équation de conservation de flux
u
F(u)

t
x
   F en multidimen
• Exemple : l’équation des ondes
u
u
 c2 2
2
t
x
2
2
s
 r

c

x
  t
s
r
 c
x
 t
avec
u

r

c

x

u
s 
t


sions  (F  flux)
u
F

t
x
avec
u   r  et F(u)   0  c   u
v
 c 0 
• On va considérer une équation conservative scalaire tout simple
u
u
 c
t
x
et la discrétisation
septembre 2011
u
t

j ,n
u nj 1  u nj
t
u
;
x

j ,n
Master 2 AMD - Méthodes
numériques pour l'astrophysique
u nj1  u nj1
2x
99
Schéma explicite FTCS
u
t

u nj 1  u nj
t
j ,n
u
;
x

u nj1  u nj1
j ,n
2x
• On injecte dans l’équation
u nj 1  u nj
t
 u nj1  u nj1 
 u nj1  u nj1 
n 1
n
  u j  u j  ct 

 c



 2x 
 2x 
• On obtient un schéma explicite dit FCTS (Forward Time Centered
Space)
• Est-il stable ? On cherche les modes propres….
 n1eikjx   neikjx
t
  neik ( j 1) x   neik ( j 1) x 
ct
   (k )  1  i
 c
sin kx
2

x

x


• |ξ(k)| > 1 pour tout k  FCTS est toujours instable !!!
Il faut trouver mieux…
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
100
Schéma de Lax
• On remplace ujn par (uj+1n+uj-1n)/2 dans la dérivée temporelle
u
t

2u nj 1  u nj1  u nj1
j ,n
2t
u
;
x
• Stabilité ? On trouve

u nj1  u nj1
j ,n
2x
u
 (k )  cos kx  i
• On veut |ξ(k)| ≤ 1 pour tout k 
n 1
j

u nj1  u nj1
2

ct n

u j 1  u nj1 
2x
ct
sin kx
x
c t
1
x
• On appelle cela la condition de Courant.
Interprétation : On calcule ujn+1 (en x=xj) à l’aide de uj+1n et uj-1n (en
x=xj-1 en x=xj+1 ). L’information se propage à la vitesse maximale ±c.
Les points uj+1n et uj-1n doivent être en dehors de zone liée à ujn+1.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
101
Schéma de Lax
• Dissipation
 c 2 t 2  2
 sin kx
 (k )  1  1 
2 
x 

• Si |c|Δt < Δx, |ξ(k)| < 1  L’amplitude des modes décroît.
Il y a de la viscosité numérique ! Interprétation :
2
 u nj1  u nj1  1  u nj1  u nj1  2u nj 
u nj 1  u nj
u
u x   2u
 
 
  c
 c 
t


2x
 2


t


t
x
2
2

t

x


terme de diffusion
• Tout se passe comme si on avait ajouté un terme de
diffusion (parabolique) à l’équation initiale…
• Dans la pratique l’effet est faible car |kΔx|<<1 (longueur
d’onde >> Δx).
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
102
Ordre 2 en temps : Staggered Leapfrog
• On va centrer le calcul de ∂u/∂t.
u
t

j ,n
u nj 1  u nj 1
2t
u
;
x

j ,n
u nj1  u nj1
2x
 u nj 1  u nj 1  
ct n

u j 1  u nj1 
x
• Le schéma est du second ordre spatialement et temporellement. Mais on
a besoin de l’information à tn et tn-1 pour calculer à tn+1
• Stabilité : On cherche les modes propres. On tombe sur
ct
ct
 ct

  2i
sin kx  1  0    i
sin kx  1  
sin kx 
x
x
 x

2
2
• Si |c|Δt/Δx ≤ 1 (Condition de Courant…) la racine est réelle on a |ξ| = 1
(stabilité); sinon ξ est imaginaire pur avec |ξ| > 1
• |ξ| = 1  pas de viscosité numérique ! C’est ce qui fait l’intérêt de
cette méthode.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
103
Staggered Leapfrog : Equations du second ordre
•
Dans le cas d’une équation de conservation de flux
le leapfrog s’écrit
unj 1  unj 1  
•
t n

Fj 1  Fjn1  avec
x
Fjn  F(unj )
On écrit ce schéma pour l’équation des ondes (r = c∂u/∂x, s = ∂u/∂x).
 n 1 n 1 ct n
n
rj  rj  x s j 1  s j 1 
avec

ct n
n 1
n 1
n
 sj  sj 
rj1  rj1 
x

•
u
F(u)

t
x
n
 n
u nj1  u nj1
u
c
rj  c

x
2x

j

n
n 1
n 1
 s n  u  u j  u j
 j t j
2t

n2
On remplace. La deuxième équation donne
uj
 2u nj  u nj 2
2( t )2
c
2
u nj2  2u nj  u nj2
2( x )2
ce qui est complètement équivalent à la discrétisation directe de l’équation des
ondes avec 2Δt et 2Δx. On écrira donc le leapfrog pour ces équations sous la
forme classique
u nj 1  2u nj  u nj 1
( t )
septembre 2011
2
c
2
u nj1  2u nj  u nj1
( x )2
Master 2 AMD - Méthodes
numériques pour l'astrophysique
104
Schéma en deux temps de Lax-Wendroff
• Quand les équations sont plus complexes, le leapfrog peut
être instable car il couple les points deux par deux. Dans
ce cas on utilise le schéma de Lax-Wendroff :
– On introduit des points intermédiaires (xj+1/2, tn+1/2 ). On calcule la
valeur correspondante de u par le schéma de Lax. On en tire le
unj1  unj
flux correspondant
t n
n 1 / 2
n
n 1 / 2
u j 1 / 2 
2

2x
F
j 1
 Fj   Fj 1 / 2
– On en tire les flux correspondants et on utilise le leapfrog pour
tirer ujn+1
t
unj 1  unj   Fjn11//22  Fjn11//22 
x
• Dans le cas de l’équation simple (F=cu), l’ensemble se
réduit à n1 n a n n
ct
n
n
n
uj  uj 
septembre 2011
2
u
j 1

 u j 1  a u j 1  2u j  u j 1 
Master 2 AMD - Méthodes
numériques pour l'astrophysique
avec a 
x
105
Schéma en deux temps de Lax-Wendroff
• Stabilité ? On trouve
  1  ia sin kx  a 2 1  cos kx     1  a 2 1  a 2 1  cos kx 2
• |ξ| ≤ 1  |α| ≤ 1 (Condition de Courant !).
• En général quand α ≠ 1, |ξ| < 1. Il y a de la dissipation…
• L’effet est plus faible que dans le schéma de Lax. Quand
|kΔx|<<1, on a
(kx )4
  1  a 1  a 
8
2
alors que dans le schéma de Lax on a
2
  1  1  a 2 (kx)2
• Quelle stratégie adopter ? Pour les problèmes qui se mettent sous la
forme d’une conservation d’un flux (type ondes), adopter d’abord le
leapfrog. S’il y a des problèmes d’instabilité, passer au schéma de
Lax-Wendroff.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
106
Raffinement : dérivation « face au vent »
• Certaines équations (advection…) sont sensibles aux problèmes de
transport (passage de chocs, changement d’états…). Si c>0,
l’information sur ujn+1 ne peut pas venir de unj+1 (et vice versa si
c<0). Dans ce cas on préfère calculer ∂u/∂x comme ceci :
u
x
j ,n
 u nj  u nj1

  n x n
 u j 1  u j
 x
si c  0 
u nj1  u nj1 
 au lieu de


2x 

si c < 0
• On appelle cela dérivation « face au vent » (upwind)
• Du coup la discrétisation spatiale n’est plus que du premier ordre.
Mais ça peut être plus stable dans certains cas… A utiliser quand
c’est nécessaire !
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
107
Equations de diffusion (paraboliques)
u   u 
 D 
t x  x 
•
Equations du type
•
D = coefficient de diffusion. C’est une équation de conservation de flux avec
F = -D(∂u/∂x). Mais on préfère en général utiliser des méthodes spécifiques.
2

u

u
Si D est constant, on a l’équation de la chaleur
D 2
avec une discrétisation immédiate
t
x
•
u nj 1  u nj
t
•
•
 u nj1  2u nj  u nj1 
 D

2
(

x
)


C’est un schéma explicite de type FTCS, mais plus stable…
4 Dt
( x )2
2  kx 
  1
sin 
 ;   1  t 
2
( x )
2D
 2 
Interprétation : Le pas de temps doit être plus petit que le temps de diffusion
à travers une (demi) cellule.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
108
Equations de diffusion : schémas implicites
•
•
•
•
La condition Δt ≤ (Δx)2/(2D) est très contraignante  temps de calcul
prohibitif
Il faut chercher des schémas plus stables. Solution : schémas implicites.
Exemple 1 : on évalue ∂2u/∂x2 en t = tn+1 au lieu de t = tn.
 u nj11  2u nj 1  u nj11 
u nj 1  u nj
 D

2
t
(

x
)


C’est un schéma implicite = il faut résoudre un système linéaire (tridiagonal)
pour trouver les un+1j (j=0,…,J)
 a u nj11  (1  2a )u nj 1  a u nj11  u nj
•
•
j  1,..., J  1 ; a 
Dt
( x )2
+conditions aux limites en j=0 et J-1
Stabilité :
1

   1 toujours !
k

x


1  4a sin 2 

2


Mais c’est un schéma du premier ordre  Mieux : Crank-Nicholson !
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
109
Equations de diffusion :
schéma de Crank-Nicholson
•
On prend la moyenne du schéma explicite et du schéma implicite
u nj 1  u nj
t
•
•
•
n 1
n 1
n 1
n
n
n
D  u j 1  2u j  u j 1   u j 1  2u j  u j 1 
 

2
( x )2

C’est encore un schéma implicite. Mais il est du second ordre en temps
(centré en tn+1/2)
 kx 
Stabilité :
1  2a sin 2 

2 


   1 toujours !
k

x


1  2a sin 2 

 2 
C’est le schéma recommandé pour tous les problèmes de diffusion…

a
2
u nj11  (1  a )u nj 1 
septembre 2011
a
2
u nj11 
a
2
u nj1  (1  a )u nj 
a
2
Master 2 AMD - Méthodes
numériques pour l'astrophysique
u nj1
j  1,..., J  1
110
Equations de diffusion plus complexes
•
Si D n’est pas une constante… On discrétise D
Si D(x), on appelle Dj+1/2=D(xj+1/2). Le schéma FTCS devient
u nj 1  u nj
t
D j 1 / 2 u nj1  u nj   D j 1 / 2 u nj  u nj1 

( x )2
avec la condition de stabilité
•
Crank-Nicholson devient
u nj 1  u nj
t
•

 ( x ) 2 
t  min 

j
 2 D j 1 / 2 
D j 1 / 2 u nj11  u nj 1  u nj1  u nj   D j 1 / 2 u nj 1  u nj11  u nj  u nj1 
2( x )2
Les vraies difficultés commencent quand l’équation n’est pas linéaire :
D(x,u) . Les schémas implicites deviennent non-linéaires. Si on est capable
de calculer z = ∫D(u) du, le membre de droite de l’équation = ∂2z/∂x2 qu’on
linéarise. On calcule zjn+1 en faisant un développement limité au 1er ordre.
septembre 2011
Master 2 AMD - Méthodes
numériques pour l'astrophysique
111