M thodes num riques pour l astrophysique m2 astrophysique et milieux dilu s
This presentation is the property of its rightful owner.
Sponsored Links
1 / 111

Méthodes numériques pour l’astrophysique M2 « Astrophysique et Milieux Dilués » PowerPoint PPT Presentation


  • 68 Views
  • Uploaded on
  • Presentation posted in: General

Méthodes numériques pour l’astrophysique M2 « Astrophysique et Milieux Dilués ». Hervé Beust. Laboratoire d’Astrophysique de Grenoble. Méthodes numériques pour l’astrophysique. Techniques de base Estimateurs et statistique Modélisation de données

Download Presentation

Méthodes numériques pour l’astrophysique M2 « Astrophysique et Milieux Dilués »

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Presentation Transcript


M thodes num riques pour l astrophysique m2 astrophysique et milieux dilu s

Méthodes numériques pour l’astrophysiqueM2 « Astrophysique et Milieux Dilués »

Hervé Beust

Laboratoire d’Astrophysique de Grenoble

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes num riques pour l astrophysique

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

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes num riques pour l astrophysique1

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

Master 2 AMD - Méthodes numériques pour l'astrophysique


R solution num rique d quations

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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Diff rents cas de figure

Différents cas de figure...

Cas standard

 Cas plus délicats 

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thode de dichotomie

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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Autres m thodes d ordre 1

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…

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thode de newton raphson

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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thode de newton raphson probl mes

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

Master 2 AMD - Méthodes numériques pour l'astrophysique


Am lioration newton raphson quartique

Amélioration : Newton-Raphson quartique

SUBROUTINE KEPLER0(M,E,U) 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 U = M OK = .FALSE. I=0DO WHILE (.NOT.OK)F = U-E*SIN(U)-MF1 = 1.d0-E*COS(U) F2 = U-M-FF3 = 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+1OK = ((ABS(DELA).LT.EPSI).OR.(I.GE.500))END DOEND

  • On peut améliorer la convergence si on sait calculer des dérivées d’ordre supérieur

  • Exemple : équation de KéplerM = u - e sin u

Master 2 AMD - Méthodes numériques pour l'astrophysique


En multidimensions

En multidimensions …

  • En dimensions plus que 1, c’est beaucoup plus difficile.

  • La seule méthode générale, c’est Newton-Raphson, sans garantie de succès…

  • Jxn = Matrice jacobienne en x

  • 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


Int gration num rique

Intégration numérique

  • Problème : calculer numériquement une intégrale de la forme

  • 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 pointsxi , i=1..n

  • Tout repose dans le choix des coefficientsaiet des pointsxi.

  • Les formules convergent vers l’intégrale lorsque n ∞

Master 2 AMD - Méthodes numériques pour l'astrophysique


Formules points r guli rement espac s

Formules à points régulièrement espacés

  • Formule des trapèzes :

  • Formule de Simpson :

Master 2 AMD - Méthodes numériques pour l'astrophysique


Quadratures de gauss polyn mes orthogonaux

Quadratures de Gauss – polynômes orthogonaux

  • Idée : ne plus imposer que les xi soient régulièrement espacés.

  • Choisir les xiet les ai de manière à ce que la formule d’approximation soit exacte pour des polynômes de degré leplus é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 xiet 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

  • La suite des polynômes Pn lorsque n varie est orthogonale au sens du produit scalaire suivant

Master 2 AMD - Méthodes numériques pour l'astrophysique


Quadratures de gauss polyn mes orthogonaux1

Quadratures de Gauss – polynômes orthogonaux

  • Les Pnobéissent à la récurrence suivante :

  • Les xise calculent comme les racines de Pn et les aivérifient :

Master 2 AMD - Méthodes numériques pour l'astrophysique


Quadratures de gauss polyn mes orthogonaux2

Quadratures de Gauss – polynômes orthogonaux

Cas particuliers :

