Introduction aux simulations numщriques.
Contents
1. e Utilisation de l ordre implicit double precision a h o z Il peut tre int ressant d utiliser un ordre implicit none au d part et de pr ciser ensuite le type de chacune des quantit s introduites dans le programme Si vous programmez vous vous rendrez compte tr s vite que si vous adoptez la convention pr c dente un seul fichier in clude en d but de TOUS les sous programmes et TOUTES les fonctions vous pouvez utiliser la convention par d faut du Fortran sans aucun danger et ne vous ferez JA MAIS d erreur de type qui sont l exp rience le montre les erreurs les plus courantes e Parameter nmax 100 les dimensions de TOUS les tableaux de votre programme doivent tre param tris es et la dimension doit absolument se trouver dans l include unique J insiste ne jamais crire quelque chose comme X 100 dans un pro gramme m me s il est vident que vous ne changerez jamais la dimension du tableau X e programme principal if n gt nmax stop nmax dans general data trop petit Ceci est absolument imp ratif Si une donn e d entr e conditionne la taille d un tableau faire un test de ce type e Il faut absolument adopter une politique coh rente pour le passage d une quantit d un sous programme l autre Il y a deux fa ons de faire cela Par common ou par argument Je propose de faire passer par common common qui doit se trouver dans le fichier include et pas ailleurs to2. read logic if logic eq true then do i 0 n 1 2 tab2 i exp tab1 i 1 enddo else do 1i 0 n 1 2 tab2 i 1 expression i enddo endif Noter la structure if proposition logique then else endif d utilisation tr s commune expression i n est pas ici un tableau mais une fonction d finie par l utilisateur Il s agit d un cas particulier de sous programme pour lequel l expression calcul e se r duit un seul nombre C est l quivalent utilisateur des fonctions intrins ques pr d finies vues pr c demment On peut utiliser plusieurs arguments call progi n tabi tab2 res Structure g n rale pour l appel un sous programme V rifier avec beaucoup d attention que le nombre et le type des variables pass s par argument correspondent exactement dans l ordre call et l ordre subroutine C est une source tr s importante d erreur do i 1 n write 1 5f10 5 res i enddo Ecriture sur le fichier 1 r sultat du tableau res Notez le format qui n est pas libre le format libre par d faut est obtenu en remplacant 5f10 5 par un symbole Ce format est dit format fixe 10 caract res pour repr senter un nombre flottant avec 5 chiffres apr s la virgule 4 caract res pour le signe et les chiffres avant la virgule un caract re pour le signe 7 Exemple 100 12345 1234 6789 sont deux chiffres crits avec ce format 5 signifiant qu il y a 5 nombres crits avec ce format par ligne Il e
3. H X r gt lt h k Dans ce qui suit les valeurs propres sont class es par ordre croissant On peut donc crire E HY X E M gt lt Qr k La m thode des puissances repose sur le r sultat vident que lim E H fo gt lt go sous la condition que rA lt 1 pour tout k ce qui est vrai si la constante E est choisie par exemple plus grande que toutes les valeurs propres Notons gt un vecteur quel conque non orthogonal au vecteur propre fondamental on peut donc extraire le vecteur 33 correspondant l tat fondamental par application successive de la matrice E H Application au probl me de l oscillateur anharmonique a Mettre en oeuvre l algorithme qui permet de calculer Y gt E H V gt o gt est un vecteur quelconque initial b En prenant comme vecteur initial un ensemble de coefficients tir s au hasard v rifier que le vecteur Y gt converge vers l tat fondamental Comparer avec les r sultats de la m thode de diagonalisation exacte 2 M thode de Lanczds Cette m thode g n ralise la m thode pr c dente L id e fondamentale est tr s simple Plut t que de prendre comme vecteur d essai pour H le vecteur Y gt E HY Y gt la m thode consiste prendre comme vecteur d essai le meilleur vecteur construit l aide de l ensemble de tous les vecteurs Y gt k 0 n Dit autrement la m thode consiste
4. Ising 2 dimensions 26 V Au voisinage de la transition de phase l algorithme de Swendsen Wang 28 A Algorithme d Swendsen Wang 44 6 34 due da 4e res 4 28 B Simulation du modele d Ising au voisinage de la transition 30 VI M thodes de diagonalisations exactes et m thodes de projection L oscillateur anharmonique quantique 31 FE G c a da t e aN a a E o a E a a E A E A a A i 31 B Oscillateur anharmonique quantique aoaaa 32 1 Diagonalisation exacte L an SNS En Sense 33 2 M thodes de projection m thode des puissances m thode de Lancz s 33 VII Mod le d Heisenberg antiferromagn tique 2D 36 I INTRODUCTION Le but de ce cours est d illustrer l aide de quelques exemples simples comment un ordi nateur peut tre utilis pour tudier un probl me physique complexe n admettant ni solution analytique ni approximation contr l e satisfaisante Les probl mes complexes sont videm ment nombreux en physique particuli rement dans le domaine de la mati re condens e o existe une grande diversit de comportements et de structures D un point de vue technique ce sont le plus souvent des probl mes pour lesquels les m thodes perturbatives perturba tion autour d une solution connue analytiquement ne fonctionnent pas et ou des probl mes correspondant un r gime hautement non lin aire Dans les cas les moins difficiles il s agit simplement de consid rer
5. Quand A est choisi tr s petit le rapport Ri est toujours tr s voisin de l unit et les mouvements d essai sont accept s avec une grande probabilit Malheureusement l espace de configuration est visit avec inefficacit et l estimateur des int grales convergera lentement Dans la limite oppos e il est clair que dans un espace multidimensionnel un mouvement de grande amplitude aboutira quasiment toujours un refus Il y a donc un optimum trouver Cet optimum correspond un choix de A donnant un taux d acceptation fini ni trop petit ni trop grand un taux voisin de 0 5 est souvent conseill Appliquons cet algorithme un cas simple Supposons que l on d sire valuer l int grale suivante 1 g dx exp x gx Cette int grale n a pas de valeur analytique connue pour g diff rent de 0 Nous allons utiliser l algorithme de Metropolis pour valuer cette int grale Evidemment c est la pire des m thodes utiliser pour ce type d int grale mais ceci nous donnera l occasion d exemplifier ce que nous venons dire dans un cas tr s simple On prendra comme densit I x z EXP 2 Proposition de programme Fortran program Metropolis CC Input print delta nsteps g read delta nsteps g CC Configuration initiale tir e au sort xold ranf 0 5 CC Initialisation des accumulateurs valint 0 naccep 0 CC G n ration des configurations successives 25 do k 1 n
6. o O oO O o 2N 2N 1 3N as G G O OR Ame G O G EN O R N N amp O O oO O 0O O P N D N 1 N 1 N 1 N N y y y y 1 2 3 N 1 N avec des conditions p riodiques aux bords comme indiqu sur la figure 1 La fonction de partition du syst me tant donn e par Z 8 D exp 6E 7 S S2 SN o 8 repr sente 1 kT et S 1 ou 1 Ecrire la forme analytique de l nergie moyenne lt E gt a et de la chaleur sp cifique 2sE gt prendre k 1 2 Dans le cas 4 sites N 2 crire un programme qui calcule lt E gt et c en som mant explicitement sur toutes les configurations de spin possibles Dans le cas J 1 verifier la pr sence d un maximum pour en fonction de T faites varier T entre 1 et 3 li l existence d une transition de phase pour ce syst me 3 R fl chir au probl me de la programmation pour une valeur de N arbitraire pas v ident exposer le probl me proposer et programmer ventuellement une solution et au temps de calcul associ quelle est la d pendance en N 27 4 Mettre en oeuvre l algorithme de Metropolis pour un syst me de taille arbitraire N Pour cela criver un programme qui visite les sites les uns apr s les autres tirer au hasard la valeur 1 ou 1 avec g n rateur de nombres al atoires fourni pour le spin courant puis accepter ou rejeter avec le crit re de Metropolis et faire cela un grand nombre de fois D
7. c dente M ibset M i Le i me bit de M est mis 1 et l entier r sultant est M M ibclr M i Le i me bit de M est mis z ro 2 Appliquer la m thode de puissance pour calculer l tat fondamental du syst me Regarder la d pendance en fonction du nombre de sites V rifier que l tat fondamental est un singulet c est dire est tat propre de Evaluer galement l aimantation altern e d finie par M N lt E s 57 10 gt iEA ieB o 0 gt d note l tat fondamental Discuter les r sultats 3 Calculer l tat fondamental par m thode de Lancz s 4 Calcul de la fonction de Green de spin Montrer que la transform e de Fourier de la fonction de Green peut s exprimer comme une valeur moyenne de l inverse de H Exprimer cette fonction de Green sous forme de fraction continue Discussion physique des r sultats 38
8. h i j fnconnect i j gt j 1 o i gt repr sente un l ment quelconque de la base nhstate i repr sente le nombre total d tats de la base connect s l tat i gt nconnect i j donne le num ro du ji me tat connect i gt et h j est l l ment de matrice non nul correspondant De cette mani re seuls les l ments de matrice non nuls de la matrice Hamiltonienne sont stock s Soit gt un vecteur arbitraire de composantes c dans la base on a nsttot nhstate i HIY gt amp h i j jnconnect i j gt iS j 1 o nsttot repr sente la dimension de la base Afin de faciliter la programmation on pourra repr senter un tat quelconque de la base l aide d un entier en base 2 M m My n32 n22 naD o n 0 si par exemple le site i est occup par un spin down et n 1 si le site est occup par un site up Ceci permet de boucler sur tous les tats en parcourant les entiers de 0 2 1 1 A partir d un entier donn on peut remonter l occupation des sites tr s 37 simplement en utilisant les fonctions de manipulation de bits du Fortran On a les trois fonctions l mentaires logic btest M i logic vaut true si le i me bit de l entier M vaut 1 sinon logic false Sur la machine que vous utilisez les entiers sont par d faut repr sent s sur 32 bits et le bit 0 est le bit de poids faible bit 1 dans la repr sentation pr
9. l intervalle i n sont entiers Toute d claration explicite du type d un mot annule cette convention pour le mot Noter l existence de la d claration implicit none qui limine cette convention Si vous utilisez cette convention vous devez d clarer le type de toutes les variables ou con stantes utilis es dans le programme C zone des parameters parameter nmax 100 pi 3 1415926535898 deuxpi 2 d0 pi D claration permettant de fixer quelques valeurs constantes Particuli rement int ressante pour d finir les dimensions des tableaux c d claration des tableaux dimension tabi nmax tab2 0 nmax D claration des tableaux tabl nmax signifie que vous r servez en m moire centrale nmax l ments de m moire tab1 1 tab1 2 tabl nmax 100 Plus g n ralement vous pouvez crire tab2 i j pour tab2 i tab2 i 1 tab2 j i j entiers positifs nuls ou n gatifs zone communes common c1i termi term2 D finition d une zone de m moire centrale commune plusieurs sous programmes cl est le nom de la zone commune terml term2 sont deux entiers ils ont t d clar s explicitement comme tels au d but de ce programme qui pourront tre utilis s dans les sous programmes o la zone commune a t introduite Il est imp ratif de conserver l ordre le type il faut absolument red finir le type dans les sous programmes et le nom des variables constantes ou tableaux pass s par common open 1 file re
10. la matrice H lt 6 H 6 gt 32 1 Diagonalisation exacte 1 En utilisant le formalisme des op rateurs de cr ation et d annihilation calculer H On rappelle les d finitions suivantes Op rateur de cr ation a lc 2 op rateur d annihilation a AU avec commutateur a a 1 On peut montrer assez simplement que 1 r Ya z at polz Montrer que i et en it rant cette relation trouver l expression explicite des l ments de matrice de H 2 En utilisant un programme de diagonalisation que l on vous fournira par exemple pro gramme JACOBI calculer les valeurs propres et les vecteurs propres de la matrice obtenue en tronquant la base une valeur maximum N V rifier que les niveaux d nergie d croissent en fonction de la dimension de la base utilis e et convergent partir d une certaine taille D terminer la d pendance de l effort num rique en fonction de N 3 Tracer la fonction d onde fondamentale pour diff rentes valeurs de g en utilisant les rela tions suivantes pour valuer les polyn mes d Hermite Vnbn 1 x vVn lp 41 x Ho 1 H 2r Hiel 22H fe 2n Hn 1 8 Commenter les r sultats 2 M thodes de projection m thode des puissances m thode de Lancz s 1 M thode des puissances Power method Nous voulons d terminer l tat fondamental du syst me Toute matrice r elle sym trique peut s crire sous la forme d composition spectrale
11. t d ergodicit est tr s importante car elle exprime le fait que tout tat de l espace de configuration peut tre atteint en un nombre de pas fini tous les tats sont visit s et en fait cette condition est n cessaire dans la d monstration pr c dente afin viter le cas pathologique que nous n avons pas consid r o la probabilit de transition vers un tat donn j s annule quelque soit l tat i On peut d montrer le th or me suivant qui pr cise la notion de convergence Soit f une distribution c est dire un ensemble de N nombres positifs On va crire l application de la probabilit de transition cette distribution sous la forme pre yn Paz Pfi9 J On a le r sultat suivant lim ju Si Po I vf Les diff rentes tapes de la d monstrations sont les suivantes e On associe P _ une matrice r elle sym trique d finie de la mani re suivante 1 Mi Vi yH Notons que la probabilit de transition elle m me ne d finit pas une matrice sym trique 23 e On v rifie facilement que VII est tat propre de M avec valeur propre 1 S MAIL Il j e On v rifie aussi que fo VII e On utilise la d compostion spectrale de M qui nous dit que dans la limite des grandes valeurs de n la matrice se r duit au projecteur sur le sous espace des vecteurs propres associ la plus grande valeur propre On peut montrer que la matrice M du fait de sa structure a de
12. tre obtenues par de simples quadratures On a v t T t F a t m F x X L nergie est videmment une int grale premi re du mouvement dE dt d o dx 2 u zE V x et on a donc z dx Hist LES ce qui permet d valuer la trajectoire x t b Nombre de degr s de libert plus grand que un D gt 1 L int grabilit correspond l existence de D 1 int grales du mouvement en plus de l nergie totale M thodes des sections de Poincar Cette m thode permet de d celer l existence d int grales premi res du mouvement sup pl mentaires Elle consiste visualiser des coupes de l hypersurface nergie constante Nous allons illustrer cette id e dans un cas deux degr s de libert Consid rons par ex emple la section de Poincar correspondant vs 0 avec v gt 0 choix du sens du temps S il existe une int grale premi re I en plus de l nergie on a I x y Vx Vy cte C est dire que l on peut crire E x y 0 v cte et I x y 0 vy cte ce qui implique qu il existe une relation implicite f x y 0 Donc dans le cas o le syst me une int grale suppl mentaire les sections de Poincar doivent tre des courbes continues ventuellement non connexes dans le plan x y D existence de sections de Poincar discontinues l allure chaotique s oppose donc l existence d une int grale premi re suppl mentaire Exemple de program
13. 2D La constante d change J est posi tive On a 1 ma l x y ou z avec les trois matrices de Pauli AUTS PO a ENE DA ETA Gi D finissons les op rateurs suivants 00 o Ge QU ST S LE De On v rifie facilement leurs propri t s fondamentales st1T gt 0 SIT lE S el S 1 0 EE e TS nes o gt 0 et gt 1 Une contribution l mentaire pour un lien peut se r crire ce 4 i5j zi SF 57 S7 575 36 c est dire re 1 Si gt 7l gt 33 1 S Si H gt zl W gt e 1 1 SiSi N gt 51 41 gt 7l N gt 3 1 1 si5i N gt 5I N gt IS On se placera dans la base de dimension 2 d finie par Im1 mn gt avec m f ou mi On va mettre en oeuvre un programme pour valuer l nergie du fondamental ainsi qu un certain nombre de ses propri t s Il est important de remarquer que de s v res limitations de m moire centrale apparaissent tr s vite dans ce type de probl me on a 27 tats con sid rer Il faut donc faire tr s attention ne pas utiliser des tableaux inutilement grands c est l objet de la premi re question 1 L tape fondamental de toute m thode de projection est la partie qui calcule l application de la matrice Hamiltonienne sur un vecteur donn de l espace vectoriel Con struire les trois tableaux suivants qu on pourra appeler par exemple h i j nhstate i et nconnect i j avec la signification suivante nhstate i Hji gt
14. HAL archives ouvertes Introduction aux simulations num riques Michel Caffarel To cite this version Michel Caffarel Introduction aux simulations num riques DEA Cours enseign en 1995 et 1996 au DEA Champs Particules Mati res Paris 6 7 et 11 38 pages 2006 lt cel 00092932 gt HAL Id cel 00092932 https cel archives ouvertes fr cel 00092932 Submitted on 12 Sep 2006 HAL is a multi disciplinary open access archive for the deposit and dissemination of sci entific research documents whether they are pub lished or not The documents may come from teaching and research institutions in France or abroad or from public or private research centers L archive ouverte pluridisciplinaire HAL est destin e au d p t et la diffusion de documents scientifiques de niveau recherche publi s ou non manant des tablissements d enseignement et de recherche fran ais ou trangers des laboratoires publics ou priv s Introduction aux simulations num riques Michel Caffarel Laboratoire de Chimie Th orique Tour 22 23 ler tage 4 place Jussieu 75252 Paris Cedex 05 e mail mc lct jussieu fr Avril 1996 Sommaire I Introduction 2 II Aspects pratiques des simulations 3 A Organisation physique d un ordinateur 3 1 Repr sentation des nombres ts ins ns de Ne Pepe se pas 3 2 Vitesse finie du PrOC SS NL 4 404 AUS LS RAA EE Sa Re x se 4 3 Capacite d s OPA serie
15. al et initial qui vous permet de visualiser le texte de positionner le curseur 10 o vous voulez et d effectuer un certain nombre de commandes 2 le mode dit d insertion dans lequel on entre en tapant la lettre 1 Vous pouvez alors entrer du texte partir de l endroit o se trouve le curseur Pour sortir du mode insertion il faut appuyer sur la touche Echappement Escape E Quelques conseils de programmation L exp rience montre qu il est tr s important d avoir une discipline stricte de program mation Ceci permet de gagner un temps norme et d viter de longues heures et jours de recherche de bugs qui peuvent facilement tre vit s Evidemment la bonne strat gie employer est une question tr s personnelle et en g n ral chacun est persuad d avoir la meilleure m thode n argumentez pas a ne sert rien Essayer de d velopper tr s vite votre propre strat gie et copier que ce qui vous convient chez les autres Quelques conseils faites en ce que vous voulez Il faut programmer vite L objectif n est pas de programmer la perfection et de faire la simulation ultime mais d tudier un probl me physique et d utiliser l ordinateur comme un outil et non pas comme une fin en soi Plus vous programmerez vite plus vous serez disponible pour le probl me tudier et plus vous serez psychologiquement pr ts utiliser sans h siter l ordinateur pour attaquer un probl me nou
16. ans le cas N 2 v rifier que vous retrouvez les r sultats exacts obtenus la question 2 5 Dans tout programme Monte Carlo il est IMPERATIF de calculer une estimation de l erreur sur les moyennes des quantit s calcul es estimation des fluctuations statistiques Soit Q la quantit calcul e ici Q lt E gt c estimer l erreur 0Q en utilisant la formule suivante 1 1 JMM I 2 Qi lt Q gt o M repr sente M calculs ind pendants donnant Q et lt Q gt est la moyenne des Q pour ces M calculs prendre M 10 6 Augmenter N et tudier l volution des diff rentes quantit s avec la taille sQ Remarque importante Il est fortement conseill d introduire un tableau de connection entre sites afin de faciliter la programmation de ce probl me Pour cela on introduit un tableau nvois i k i 1 N k 1 4 avec la signification suivante nvois i l donne le num ro du site de gauche pour le site num ro i nvois i 2 le num ro du site d en haut nvois i 3 de droite et nvois i 4 d en bas De cette facon les sommes doubles sur i j peuvent tre faites simplement en sommant sur les sites i et l indice k associ V AU VOISINAGE DE LA TRANSITION DE PHASE L ALGORITHME DE SWENDSEN WANG A Algorithme de Swendsen Wang L algorithme de Metropolis utilis dans la section pr c dente repose sur un mouvement d essai local la nouvelle configuration d essai est construite en tentant de modifier la valeur du sp
17. babilit non nulle apr s au plus N pas de connecter deux configurations quelconques de l espace de phase Notre deuxi me point va consister d montrer qu un tel mouvement d essai v rifie le bilan d taill Tout d abord explicitons en terme de probabilit s le mouvement d essai propos On cherche calculer la probabilt de transition P S S connectant deux configuration quelconques Il est clair que cette probabilit est non nulle que dans le cas o S et S ne diff rent que par un flip des spins sur un cluster C Pour une configuration donn e S il existe un ensemble fini de clusters susceptibles d tre retourn s N clusters un site un certain nombre de clusters de taille 2 etc je noterai cet ensemble Cs l indice S pour rappeler que cet ensemble d pend de la configuration Il est facile de se convaincre que la probabilit de transition d essai s crit IT p l Si Si l Il a p Si 5 P S S lt i jJ gt EC lt 1 9 gt us D RE CE CeCs lt i j gt ec lt i j gt 8C La notation lt i j gt C signifiant que la liaison i j appartient au cluster C diff renciant les deux configurations et S et lt i j gt OC signifiant qu il s agit d une liaison du bord de cluster qui connecte un site i du cluster un site j qui n appartient pas au cluster Montrons 29 maintenant que cette probabilit de transition d essai v rifie le bilan d taill Pour cela
18. cul exact En revanche d s que des r els sont utilis s une erreur dite de pr cision finie est introduite Cette erreur peut toujours tre r duite en repr sentant les r els l aide de plusieurs mots Par exemple utiliser 2 mots s appelle travailler en double precision 4 mots en quadruple etc On vite g n ralement d utiliser des pr cisions trop grandes parce que les temps de calcul augmentent de mani re importante En fait la double pr cision sur huit octets sur des machines destin es au calcul num rique est tr s souvent cabl e c est dire que le processeur est physiquement pr par travailler en double pr cision et les temps de r ponse ne sont pas significativement plus grands qu en simple pr cision Au del il s agit d une solution logicielle c est dire qu un programme a t crit pour juxtaposer des mots ensemble Cette solution est donc beaucoup plus lente Il est important de souligner qu on peut toujours soi m me juxtaposer autant de mots qu on veut pour repr senter un r el et atteindre ainsi une pr cision arbitraire notons que ceci est dej fait pour vous dans les programmes de calcul symbolique comme Maple ou Mathematica Ceci peut tre tr s utile dans des applications o une tr s grande pr cision est essentielle ou pour v rifier la stabilit num rique d un algorithme Ces notions sont illustr es dans les exercices pr sent s dans la section I F 2 V
19. diagonaliser H dans la base Y gt H V gt H Y gt H Y gt L espace vectoriel sous tendu par cet ensemble est quelque fois appel l espace de Krylov associ Y et H On va voir que la convergence dans un tel espace est norm ment augment e D un point de vue pratique il est souhaitable de disposer d une base orthonormale on va donc effectuer une orthonormalisation syst matique des vecteurs construits Partant de Y gt notre vecteur arbitraire de depart on construit le premier vecteur de Lanczds normalis Y gt V lt V IY gt L tape suivante consiste calculer H uo gt On construit ensuite le second vecteur de uo gt Lancz s u1 gt l aide de ce vecteur et du vecteur de Lancz s pr c dent uo gt Ba u1 gt H uo gt uo gt a est d termin en imposant que u gt soit orthogonal uo gt et B2 d termin pour que le vecteur u1 gt soit de norme unit ai lt uo H uo gt Pa y lt uo H uo E uol H uo gt A lit ration suivante on crit le nouveau vecteur sous la forme Ba lu gt H u gt a u1 gt b uo gt Noter que par construction le nouveau vecteur est orthogonal uo gt a est obtenu en imposant lorthogonalit u gt et 83 pour que le vecteur soit de norme unit as lt u H u gt Finalement on arrive au cas g n ral Biri ui gt 5 H ui 1 gt Glisse 1 1 N 34 avec ai l
20. e s crit dans la base de tous les vecteurs Lanczds pr c dents on peut donc arr ter le processus it ratif l espace de Krylov obtenu est ferm La convergence en fonction du nombre d it rations est au moins exponentiellement rapide Un argument heuristique peut nous en convaincre une it ration n l tat fonda mental obtenu peut s crire sous la forme Vo X cH Y gt i 0 L nergie associ e cette fonction d onde est optimale principe variationnel par rapport au coefficient c On a donc E Yo lt E 35 o pi i 0 _ 1 gt H Y gt 1 Pour n assez grand on voit donc que E Yo est major e par une quantit qui converge exponentiellement vite vers l nergie de l tat fondamental Probl me 1 En partant du m me vecteur initial que pr c demment construire la matrice tridiag onale de Lanczos 2 Diagonaliser cette matrice avec la m thode de Jacobi Etudier la convergence des nergies en fonction de n et discuter l apparition d instabilit associ es ce qu on appelle souvent des solutions fant mes On fera cette tude en simple double et quadruple pr cision VII MOD LE D HEISENBERG ANTIFERROMAGN TIQUE 2D L Hamiltonien est d fini de la mani re suivante H N JSS lt i 3 gt sur un r seau bidimensionnel carr avec conditions aux limites p riodiques identique au r seau utilis pr c demment pour le mod le d Ising
21. er les convergences en fonction du temps calcul obtenue avec l algorithme local et l algorithme de cluster VI M THODES DE DIAGONALISATIONS EXACTES ET M THODES DE PROJECTION L OSCILLATEUR ANHARMONIQUE QUANTIQUE A Th orie Le probl me de la diagonalisation de matrices est pr sent dans la plupart des simula tions des syst mes quantiques On sait que l tat d un syst me quantique est d crit par un vecteur d un espace vectoriel espace des tats Dans le cas d un tat pur ce vecteur est une combinaison d tats propres de l op rateur Hamiltonien Dans le cas d un m lange statistique m canique statistique quantique les valeurs moyennes des grandeurs physiques sont repr sent es par un m lange pond r des valeurs moyennes associ es aux tats propres du syst me La d termination des tats propres de l op rateur Hamiltonien est donc au centre de toute tude quantique Il est important de souligner qu en derni re analyse la m canique quantique d un syst me complexe peut toujours se ramener la r solution d un probl me d alg bre lin aire l mentaire dans un espace vectoriel de dimension N finie la seule difficult pratique tant bien s r d valuer la limite N Soit M une matrice NxN On cherche r soudre le probl me aux valeurs propres suivant Mu hu k 1 N On va supposer par simplicit que les l ments de matrice M sont r els g n ralement c est le cas da
22. etc I x est une densit de probabilit c est dire Il x gt 0 et de He 1 L algorithme le plus utilis pour calculer de telles int grales est l algorithme de Metropolis Cet algorithme est un algorithme probabiliste qui permet d engendrer des configurations x selon la densit I x Ainsi l int grale I peut tre valu e sous la forme P I lim LS F0 Po P par En pratique c p y ls f x our P assez grand p 2 a TE g o z p est un ensemble de P configurations construites avec l algorithme de Metropo lis et o F repr sente l erreur r siduelle th or me de la limite centrale loi des grands nombres On va voir qu une caract ristique essentielle de l algorithme est qu il ne requiert que la connaissance de probabilit s relatives c est dire Ie Ceci est particuli rement int ressant car les int grales valuer mettent souvent en jeu une probabilit de la forme exp x ie SPELE Fdzexp o a 21 o x est connue explicitement la norme f dx exp x tant inconnue int grale mul tidimensionnelle compliqu e Exemple ensemble canonique classique o 5H H tant la fonction nergie totale du syst me B Algorithme de Metropolis Sans perdre de g n ralit nous nous placerons dans le cas d un espace de configuration ayant un nombre fini N de configurations possibles e D f 1 Probabilit IL gt 0 i 1 N e D f 2 Probabili
23. fa on Une r gle tr s importante est que le temps que vous passez v rifier votre programme doit tre beaucoup plus long que le temps que vous avez pass l crire 11 Il faut programmer de mani re syst matique Choisissez la m thode que vous voulez mais soyez syst matique Par exemple Programme Fortran program general include general data dimension x nmax etc print n7 end CC Ici vous donnez des explications d taill es de ce qui est fait CC dans le sous programme progi Il faut absolument donner CC TOUS les d tails formule exacte conventions multiples r f rences CC dans la litt rature etc CC sinon a ne sert absolument rien cc CC input CC kkk CC liste des quantit s utilis es comme input et non modifi es cc CC output CC kko k CC liste des quantit s calcul es et modifi es cc subroutine progi include general data Fichier general data ins r par include dans TOUTES les subroutines et fonctions implicit double precision a h o 2 Parameter pour dimensions de TOUS les tableaux du programme parameter nmax 100 common c1i x1 nmax 12 Commentaires e L ordre include permet d introduire le fichier ayant le nom general data l endroit o l ordre include se trouve Il faut absolument introduire cet include dans TOUS les sous programmes et fonctions du programme Ceci vous vitera un nombre incalculables d erreurs
24. il faut calculer le rapport IL PUS 55 II 1 p s S PS gt S Lire O l aisee O O PSS M psi s plsi S lt i j gt EC lt i j gt EC Ceci tant vrai parce que les normalisations sont identiques dans les deux cas Apr s sim plification on obtient P S 5 1 pl Si 5 PSS Il TE lt i j gt EC exp JO Min 0 29J6 5 Min 0 28J6 5 lt i j gt EC exp D 28JS S lt i j gt EC la derni re galit r sultant de la relation Min 0 X Min 0 X X Par ailleurs on a exp BJ X SiS P S lt i j gt P S exp 8J X Si S lt 1 3 gt exp 8J D SiS BJ D S lt i j gt EC lt i j gt EC exp J 28 55 lt i J gt E3C Le bilan d taill est donc v rifi B Simulation du mod le d Ising au voisinage de la transition 1 En utilisant le programme bas sur l algorithme de Metropolis local observez le ph nom ne de ralentissement critique Pour cela on regardera la convergence en fonction du nombre de pas Monte Carlo de l estimateur de l nergie moyenne et de la chaleur sp ci fique au point critique 2 Mettre en oeuvre l algorithme de Swendsen Wang On programmera les estimateurs 30 de l nergie moyenne de la chaleur sp cifique ainsi que de la taille moyenne des agr gats form s une temp rature donn e 3 Analyser la d pendance de la taille des agr gats en fonction de la temp rature 4 Au point critique compar
25. important de souligner que lorsque le mouvement est refus l ancienne configuration doit bien tre consid r e comme une nouvelle configuration en particulier sa contribution aux int grands doit tre ajout e comme toute autre nouvelle configuration Algorithme de Metropolis classique J appelle ici algorithme de Metropolis classique l algorithme qui est g n ralement pr sent dans les ouvrages et qui correspond en fait un choix particulier de probabilit de transition d essai P i j Constante 24 pour un sous ensemble d tats j voisins dans l espace de configurations notion de voisinage i d pend du probl me Pour cet algorithme on a donc R T Cet algorithme est g n ralement exprim pour un espace de configuration continu No tons 1 2 Z4 un l ment de l espace des tats la probabilit de transition s crit formellement TRY Ad o O1 est la fonction caract ristique de l hypercube centr en 0 et de cot de longueur A En pratique cela veut dire que les mouvements d essai sont effectu s de la mani re suivante yi z Al amp 0 5 i 1 d o sont des nombres al atoires uniformes tir s entre 0 et 1 A est une quantit posi tive a priori arbitraire comme le choix de toute probabilit de transition d essai Cepen dant A d termine la vitesse de convergence vers la densit stationnaire Prenons les deux cas limites possibles pour A
26. in en un seul site Or on sait qu il existe une temp rature critique dans ce syst me correspondant une transition de phase du second ordre Au point critique la longueur de correlation d finie par la relation suivante lt Si lt Si gt S lt S gt gt exp 1 7 grand o lt gt d signe la valeur moyenne thermodynamique correspondant l ensemble canon ique devient infinie Physiquement cela signifie qu il existe des corr lations tr s longue port e dans le syst me et que des agr gats de spins identiques et de tailles arbitraires se forment Il est clair qu au point critique les valeurs moyennes obtenues partir d un al gorithme qui visite l espace des tats par des modifications locales est vou l chec Ce ph nom ne est connu sous le nom de ralentissement critique critical slowing down 28 Ce qu il faut c est imaginer des mouvements d essai qui incorporent cet l ment physique fondamental du syst me dans son r gime critique Nous allons pr senter une telle solution Il s agit d une variante du travail original de Swendsen et Wang R H Swendsen and J S Wang Phys Rev Lett 58 86 1987 Une tape l mentaire de cet algorithme correspond la construction d un agr gat ou cluster de N spins identiques dans le syst me 1 lt N lt N N tant le nombre de sites total e Tirage d un site quelconque du syst me Ce si
27. it il est important de souligner qu il n est pas n cessaire de conna tre grand chose sur l organisation physique des ordinateurs pour les utiliser de mani re satisfaisante n oublions pas que nous ne visons pas la simulation ultime Le seul point tr s important r aliser est ce qu on peut appeler les aspects finis de l ordinateur Un ordinateur est compos d un processeur qui travaille une vitesse finie l aide d une repr sentation finie des objets qu il manipule entiers r els etc Il est associ une m moire centrale qui a une capacit de stockage finie et laquelle il acc de en temps fini Un processeur peut aussi acc der un disque dur m moire dure avec videmment l aussi une capacit et un temps d acc s finis bien sup rieurs d ailleurs Ces aspects tant tr s importants ils seront illustr s titre d exercices sur des exemples concrets dans la section ILF 1 Repr sentation des nombres On appelle g n ralement mot l entit de base manipul par le processeur Pour la plupart des stations de travail le mot est compos de 4 octets c est dire 32 bits le bit est l l ment l mentaire ultime et il peut prendre que les valeurs 0 ou 1 Les entiers sont g n ralement par d faut repr sent s par un mot et les r els galement Dans ce cas on parle de calculs en r els simple pr cision Un calcul effectu l aide d entiers uniquement est un cal
28. itesse finie du processeur Les machines modernes op rent avec une puissance qui leur permet d effectuer de l ordre de quelques millions de multiplications r elles la seconde 3 Capacit de stockage Typiquement sur les stations de travail la m moire centrale peut stocker quelques dizaines de millions de mots Les plus grandes m moires peuvent atteindre quelques milliards de mots La capacit de stockage sur disque dur est virtuellement infinie on peut juxtaposer les disques Le vrai probl me n est pas la capacit des disques mais le temps d acc s l information stock e voir discussion qui suit 4 Temps d acc s finis l information en m moire Une caract ristique tr s importante des m moires est le temps qu il faut pour acc der un mot stock Les temps d acc s la m moire centrale sont beaucoup plus courts que ceux associ s la m moire dure Dans l exercice de la section IL F 3 on valuera tr s grossi rement cette diff rence Un facteur cent est un ordre de grandeur typique B Quelques l ments de Fortran Le Fortran est le langage de programmation de la communaut scientifique interna tionale M me si ce n est pas et de beaucoup le meilleur langage de programmation c est celui que vous retrouverez partout dans les laboratoires de recherche dans les banques de programmes etc Il est donc indispensable de l apprendre C est un langage tr s simple et quelques heures avec
29. l nergie E variant de 1 2 Observer la modification brutale de la section en fonction de l nergie correspondant au passage d un r gime quasi int grable un r gime de chaos fort Expliquer comment on peut comprendre ce ph nom ne en raisonnant sur le nombre d int grales premi res du mouvement 7 Refaire la m me tude pour un syst me de deux oscillateurs coupl s d crit par le poten tiel her eds 3 31 322 i Vey 58 5Y ty y aey 5EY 4 20 Discuter les r sultats que vous obtenez IT On se propose maintenant d tudier un syst me non conservatif repr sentant un os cillateur entretenu et d crit par l quation dite de van der Pol do d0 EERE 8A dE e 0 Le 0 5 1 Ecrire un programme qui int gre les quations du mouvement 2 Montrer qu il existe un attracteur cycle limite pour les trajectoires dans l espace des phases 0 IV M THODES DE MONTE CARLO TRANSITION DE PHASE DANS LE MOD LE D ISING 2 D A M thodes de Monte Carlo Le coeur des m thodes de Monte Carlo repose sur l valuation num rique d int grales multidimensionnelles de la forme 1 f dall f x x espace de configuration ou espace des tats Le signe int grale repr sente ici aussi bien une int grale qu une somme discr te Donnons quelques exemples E R RN o N est un nombre de particules E 2 o N est un nombre de sites pour un syst me deux tats par site
30. lations honorables sans beaucoup d efforts Il faut souligner que les simulations correspondant ce qu on appelle l tat de l art state of the art ne sont g n ralement que des simulations dont les algorithmes ont t pouss s leurs limites ul times et ne repr sentent malheureusement bien souvent qu un travail fastidieux d analyse num rique n apportant pas grand chose de plus du point de vue de la connaisance physique du probl me Pour donner un exemple pr cis crire un programme qui permet de simuler un mod le de Hubbard pour un r seau de 12 sites le mod le de Hubbard est un des mod les f tiches de la communaut des fermions fortement corr l s peut se faire en quelques jours atteindre les 16 sites n cessite certainement quelques longues semaines mois de travail la limite thermodynamique correspondant videmment un nombre de sites de l ordre de 10 Un autre point important est de bien prendre en compte les processus physiques sous jacents lors de la simulation Ce principe est la base de toute simulation efficace Chaque probl me est cependant tr s sp cifique de ce point de vue et seule la mise en oeuvre de cette id e sur des probl mes concrets permet de s en convaincre Le cas du mod le d Ising trait avec l algorithme de Swendsen Wang en est une tr s bonne illustration voir Section V II ASPECTS PRATIQUES DES SIMULATIONS A Organisation physique d un ordinateur En fa
31. m moire centrale et un l ment de m moire dure Premier programme program tempsi print n7 read n open 1 file ecriture rewindi do i i n 15 do j 1 10000 write 1 j enddo rewindi enddo end Deuxi me programme program temps2 parameter nmax 10000 dimension x nmax print n7 read n do i i n do j 1 10000 x j j enddo enddo end III DYNAMIQUE MOL CULAIRE CHAOS CLASSIQUE A G n ralit s 1 Equations diff rentielles Une quation diff rentielle g n rale peut toujours se r duire un syst me de N quations diff rentielles du 1 ordre coupl es ww f t yn yN i 1 N plus des conditions aux limites qui jouent un r le important On peut distinguer deux classes de conditions aux limites i Conditions aux valeurs initiales Exemple les valeurs de 41 0 yn 0 sont impos es ii Conditions aux limites Exemple les valeurs de 41 0 yx 0 yr 1 t yn t sont impos s Il existe de nombreuses m thodes pour int grer le syst me pr c dent Runge Kutta ex trapolation de Richardson etc voir Numerical Recipes Cependant toutes les m thodes sont construites partir de la m me id e l mentaire Illustrons cette id e avec le syst me le plus simple qu on puisse imaginer On veut r soudre l quation y t J y t 16 On utilise une repr sentation discr tis e de la d riv e ainsi que l quation donnant l expression de la d
32. me Fortran 18 program dynmol CC Input print Tau Nsteps read tau nsteps print x0 y0 vx0 vy0 read xold yold vxold vyold CC Calcul de l nergie initiale e0 0 5 vxold k 2 vyold 2 v xold yold print Energie e0 CC Calcul de l acc l ration initiale call computea xold yold axold ayold CC Construction it rative de la trajectoire do k 1i nsteps xnew xoldttau vxold 0 5 tau x2 axold ynew yoldttauxvyold 0 5 tau x2 ayold call computea xnew ynew axnew aynew vxnew vxold 0 5 tau x axold axnew vynew vyold 0 5 taux ayold aynew e0 0 5 vxnew 2 vynew 2 v xnew ynew print t x y energie print float k xtau xnew ynew e0 xold xnew yold ynew vxold vxnew vyold vynew axold axnew ayold aynew enddo print Simulation terminee end subroutine computea x y ax ay ax x 2 y k 2 x ay y 2 x k 2 y end function v x y v 0 5 x k 2 ykk2 x k2 y k2 end B Apparition du chaos dans un syst me d oscillateurs coupl s I On se propose d observer le d veloppement du chaos pour un syst me de deux oscilla teurs coupl s d crit par le potentiel Ces V r y se a 1 19 l nergie tant donn e par E ei vi V x y masses 1 2 1 Ecrivez un programme qui met en oeuvre l algorithme de Verlet pour int grer les quations classiques du mouvement on travaillera en double pr cision On rappelle que l algorithme s crit SAND Le al S z NE T vlt
33. n 2 avec zo 1l et x 1 6 On peut montrer facilement que la solution est En 1 6 Ecrire un programme qui calcule 1 It rativement la solution de l quation pr c dente 2 Directement la solution connue 3 Analyser l volution de l erreur en fonction de n Faites cette tude en simple pr cision en double pr cision et en quadruple pr cision la quadruple pr cision peut tre obtenue sur votre machine IBM en compilant le programme crit en double precision l aide de l option de compilation suivante xlf qautodbl dblpad programme f Une solution possible program simpleprecision 14 print n read n if n le 1 stop n plus grand que 1 SVP x0 1 x1 1 6 do i 1 n 1 x2 37 x1 6 x0 x0 x1 x1 x2 enddo print n n print xn calcule iterativement x2 print xn solution exacte 1 6 n end program doubleprecision implicit double precision a h o z print n read n if n le 1 stop n plus grand que 1 SVP x0 1 d0 x1 1 d0 6 d0 do i i n i x2 37 d0 x1 6 d0 x0 x0 x1 x1 x2 enddo print n n print xn calcule iterativement x2 print xn solution exacte 1 4d0 6 d0 n end 3 Temps d acc s finis la m moire centrale et la m moire dure En faisant tourner les deux programmes suivants et en utilisant la commande time valuer grossi rement le rapport entre le temps d acc s moyen un l ment de
34. ns les applications On se placera aussi dans le cas o la matrice est sym trique cas physique standard pour une matrice repr sentant une observable physique La matrice M est donc hermitique ou auto adjointe et toutes ses valeurs propres sont r elles Il existe de nombreuses m thodes pour diagonaliser une matrice Elles d pendent du type de matrice r elle sym trique r elle asym trique tridiagonale matrice avec une structure de bande un petit nombre d l ments non nuls au dessus et en dessous de la diagonale principale matrice pleine etc et du fait que vous vouliez ou non tout ou partie des valeurs propres et ou des vecteurs propres On entre ici tr s rapidement dans une situation critique d un point de vue de la simu lation num rique Jusqu maintenant dynamique mol culaire Monte Carlo la quantit d information stocker en m moire centrale ne posait aucun probl me Par exemple dans le cas de la dynamique mol culaire si on note N le nombre de particules et D la dimension de l espace il faut stocker 2ND coordonn es c est dire qu il faut essentiellement stocker 31 2ND mots en m moire centrale En g n ral l essentiel de la physique correspondant la limite thermodynamique N oo est atteinte pour des valeurs de N telles que 2ND lt lt taille de la m moire centrale Dans le cas de la dynamique mol culaire l effort num rique porte essentiellement sur le temps d ex cution
35. qui est proportionnel au nombre de pas l mentaires effectu s Topu Te No Dans le cas des probl mes de diagonalisations exactes la limite de la possibilit de stockage de la machine est tr s vite atteinte Typiquement deux situations se pr sentent e M thodes de diagonalisations exactes calcul complet de toutes les valeurs propres et vecteurs propres de la matrice Effort num rique M moire centrale N Topu N e M thodes de projection Lanczds Davidson Olsen etc calcul de un ou quelques tats proches de l tat fondamental Effort num rique Topv 100N M moire centrale 3N pour des matrices creuses hamiltoniens locaux Illustrons ces deux types d approche sur un probl me quantique simple quoique non soluble B Oscillateur anharmonique quantique On va se proposer de r soudre le probl me d une particule dans un potentiel anhar monique une dimension ressort anharmonique quantique Il s agit de r soudre le prob l me suivant Hg Eigi avec f p7 1 H pre L 2 4 Data Die Pour cela on crit les solutions sur une base compl te discr te de fonctions l mentaires N 00 N x lim X criyila o l ensemble des fonctions 4 est solution du probl me au valeurs propres de l oscillateur harmonique et forme donc un ensemble complet 2 E a e de pale DH a o H x sont les polyn mes d Hermite Pour r soudre notre probl me il faut trouver le spectre de
36. r t zH a t 7 3 e hs AE V V o r x y U Vz vy d Ss Fe et T le pas temporel choisi 2 Construire une trajectoire en prenant des conditions initiales x y v al a toires vfo tant choisi de mani re imposer la valeur de l nergie E d sir e v XE VO y0 jui 3 Analyser le comportement num rique de l nergie en fonction du nombre de pas et du pas temporel Ecrire un syst me d quations pour les quations du mouvement admet tant une erreur d int gration un ordre de grandeur plus petit en T Retrouver ce r sultat num riquement on illustrera ces resultats l aide de courbes Quelle est la valeur de 7 optimale pour un calcul en double precision 4 Tracer une trajectoire x t y t l aide du programme graphique mis votre dis position 5 On se propose maintenant de r aliser une section de Poincar c est dire une coupe particuli re de l hypersurface ici 3 dimensions nergie constante On consid rera la coupe d finie par vs 0 avec v gt 0 choix du sens du temps le long de la trajectoire En pratique on crira dans un fichier les points x y de la trajectoire v rifiant vs lt et vy gt 0 Afin de disposer d un nombre de points suffisant on prendra les param tres suivants T 0001 e 0 0001 et quelques millions de pas d int gration 6 Visualiser les sections de Poincar pour diff rentes valeurs de
37. riv e y t 7 y t T f y t t OCT et on peut donc crire ylt T y t rflylt t O r Si y 0 est connue alors la trajectoire solution peut tre construite de mani re it rative et sans stockage d information On peut augmenter la pr cision en diminuant r et ou en faisant intervenir autant de d riv es d ordre sup rieurs que l on d sire Par exemple si on veut pousser la pr cision un ordre sup rieur T y t 7 y t ry Zy Or c est dire A ult 7 y t TFll t FUO t Or et ainsi de suite 2 Dynamique mol culaire On souhaite int grer les quations de Newton c est dire 2DN quations diff rentielles coupl es o D est la dimension de l espace et N le nombre de particules aap i 1 N at TS pesi II 3h avec l ensemble de conditions initiales 0 0 0 1 N Algorithme du 1 ordre F t 7 F t r t O r i 1 N Bilt T T t rdi t O r ou sa g n ralisation un ordre quelconque Algorithme de Verlet ou leapfrog saut de grenouille RT R Talt iral SIN dt T vi t dit amp t 7T V rifier qu il s agit d un algorithme du 2 ordre 17 3 Syst mes Hamiltoniens conservatifs et chaos a Nombre de degr s de libert gal un D 1 Un syst me Hamiltonien un degr de libert est int grable c est dire les trajectoires peuvent
38. s se NS NT ANA a a Ea aa 4 4 Temps d acc s finis l information en m moire 4 B Quelques l ments de Fortran 4 C Quelques l ments d Unix aoaaa a 8 I Systeme de fichiers eoii n is dah Do AR a N nt Ne dent 9 2A Process maa a a a Made ARE n r aaa a eia a aei a 10 D H f diteur detecte Ve e eurek a o a aaa M AT de aa dt a 10 E Quelques conseils de programmation ooo a a 11 F Aspects finis de l ordinateur uses nat Data A dan sen 13 12 Vit sse iile n r e gan g a SE MES MS SE a Ea 13 2 y Pre isiom Miena eie opr a a A E A E E A a a A 14 3 Temps d acces finis la m moire centrale et la m moire dure 15 II Dynamique mol culaire Chaos classique 16 A MOCH TANR S Le dE NN D DA A a E Meet ae db du 16 1 Equations diff rentielles 3 2 s sigma org unes nent nest 16 2 Dynaniquemol calair 22e La ed e tue se Matte e Mare 17 3 Syst mes Hamiltoniens conservatifs et chaos 18 B Apparition du chaos dans un syst me d oscillateurs coupl s 19 Document crit l occasion d un cours donn dans le cadre du DEA Champs Particules Mati res IV M thodes de Monte Carlo Transition de phase dans le mod le d Ising 2 D 21 A M thodes de Monte Carlo 80 25e 4 et As Note gra a 21 B Algorithme de Metropolis oaoa aa aaa de 28 RACE RAS 22 C mpl mentation pratique ooa a 24 D Simulation du modele d
39. s valeurs propres inf rieures ou gales en module 1 et que dans le cas o la densit II ne s annule pas le sous espace associ n est pas d g n r En cons quence PRET o c est un coefficient qui repr sente le recouvrement de la distribution f 0 VTI et du vecteur propre VII de la matrice M C Impl mentation pratique On se donne P une probabilit de transition appel e dans la suite probabilit de transition d essai On suppose que l on sait chantillonner cette probabilit c est dire tant donn e une configuration i tirer des configurations al atoires j distribu es selon la loi PF Comment passer d une configuration la suivante selon la probabilit de transition P d finie pr cedemment Les tapes sont les suivantes e Soit i la configuration courante On effectue un mouvement d essai trial move dans l espace des tats en tirant une nouvelle configuration j selon P ij Mj Pii 7 b A 7 A H e On calcule le rapport Ri 4 Le mouvement d essai est accept c est dire t ic la nouvelle configuration est j avec probabilit q Min 1 R Si le mouvement est refus alors la nouvelle configuration est l ancienne c est dire la configuration i Accepter avec probabilit q est r alis en comparant q un nombre al atoire uniforme tir entre 0 et 1 Si q gt le mouvement est accept sinon le mouvement est refus Attention Il est
40. steps xnew xold ranf 0 5 xdelta ratio proba xnew proba xold if ratio ge ranf then naccep naccep i else xnew xold endif valint valint f xnew g xold xnew enddo valint valint float nsteps print Integrale valint tauxaccep float naccep float nsteps print Taux acceptation tauxaccep print Simulation terminee end function proba x proba exp x 2 end function f x g f exp g x 4 end Ici ranf d signe une fonction intrins que d livrant un nombre al atoire compris entre 0 et 1 avec une distribution uniforme En g n ral les machines disposent toujours d une fonction de ce type dans leur biblioth que math matique de base videmment le nom de la fonction peut tre diff rent d une machine l autre D Simulation du mod le d Ising 2 dimensions On se propose d crire un programme pour calculer l nergie moyenne et la chaleur sp cifique du mod le d Ising deux dimensions Le modele est d crit par la fonction nergie suivante E _J Y S S 6 lt i 3 gt o les variables d Ising prennent les valeurs 1 ou 1 lt 4 3 gt signifiant que la sommation est prise sur les 4 plus proches voisins J est une constante positive int grale d change On consid rera un r seau carr compos de NV sites avec les notations suivantes 26 N 1 N 1 N N N O Ge De Or es dl TE N 1 N Nas G S OG p NH N 1 2N N4
41. sultats rewindi Le fichier est associ Open Permet d ouvrir un fichier dont le nom est crit entre l unit logique num ro 1 Une fois ouvert on pourra crire ou lire des donn es sur ce fichier par des ordres du type read 1 ou write 1 On peut ouvrir autant de fichiers qu on veut avec des num ros diff rents Rewind 1 est un ordre permettant de se positionner en d but de fichier dont l unit logique est 1 Chaque fois qu un l ment est lu ou crit dans le fichier un curseur de positionnement est incr ment termi 10 term2 20 print Entrez la valeur de n read n print ordre d criture le plus simple sur l cran signifie l unit 6 ou cran par d faut Mettre n importe quel texte entre read n ordre de lecture le plus simple sur l cran signifie l unit 5 ou cran par d faut do i i n tabi i sqrt float i enddo Les structures de boucle jouent un r le tr s important Forme plus g n rale do i idebut ifin increment o l incr ment est un entier positif ou n gatif Par d faut l incr ment est 1 sqrt ou float fonctions intrins ques pr d finies et charg es lors de l tape de load chargement Il existe de nombreuses fonctions intrins ques pr d finies sqrt cos sin exp float alog etc et leurs versions en double ou m me quadruple extended pr cision print Entrez la valeur de logic
42. t de transition ou matrice stochastique P i Pis 20 N ii gt Pi 1 ind pendant de i j e D f 3 Probabilit de transition ergodique Vio Vi il existe une probabilit non nulle qu apr s un nombre fini d application de la probabilit de transition partant de l tat o on arrive l tat 2 Enon ons maintenant l algorithme de Metropolis Soit P une probabilit de transition ergodique alors Pi d finie de la mani re suiv ante P P Min 1 Ri ji i j Frs FF F 2 1 Min 1 Rx j i I P avec Ri ee t43 est une probabilit de transition ergodique admettant Il comme probabilit stationnaire D monstration 1 Probabilit de transition F gt 0 vident gt 2 Ps D Pr Po gt Pr Pr zr 22 2 Probabilit stationnaire Il faut montrer que ILF Il Pour cela nous allons montrer tout d abord que le couple P 1 v rifie la condition dite de bilan d taill ILP IP V 5 D monstration i j galit vidente e iZj le rapport des deux membres de l galit pr c dente vaut IEP R Min 1 R IPs Min 1 Ry en notant que R 1 R et en distinguant les deux cas correspondant R gt let Ry lt 1 on v rifie facilement que ce rapport vaut un Finalement en utilisant la propri t de bilan d taill on a SILP P Il i i ce qui montre que Il est probabilit stationnaire La propri
43. t uil H lu gt et La matrice repr sentant H dans cette base est donc tridiagonale B2 az Bz 0 0 63 as Ha 0 B 0 Qi 0 Biyi 0 o ba 0 On an et il existe des m thodes performantes pour diagonaliser les matrices tridiagonales algo rithme QL par exemple voir Numerical Recipes Faisons un certain nombres de remarques concernant cet algorithme A priori si on s int resse essentiellement la d termination de l tat fondamental le vecteur initial Y gt sera choisi le plus proche possible de l tat fondamental recherch A la premi re it ration n 1 nergie obtenue sera l nergie variationnelle associ e au vecteur d essai c est dire g Lue v zut exacte A la seconde it ration n 2 le vecteur u gt est un facteur pr s le vecteur H Y gt E Y gt Si Y gt est le vecteur exact u1 gt est nul 8 0 et le processus est arr t on a construit un espace qui contient la solution exacte l espace vectoriel est stable ou ferm sous l action de H Dans le cas contraire u gt repr sente la direction privil gi e pour calculer le r sidu entre l tat variationnel et l tat exact En continuant ainsi on se rapproche tr s vite de la solution exacte En it rant 5 va essentiellement devenir de plus en plus petit Lorsque 5 s annule cela signifie que H appliqu au dernier vecteur Lancz s Cette nergie peut donc tre d ja une tr s bonne approximation de l nergi
44. te not i est le premier site du cluster en cours de construction e Les sites voisins not s j de ce site i sont incorpor s au cluster avec probabilit pl Si Sj 1 exp Min 0 8J on dit aussi que la liaison i j est activ e avec probabilit p Plus simplement si le spin du site j voisin de i a une valeur diff rente de celui de i le site j n est pas ajout au cluster Dans le cas contraire le site j est ajout avec probabilit 1 exp 8J e Cette proc dure d agr gation de sites au cluster est it r e Les liaisons avec les voisins ext rieurs au cluster de chacun des sites ajout s sont activ es ou non avec la probabilt d finie pr c demment La proc dure est it r e jusqu ce que le cluster se termine ceci est assur puisqu videmment la taille du cluster ne peut d passer le nombre de sites e Finalement tous les spins du cluster de N sites ainsi form sont retourn s on dit qu on fait un flip c est dire S S pour i C Le premier point important remarquer est qu un tel mouvement l mentaire dans l espace de configuration correspond un processus ergodique En effet notons S S SN et S 51 SN deux configurations quelconques de l espace de configu ration Il est clair que lors d un tape l mentaire il existe une probabilit non nulle de construire un cluster de taille unit compos d un site quelconque du syst me On a donc une pro
45. tion d un fichier ex cutable issu du compilateur Fortran xlf fichier f xlf est le nom du compilateur Fortran sur la machine que vous utiliserez IBM Risc Ce nom peut videmment varier d une machine l autre Noter que tout fichier soumis au compilateur Fortran doit avoir un nom se terminant par f xlf cr e un fichier qu il appelle a out par d faut a out est une commande ex cutable a out lance le programme La commande ps permet d avoir quelques formations sur les proces sus en cours d ex cution La commande kill permet de tuer un processus La commande kill doit tre suivie du num ro du processus tuer donn par la commande ps Noter aussi l existence de la commande time qui permet de savoir le temps total d ex cution de la commande concern e Par exemple time a out vous permet de conna tre le temps d ex cution de votre programme D L diteur de texte vi Pour crire un texte dans un fichier il faut un diteur de texte Il existe de nombreux diteurs de texte plus ou moins sophistiqu s Sous Unix il existe un diteur standard que vous retrouverez sous tous les syst mes partout dans le monde Il s agit de l diteur de texte vi qu il faut donc mon avis absolument conna tre Avec ce document est joint un document complet sur vi Je ne mentionnerai donc que le point important suivant Il existe essentiellement deux modes sous vi 1 Le mode norm
46. un grand nombre de degr s de libert en interaction Ce cours n est ni un cours d informatique ni un cours d analyse num rique Il existe de nombreux ouvrages sur ces aspects l et on peut s y reporter seul avec profit L objet de ce cours est plut t d illustrer comment des simulations simples mettre en oeuvre peu vent conduire des r sultats difficiles obtenir autrement On insistera sur les contraintes fondamentales qui conditionnent le choix de la m thode utiliser pour un probl me sp ci fique Un point important est de savoir valuer la d pendance de l effort num rique d une m thode donn e en fonction du parametre critique de la simulation Ici le param tre cri tique sera d fini comme le param tre conditionnant la quantit de m moire et le temps d ex cution n cessaire pour effectuer la simulation une valeur infinie du param tre corre spondant en g n ral une simulation exacte Concr tement ce param tre peut tre un nombre de particules un nombre de pas l mentaires de dynamique mol culaire etc Il est videmment important d utiliser la m thode optimale pour un probl me donn Choisir la mauvaise m thode conduit le plus souvent ne pas pouvoir en pratique effectuer la simula tion quantit de m moire requise et ou temps d ex cution d passant les limites physiques de l ordinateur En revanche une fois la bonne m thode choisie il est souvent possi ble d effectuer des simu
47. un livre quelconque sur le sujet suffiront rendre expert n importe quelle personne motiv e Cet effort pr liminaire est videmment indispensable si vous voulez profiter dans les meilleures conditions de ce cours Je pr sente maintenant un programme fictif surtout ne chercher aucune signification la t che accomplie utilisant quelques unes des principales conventions et structures du Fortran J ins re dans le texte les commentaires correspondants program ExEmPle Le texte principal d un programme Fortran except pour les signes de commentaires les tiquettes et les symboles pour indiquer les lignes suites commence la 7i me colonne et finit la 72i eme colonne M fiez vous de ne pas d passer cette 72i me colonne ExEmPle est le nom du programme principal Notez que le Fortran ne distingue pas entre les caract res crits en minuscule ou en majuscule Les mots ExEmPle Exemple exemple etc sont donc quivalents en Fortran zone type des variables Ceci est une ligne de commentaires qui n est pas prise en compte C est le premier caract re de la ligne qui d finit une ligne commentaire On peut utiliser indiff remment les caract res ou c integer termi term2 real icase rapid double precision pi deuxpi logical logic D claration des types de quelques variables En Fortran par d faut tous les mots dont le nom commence par une lettre a h sont consid r s comme r els les autres dans
48. ute quantit qui est commune plus de deux sous programmes Les autres quantit s doivent passer par argument entre le programme appelant et le programme appel F Aspects finis de l ordinateur Il est tr s important de bien comprendre les aspects finis de la machine Nous allons les illustrer sur des exemples simples 1 Vitesse finie On se propose d estimer l ordre de grandeur de la vitesse du processeur Pour cela on va faire tourner un petit programme calculant le nombre m l aide de la s rie infinie suivante 13 ce n est pas la s rie la plus rapidement convergente pour calculer xl r lim D 6 n o k2 Ecrivez un programme qui permet de calculer 7 calculer le nombre d op rations l mentaires flottantes effectu es dans votre programme pour un n maximum donn Faites tourner le programme et mesurez les temps d ex cution l aide de la commande time Une solution possible program calcpi implicit double precision a h o z print n read n nsort n 10 pi 0 d0 do i 1 n pi pi 1 d0 dfloat i 2 if mod i nsort eq 0 print i dsqrt pi 6 d0 enddo end Effectuer les calculs avec ou sans l ordre de sortie Ordres de grandeur 2 Pr cision finie Les calculs avec des entiers sont exacts sur ordinateur Ceci n est pas le cas avec des r els On se propose de calculer la solution g n rale de l quation aux diff rences finies suivante Tn 37tn 1 6 Z
49. veau Si vous devez prendre des semaines pour mettre en place votre algorithme vous abandonnerez tr s vite Il faut programmer juste Programmer juste c est mon avis d velopper une strat gie sophistiqu e de VERIFICA TION de votre programme Pour cela les points suivants me paraissent importants i Avant d crire quoi que ce soit murissez l algorithme dans votre t te aussi longtemps que la structure g n rale n est pas claire Ne vous jettez pas sur la machine ii D composer votre programme en autant d l ments l mentaires que possible Tester chacune des sous parties s par ment Tester les sous programmes in situ c est dire dans votre programme principal Il arrive souvent que le sous programme marche quand il est trait s parement mais donne des r sultats faux dans le programme principal par exemple variables globales qui interf rent ili Ecrivez un seul programme viter les multiples versions monprog f monprogl f monprog2 f monprog2a f etc Eviter de traiter les cas particuliers s par ment passer au cas le g n ral tout de suite Ne multiplier pas les commons les include les fichiers d entr e et de sortie iv V rifiez tous les cas limites connus exactement que votre programme g n ral peut produire Essayer de comprendre chaque d cimale une erreur relative de 1076 sur un r sultat en double pr cision est toujours instructive N oubliez pas que le chaos a t d couvert de cette
50. xiste de nombreux formats voir livres close 1 end Ordre de fermeture du fichier correspondant l unit logique 1 puis fin du programme prin cipal function expression k sum 0 do i 1 k sum sumtfloat i x xs enddo sum sum float k expression sum end subroutine progi n x y 2z integer termi term2 dimension x n y n 1 z n dimension xloc 10 common c1i termi term2 do i i n z i x i y i endo do 1i 1 10 z i z i xztab i enddo facteur float termi xterm2 do i 1 n z i Zz i facteur enddo end Sous programme progl noter comment les dimensions des tableaux pass s par argument sont d clar s Le tableau xloc 10 est un tableau local Ses dimensions ne peuvent pas tre d finies par une variable pass e par argument variable non locale C Quelques l ments d Unix Tout comme le Fortran le syst me d exploitation SE Unix est devenu une quasi norme dans la communaut scientifique Unix est un SE tr s simple utiliser En fait l id e fondamentale d Unix est que tout objet trait est consid r soit comme un fichier un fichier Fortran une commande ex cutable du syst me l cran une imprimante etc soit comme un processus commande en cours d x cution 1 Syst me de fichiers bin dev etc usr tmp unix boot vous paulot bob phle N tous vos fichiers Le Home Directory est le directory ensemble de fichiers et de directories dans lequel vous vous trouve
51. z quand vous vous loggez connectez la machine Pour savoir o vous tes pwd usr vous Le signe est ce qu on appelle le prompteur Le signe repr sentant le prompteur peut changer d une machine l autre et il peut tre red fini par l utilisateur Dans votre drectory vous pouvez faire ce que vous voulez mais videmment pas ailleurs Pour monter d un niveau dans l arbre des fichiers cd pwd usr Pour descendre d un niveau cd vous pwd usr vous vous tes revenus chez vous Vous pouvez aller partout dans le syst me Pour conna tre la liste des fichiers ou des directories l endroit o vous vous trouvez faire la commande ls En fait pour plus de d tails sur les fichiers ls 1 Pour conna tre le mode d emploi de toute commande Unix vous pouvez utiliser le manuel en ligne en tapant man suivi du nom de la commande ls l permet de savoir en autres choses si on a affaire un fichier ou un directory Copie d un fichier cp ficl fic2 Changement de nom mv ficl fic2 Destruction d un fichier rm ficl etc voir livre quelconque sur Unix 2 Processus Pour lancer un processus on lance une commande qui n est rien d autre qu un fichier ex cutable Le fichier ex cutable peut tre soit un fichier ex cutable syst me comme ls cp mv etc ou un fichier ex cutable que vous avez fabriqu L exemple le plus courant pour nous sera la cr a
Download Pdf Manuals
Related Search