Dans le cas de Gauss-Legendre, on a en outre

Et pour tout intervalle (a,b) on a

Master 2 AMD - Méthodes numériques pour l'astrophysique


Exemple r solu gauss legendre

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

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


Int gration multidimensionnelle

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é

  • Dans les autres cas : Intégration Monte-Carlo

Master 2 AMD - Méthodes numériques pour l'astrophysique


Int gration monte carlo

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

  • On tire au hasard un nombre n de points dans l’ensemble et on approxime la moyenne par

  • L’erreur moyenne sur l’intégrale esta convergence en n-1/2, lente

  • 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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Minimisation de fonctions

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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes en dimension 1

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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thode interpolation par paraboles inverses

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

  • Le nouvel encadrement est (a,d,b) ou (b,d,c).

  • La méthode de Brent combine les deux méthodes précédentes.

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes multidimensionnelles

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).aOn 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.aLes directions n et n’ doivent être conjuguées.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Directions congugu es

Directions conguguées

  • On se place en un point P

  • Le changement de gradient dans un déplacement infinitésimal vaut

  • 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.

  • Les directions u et v sont alors dites conjuguées.

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thode de powell

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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes avec gradient

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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thode de fletcher reeves polak ribiere

Méthode de Fletcher-Reeves-Polak-Ribiere

  • On part d’un point x0, on fixe une direction u0, et on calcule

  • A chaque étape on dispose de

  • On minimise le long de ui. On appelle xi+1 le nouveau point trouvé. On calcule

  • On calcule ui+1 comme

  • 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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes num riques pour l astrophysique2

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

Master 2 AMD - Méthodes numériques pour l'astrophysique


Estimateurs et statistiques

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

  • On dit que l’estimateur est convergent si pour tout ε>0, on a

  • La loi faible des grands nombre montre que Xn est un estimateur convergent de l’espérance.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Propri t s des estimateurs

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 θ.

  • On pose Y=ln X. On a alors E(Y)=ln θ – 1 (calcul). Doncest un estimateur convergent de ln θ – 1 . Et ensuiteest un autre estimateur convergent de θ

  • On peut en créer beaucoup d’autres. Comment sélectionner le meilleur ?

Master 2 AMD - Méthodes numériques pour l'astrophysique


Propri t s des estimateurs1

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 Tnest 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∞.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Exemples

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/naTn 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)

Master 2 AMD - Méthodes numériques pour l'astrophysique


Convergence des estimateurs

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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Estimateurs standards

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)

Master 2 AMD - Méthodes numériques pour l'astrophysique


Recherche d estimateurs

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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Estimation par ajustement

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 :

  • Distance de Kolmogorov-Smirnov = comparaison entre les fonctions de répartition théoriques (F) et empiriques (Fe)

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes num riques pour l astrophysique3

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

Master 2 AMD - Méthodes numériques pour l'astrophysique


Mod lisation de donn es par r gression

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

  • Toute la difficulté consiste à

    • Trouver les paramètres réalisant le minimum

    • Estimer l’incertitude sur les paramètres.

Master 2 AMD - Méthodes numériques pour l'astrophysique


R gression lin aire

Régression linéaire

  • But : modéliser des données par une relation linéairey = ax+b

  • Méthode : Minimiser le χ2pour trouver a et b

Master 2 AMD - Méthodes numériques pour l'astrophysique


R gression lin aire 2

Régression linéaire (2)

  • Erreur sur a et b :

    • Var(yi)=σi2 Var(aiyi)=ai2 σi2(Somme de variables aléatoires indépendantes)

Master 2 AMD - Méthodes numériques pour l'astrophysique


Mod les lin aires de moindres carr es

Modèles linéaires de moindres carrées

  • On cherche un modèle qui dépend linéairement des paramètresy = f(α1,…,αk; x) = α1X1(x)+∙∙∙+αkXk(x)où les X1,…,Xk sont des fonctions données à l’avance

  • On va chercher à minimiseren disant ∂χ2/∂αi=0 pour i=1…k

  • On introduit la matrice A, n×k, et le vecteur b de longueur n tels que

Master 2 AMD - Méthodes numériques pour l'astrophysique


Mod les lin aires de moindres carr s

Modèles linéaires de moindres carrés

  • ∂χ2/∂αj=0 pour j=1…k revient à

    …autrement dit résoudre un système linéaire ! (équations normales)

  • On appelle fij = [c]-1ij. On a

Master 2 AMD - Méthodes numériques pour l'astrophysique


Mod les non lin aires

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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thode de levenberg marquardt

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

    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 on part d’un vecteur α0 , on a

  • 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 gradient

  • Idée de la méthode : Alterner entre les deux formules

  • Problème : Il faut la matrice H = matrice Hessienne

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thode de levenberg marquardt 2

Méthode de Levenberg-Marquardt (2)

  • Problème : Les termes cjldépendent des dérivées secondes de f qui ne sont pas forcément accessibles

  • Souvent on laisse tomber les termes en question…

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thode de levenberg marquardt 3

Méthode de Levenberg-Marquardt (3)

  • Deux questions :

    • Comment décider de l’alternance entre les deux formules (1) et (2) ?

    • Comment fixer la constante dans la formule (2) ?

  • Levenberg-Marquardt :

    • 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

    • Quand λ>>1, (3) revient à (2); quand λ<<1, (3) équivaut à (1)

  • Recette :

    • 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 da. On calcule χ2(α+da).

    • Si χ2(α+da) > χ2(α), on est loin du minimum  on multiplie λ par 10 et on revient à 2.

    • Si χ2(α+da) < χ2(α), on approche  On divise λ par 10, on remplace αparα+da, et on revient à 2.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Intervalles de confiance des param tres

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…

Master 2 AMD - Méthodes numériques pour l'astrophysique


Avec la matrice des covariances

Avec la matrice des covariances…

  • Si les erreurs sont gaussiennes, alors on a 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…

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes statistiques ou quand les erreurs ne sont pas gaussiennes

Méthodes statistiquesou 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).

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thode mcmc monte carlo cha nes de markov

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 , on connaît , probabilité d’observer y sachant . On cherche

  • Par le théorème de Bayes

  • On suppose connaître (prior).

  • Une chaîne de Markov = une suite de modèles i, où chaque i+1se déduit de ivia une probabilité de transition . A la fin, la distribution des i échantillonne

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thode mcmc 2

Méthode MCMC (2)

  • Définition : La chaîne de Markov est réversible si

  • Théorème : Si la chaîne est réversible, apériodique et irréductible, ell converge vers

  • Problème : Il n’est pas facile de trouver une probabilité de transition réversible. On peut en construire une à partir d’une autre qui ne l’est pas, par l’algorithme de Metropolis-Hastings

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thode mcmc 3

Méthode MCMC (3)

  • Méthode à chaque pas :

    • Connaissant i, on tire au hasard un i+1en suivant

    • On calcule 2(i+1) et 2(i).

    • On détermine le rapport

    • On tire un nombre u au hasard entre 0 et 1.

    • On calcule . Si u<, on accepte le pas, sinon on revient à i.

  • Souvent, on prend qtr symétrique dans une direction v. On fait tourner les directions.

  • On peut prendre qtrnon 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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thode mcmc 4

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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Probl mes de minima locaux

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 non-liné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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Recuit simul

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 T0est é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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Algorithmes g n tiques

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 :

    • 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 (Levenberg-Marquardt…) sur ma population et j’obtiens ma génération suivante. Puis retour en 1…

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes num riques pour l astrophysique4

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

Master 2 AMD - Méthodes numériques pour l'astrophysique


R solution num rique d quations diff rentielles

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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations diff rentielles ordinaires

Equations différentielles ordinaires

  • C’est résoudre le problème

  • Remarques :

    • 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 :

    • avec une condition aux limites de la formey(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)

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes de type runge kutta

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)

  • 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

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes de type runge kutta 2

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)

  • Runge-Kutta d’ordre 4: yn+1=y(xn+1)+O(h5)

  • Remarque : Si le système est vectoriel les ki’s sont des vecteurs…

Master 2 AMD - Méthodes numériques pour l'astrophysique


Contr le du pas

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

  • Stratégie :

    • Si Δ>Δ0, on recommence en diminuant le pas. On prend

    • Si Δ<Δ0, on peut augmenter le pas en posant

Master 2 AMD - Méthodes numériques pour l'astrophysique


Contr le du pas 2

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

  • Si un pas est bon, quel valeur choisir ? y1 ou y2 ? Une combinaison…

  • 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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Contr le du pas 3 felhberg cash carp

Contrôle du pas (3) : Felhberg-Cash-carp

  • La formule d’ordre 5 est de la forme

  • La formule d’ordre 4 utilise les mêmes coefficients aiet bij, mais des coefficients c′i différents

Master 2 AMD - Méthodes numériques pour l'astrophysique


Contr le du pas 4 felhberg cash carp

Contrôle du pas (4) : Felhberg-Cash-carp

  • L’estimation de l’erreur se fait au moyen de

  • Δ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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thode de bulirsch et stoer

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 deh (ou n=H/h) et extrapoler cette fonction à h→0 (ou n→∞).

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thode du point m dian modifi e ou diff rences centr es

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 Bulirsh-Stoer. On veut aller de x0 à x0+H en n pas de h=H/n.

  • 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 :

Master 2 AMD - Méthodes numériques pour l'astrophysique


Extrapoler

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é k-1 (interpolation de Lagrange) qui passe par les k points de mesure que l’on calcule de la manière suivante, et on évalue en X=0.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Extrapoler 2

Extrapoler (2) ?

  • Extrapolation rationnelle : On utilise des fractions rationnelles au lieu de polynômes. Les résultats sont en général meilleurs

  • 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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations mal conditionn es stiff

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’Eulerqui 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)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

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations mal conditionn es stiff1

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

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes implicites

Méthodes implicites

  • Exemple 1toujours stable !

  • Exemple 2toujours 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

    • Il est difficile de définir des schémas implicites d’ordre supérieur Méthodes semi-implicites de Rosenbrock

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes semi implicites

Méthodes semi-implicites

  • On approxime…

  • 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 supérieur

  • 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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Codes n corps

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

Master 2 AMD - Méthodes numériques pour l'astrophysique


Int gration symplectique

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) :

  • Pour toute quantité q(xi,pi)

  • La solution donnant q(t) à partir de q(t-t) est

Master 2 AMD - Méthodes numériques pour l'astrophysique


Int gration symplectique 2

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 exactementH ou un autre Hamiltonien voisin de H. Du coup, l’erreur faite sur Hest bornée

Master 2 AMD - Méthodes numériques pour l'astrophysique


Int gration symplectique 3

Intégration symplectique (3)

  • Hypothèses de base de l’intégration symplectique : H = HA+HB(F = A+B) où on sait intégrer HAet HB(on sait calculer exp(τA) et exp(τB))

  • Souvent, on suppose en plus HB<<HA (HB/HA ~ ε <<1)

  • On a alors 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

Master 2 AMD - Méthodes numériques pour l'astrophysique


Int gration symplectique 4

Intégration symplectique (4)

  • On a aussia On obtient un intégrateur symplectique d’ordre 2 en intégrant 1) HB pendant t/2,2) HApendant t, 3) HBpendant t/2.Dans ce cas,

  • Il y a aussi des méthodes d’ordre 4, 6, 8…

  • En fait, compte tenu de HB/HA~ε on a mêmea On peut intégrer avec un grand pas de temps

Master 2 AMD - Méthodes numériques pour l'astrophysique


Int grateurs symplectiques pour codes n corps

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)

Master 2 AMD - Méthodes numériques pour l'astrophysique


Les probl mes de rencontres proches

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 RMVSLevison & Duncan)

    • On modifie le découpageHA+HB le temps de la rencontre proche. On est toujours symplectique mais le temps de calcul est plus long (code MERCURYChambers)

    • On découpe les potentiels en cercles concentriques avec un pas de temps adapté à chaque tranche (Code SyMBADuncan, Lee & Levison).

Master 2 AMD - Méthodes numériques pour l'astrophysique


Codes autogravitants

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 ~108particules

    • 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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Codes autogravitants 2

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 commeC’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 :

    • 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…)

Master 2 AMD - Méthodes numériques pour l'astrophysique


Conditions aux limites en plusieurs points

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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes de tir

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 n1conditions 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 Newton-Raphson. 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 !

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes de relaxation

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

    • Et on remplace l’équation par

    • Et les conditions aux limites s’écrivent

  • 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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes de relaxation 2

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

    • Et on remplace l’équation par

  • Ca permet de diminuer la dimension du système d’équations.

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes de relaxation r solution

Méthodes de relaxation : résolution

  • On doit résoudre l’équation aux différences finies:

  • 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  La matrice JF est bloc-diagonale ! Il faut tenir compte de cette forme pour résoudre le système (élimination Gaussienne bloc par bloc) !

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes de relaxation plus loin

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 tine 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…

Master 2 AMD - Méthodes numériques pour l'astrophysique


M thodes num riques pour l astrophysique5

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…

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations aux d riv es partielles pdes

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

Master 2 AMD - Méthodes numériques pour l'astrophysique


Types d quations aux d riv es partielles

Types d’équations aux dérivées partielles

  • Equation hyperbolique B2-4AC>0 . Exemple l’équation des ondes (c = vitesse de propagation)

  • Equation parabolique B2-4AC=0 . Exemple l’équation de diffusion (D = coefficient de diffusion)

  • Equation elliptique B2-4AC<0 . Exemple l’équation de Poisson (ρ = terme source)

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations elliptiques probl mes de conditions aux limites

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é.

  • On approxime

  • L’équation de Poisson se réduit à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

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations elliptiques r solution du syst me

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

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations elliptiques m thodes de relaxation

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

    équivalente au système A.u=b. Idée : On décompose A en E-F où E est facilement inversible (diagonale ou tridiagonale)

  • On part d’un choix initial u(0) et itère l’équation par une procédure de point fixe:

  • 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

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations elliptiques m thodes de relaxation1

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 matriceE-1.F. Plus la grille est grande, plus ρs est proche de 1  convergence lente !Typiquement,

  • Solution : méthode de surrelaxation (SOR). On part de Gauss-Seidelavec w(k) = vecteur résidu, et on remplace par avec 1<ω<2

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations elliptiques m thode de surrelaxation sor

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

  • Dans ce cas le rayon spectral de la méthode SOR est

  • 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)

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations elliptiques m thode de surrelaxation sor formules pratiques

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

  • Le résidu wj,l s’écrit

  • La formule d’itération à condition de prendre les variables dans la bon ordre pour faire l’inversion de L+DDO 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 DOEND 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

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations de propagation hyperboliques

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

  • 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

  • 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)…

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations de conservation de flux

Equations de conservation de flux

  • La plupart des équations de propagation peuvent se réécrire comme des équation de conservation de flux

  • Exemple : l’équation des ondes

  • On va considérer une équation conservative scalaire tout simpleet la discrétisation

Master 2 AMD - Méthodes numériques pour l'astrophysique


Sch ma explicite ftcs

Schéma explicite FTCS

  • On injecte dans l’équation

  • On obtient un schéma explicite dit FCTS (Forward Time Centered Space)

  • Est-il stable ? On cherche les modes propres….

  • |ξ(k)| > 1 pour tout k FCTS est toujours instable !!! Il faut trouver mieux…

Master 2 AMD - Méthodes numériques pour l'astrophysique


Sch ma de lax

Schéma de Lax

  • On remplace ujn par (uj+1n+uj-1n)/2 dans la dérivée temporelle

  • Stabilité ? On trouve

  • On veut |ξ(k)| ≤ 1 pour tout k

  • On appelle cela la condition de Courant.Interprétation : On calcule ujn+1 (en x=xj) à l’aide de uj+1net uj-1n(en x=xj-1 en x=xj+1 ). L’information se propage à la vitesse maximale ±c. Les points uj+1net uj-1ndoivent être en dehors de zone liée à ujn+1.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Sch ma de lax1

Schéma de Lax

  • Dissipation

  • Si |c|Δt < Δx, |ξ(k)| < 1 L’amplitude des modes décroît. Il y a de la viscosité numérique ! Interprétation :

  • 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).

Master 2 AMD - Méthodes numériques pour l'astrophysique


Ordre 2 en temps staggered leapfrog

Ordre 2 en temps : Staggered Leapfrog

  • On va centrer le calcul de ∂u/∂t.

  • Le schéma est du second ordre spatialement et temporellement. Mais on a besoin de l’information à tn et tn-1pour calculer à tn+1

  • Stabilité : On cherche les modes propres. On tombe sur

  • 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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Staggered leapfrog equations du second ordre

Staggered Leapfrog : Equations du second ordre

  • Dans le cas d’une équation de conservation de fluxle leapfrog s’écrit

  • On écrit ce schéma pour l’équation des ondes (r = c∂u/∂x, s = ∂u/∂x).

  • On remplace. La deuxième équation donne

    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

Master 2 AMD - Méthodes numériques pour l'astrophysique


Sch ma en deux temps de lax wendroff

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 flux correspondant

    • On en tire les flux correspondants et on utilise le leapfrog pour tirer ujn+1

  • Dans le cas de l’équation simple (F=cu), l’ensemble se réduit à

Master 2 AMD - Méthodes numériques pour l'astrophysique


Sch ma en deux temps de lax wendroff1

Schéma en deux temps de Lax-Wendroff

  • Stabilité ? On trouve

  • |ξ| ≤ 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 alors que dans le schéma de Lax on a

  • 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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Raffinement d rivation face au vent

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 :

  • 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 !

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations de diffusion paraboliques

Equations de diffusion (paraboliques)

  • 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.

  • Si D est constant, on a l’équation de la chaleuravec une discrétisation immédiate

  • C’est un schéma explicite de type FTCS, mais plus stable…

  • Interprétation : Le pas de temps doit être plus petit que le temps de diffusion à travers une (demi) cellule.

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations de diffusion sch mas implicites

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.

  • C’est un schéma implicite = il faut résoudre un système linéaire (tridiagonal) pour trouver les un+1j(j=0,…,J)+conditions aux limites en j=0 et J-1

  • Stabilité :

  • Mais c’est un schéma du premier ordre  Mieux : Crank-Nicholson !

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations de diffusion sch ma de crank nicholson

Equations de diffusion : schéma de Crank-Nicholson

  • On prend la moyenne du schéma explicite et du schéma implicite

  • C’est encore un schéma implicite. Mais il est du second ordre en temps (centré en tn+1/2)

  • Stabilité :

  • C’est le schéma recommandé pour tous les problèmes de diffusion…

Master 2 AMD - Méthodes numériques pour l'astrophysique


Equations de diffusion plus complexes

Equations de diffusion plus complexes

  • Si D n’est pas une constante… On discrétise DSi D(x), on appelle Dj+1/2=D(xj+1/2). Le schéma FTCS devient

    avec la condition de stabilité

  • Crank-Nicholson devient

  • 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.

Master 2 AMD - Méthodes numériques pour l'astrophysique


  • Login