计算物理ComputationalPhysics计算物理 (9).pdf

《计算物理ComputationalPhysics计算物理 (9).pdf》由会员分享，可在线阅读，更多相关《计算物理ComputationalPhysics计算物理 (9).pdf（43页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、Computational physiCsDice骰子Monte Carlo simulations Sampling and integration Darts method to calculate p Random-number generators The Metropolis algorithm Ising modelWhat is Monte Carlo?Monte-Carlo,MonacoMonte Carlo Casino Monte Carlo method:computational statisticsBirthday:World War II,1940sFather:S

2、.UlamBirth Place:neutron reaction in atom bomb,Manhattan ProjectStanislaw Ulam(1909-1984)1.Teller-Ulam configuration2.Yu MinconfigurationYu Min1926-2019PrincipleWhy is it called Monte Carlo?Physical Problems Simulated by Random Processes Numerical solutions Sampling and integration If we want to fin

3、d the numerical value of the integralDivide the region 0,1 evenly into M slices with x0=0 and xM=1,and then the integral can be approximated aswhich is equivalent to sampling from a set of points x1,x2,.,xM in the region0,1 with an equal weight,in this case,1,at each point.We can also select xn with

4、 n=1,2,.,M from a uniform random number generator in the region 0,1 to accomplish the same goal.If M is very large,we would expect xn to be a set of numbers uniformly distributed in the region 0,1 with fluctuations proportional to 1/sqrt(M).Then the integral can be approximated by the averagewhere x

5、n is a set of M points generated from a uniform random number generator in the region 0,1.Example In order to demonstrate the algorithm clearly,let us take a very simple integrand f(x)=x2 from 0 to 1.The exact result of the integral is 1/3.The following program implements such a sampling.8.1.MC.cpp

6、8.rand.hA history of p(pie):TimeWhoValue2000-1850 BCEgyptian阿美斯纸草书256/81=3.160493.1900 BCBabylonian25/8=3.125900 BCIndian(Shatapatha Brahmana)339/108=3.13888600 BCHebrew希伯来圣经3250 BC3.14163491100 BC?Chinese周髀算经3径一而周三TimeWhoValue20 BCVitruvius3.12550-23 BC刘歆3.1547130 AD?张衡3.146551150 AD3.141666250 AD王

7、蕃3.155555263 AD刘徽3.141024 p 3.142104400 AD何承天111035/35329=3.142885.480 AD祖冲之3.1415926 p 3.1415927A record for one thousand years!Zu Chongzhi Born in Nanjing!(南朝刘宋）大明历 p=335/113(祖率)The record holds for about 1000 years.Broken by Jamshid Masud Al Kashi(Arab)in 1424.The history of science,is a flag of

8、human civilization!World-leading scientists were very rare in the Chinese long history!Chinese name in Moon:石申、张衡、祖冲之、郭守敬、万户Chinese name in Mercury:李清照The first rocket man!We have largest population!We have one of longest history and unbroken culture!We have the 2nd GPD now!We have the 1st-3rd golde

9、n medals in Olympic games!We need more Scientists&philosophers!We need some people looking at sky!Random-number generatorsNot a number but a sequence with random numbersRandom number generators:Real random number(RRN)generators:Quantum processes like nuclear decay Pseudo-random number(PRN)generators

10、:Lottery,computational PRN:depends on the initial one,algorithm,and word lengthGood PRNBad PRNThe rand function returns a pseudorandom integer in the range 0 to RAND_MAX.Use the srand function to seed the pseudorandom-number generator before calling rand.Or you will always get the same number sequen

11、ce!PRN function in C/C+void srand(unsigned int seed);int rand(void)1.How to generate a real PRN in 0-1?2.How to generate an integer PRN in 0-6?3.How to generate a binary PRN-1/1?rand()/32767 1.0*rand()/32767 1.0*rand()/RAND_MAX8.3.Rand.cppImportance sampling Is there any way to increase the accuracy

12、 for some specific types of integrands?If F(R)is smooth and close to constant,the accuracy from the Monte Carlo quadrature would be much higher.In many cases,the function F(R)is not a smooth function.The idea of importance sampling introduced by Metropolis et al.(1953)is to sample the points from a

13、nonuniform distribution.Nicholas Metropolis(1915-1999)The algorithm by Metropolis(1953)has been cited as among the top 10 algorithms having the greatest influence on the development and practice of science and engineering in the 20th century.Metropolis algorithm)()()(11iMiiiRWRWRFMSwhere M is the to

14、tal number of points sampled according to the distribution function W(R).Rewrite the integral as:If a distribution function W(R)can mimic the drastic changes in F(R),we should expect a much faster convergence withwhere W(R)is positive definite and satisfies the normalization condition.G(R)=F(R)/W(R)

15、.We can imagine a statistical process that leads to an equilibrium distribution W(R)and the integral S is merely a statistical average of G(R).This can be compared with the canonical ensemble average.where A(R)denotes the physical quantities to be averaged.The probability or distribution function W(

16、R)is given by:with U(R)being the potential energy of the system for a given configuration R.Here kB is the Boltzmann constant and T is the temperature of the system.The selection of the sampling points is viewed as a Markov process.In equilibrium,the values of the distribution function at different

17、points of the phase space are related bywhere T(R R)is the transition rate from the state characterized by R to the state characterized by R.This is usually referred to as detailed balance in statistical mechanics.Metropolis algorithmNow the points are no longer sampled randomly but rather by follow

18、ing the Markov chain.The transition from one point R to another point R is accepted if the ratio of the transition rates satisfieswhere wi is a uniform random number in the region 0,1.Outline of the stepsWe first randomly select a configuration R0 inside the specified domain D.Then W(R0)is evaluated

19、.A new configuration R1 is tried with R1=R0+DR,where DR is between h,h.The actual value of h is determined from the desired accepting rate.In practice,h is commonly chosen so that the accepting rate of moves is around 50%.We can evaluate the probabilityComparing p with a uniform random number wi in

20、the region 0,1.If p wi,the new configuration is accepted;otherwise,the old configuration is assumed to be the new configuration.This procedure is repeated and the physical quantity A(Rk)for k=n1,n1+n0,.,n1+(M-1)n0 is evaluated.The numerical result of the integral is then given byNote that the first

21、n1 steps are used to remove the influence of the initial selection.The data are taken n0 steps apart to avoid high correlation between the data points.In most cases,the distribution function W(R)varies by several orders of magnitude,whereas A(R)stays smooth or nearly constant.This sampling by import

22、ance is much more efficient than the direct,random sampling presented in the preceding section.with f(x)=x2.We can choose the distribution function aswhich is positive definite.The normalization constant Z is given by:Now we would like to compare this procedure numerically with the direct,random sam

23、pling.We still consider the integralThen the corresponding function g(x)=f(x)/W(x)is given by:Now we are ready to put all of these into a program that is a realization of the Metropolis algorithm for the integral specified.8.4.Metropolis.cppHomework 1.Generate random integer numbers:0,1,2 with the p

24、ossibility 25%,50%,25%,respectively.2.Generate random numbers from 0 to 2 with the given distribution:exp-(x-1)2Applications in statistical physicsIn this section,we discuss some applications of the Metropolis algorithm in statistical physics.We will use the Ising model as the illustrative example.T

25、he Ising model is a discrete lattice system.Ising modelWilhelm Lenz 1888-1957Ernst Ising,1900199819201924 Ernst Ising was born in Cologne in 1900.After school,he studied physics and mathematics at the University of Gttingen and University of Hamburg.In 1922,he began researching ferromagnetism under

26、the guidance of Wilhelm Lenz.He earned a Ph.D in physics from the University of Hamburg in 1924 when he published his doctoral thesis(an excerpt or a summary of his doctoral thesis was published as an article in a scientific journal in 1925 and this has led many to believe that he published his full

27、 thesis in 1925).His doctoral thesis studied a problem suggested by his teacher,Wilhelm Lenz.He investigated the special case of a linear chain of magnetic moments,which are only able to take two positions,up and down,and which are coupled by interactions between nearest neighbors.Mainly through fol

28、lowing studies by Rudolf Peierls,Hendrik Kramers,Gregory Wannier and Lars Onsager the model proved to be successful explaining phase transitions between ferromagnetic and paramagnetic states.After earning his doctorate,Ernst Ising worked for a short time in business before becoming a teacher,in Sale

29、m,Strausberg and Crossen,among other places.In 1930,he married the economist Dr.Johanna Ehmer(February 2,1902-February 2,2012).As a young German-Jewish scientist,Ising was barred from teaching and researching when Hitler came to power in 1933.In 1934,he found a position,first as a teacher and then a

30、s headmaster,at a Jewish school in Caputh near Potsdam for Jewish students who had been thrown out of public schools.Ernst and his wife Dr.Johanna Ising,lived in Caputh near the famous summer residence of the Einstein family.In 1938,the school in Caputh was destroyed by the Nazis,and in 1939 the Isi

31、ngs fled to Luxembourg,where Ising earned money as a shepherd and railroad worker.After the German Wehrmacht occupied Luxembourg,Ernst Ising was forced to work for the army.In 1947,the Ising family emigrated to the United States.Though he became Professor of Physics at Bradley University in Peoria,I

32、llinois,he never published again.Ising died at his home in Peoria in 1998,just one day after his 98th birthday.Spin:si=1Ising,ijijii jiHJ S ShS Ising model HamiltonianAny physical observationsBsBTTksHsATksHsA/)(exp)(/)(exp)(s=s1,s2,si,Total states:2NGrains on chessboardcanonical ensemble partition f

33、unction The Ising model was used historically to study magnetic phase transitions.The magnetization is defined aswhich is a function of the temperature and external magnetic field.For the h=0 case,there is a critical temperature TC that separates different phases of the system.For example,the system

34、 is ferromagnetic if T Tc,and unstable if T=Tc.The complete plot of T,m,and h forms the so-called phase diagram.Another interesting application of the Ising model is that it is also a generic model for binary lattices:that is,two types of particles can occupy the lattice sites with two on-site energ

35、ies which differ by 2h.So the results obtained from the study of the Ising model apply also to other systems in its class,such as binary lattices.Other quantities of interest include,but are not limited to,the total energy E=and the specific heatIn order to simulate the Ising model,for example,in th

36、e calculation of m:where s=S/N,with S=si being the total spin of a specific configuration labeled by and H being the corresponding Hamiltonian(energy).The above summation is over all the possible configurations.Here Z is the partition function given by The average of a physical quantity,such as the

37、magnetization,can be obtained fromwith =1,2,.,M indicating the configurations sampled according to the distribution function Start from a given spin configuration.Choose a spin.Change it.Calculate the energy change DE=new-old.Compare the possibility p=exp-DE/kBT with a PRN.If pPRN,reject the change,

38、else,accept the new configuration.Repeat the above processes for many steps.Perform measurements and calculate average.Outline of the stepsLars Onsager(November 27,1903 October 5,1976)was a Norwegian-born American physical chemist and theoretical physicist,winner of the 1968 Nobel Prize in Chemistry

39、.He held the Gibbs Professorship of Theoretical Chemistry at Yale University.During the 1940s,Onsager studied the statistical-mechanical theory of phase transitions in solids.He obtained the exact solution for the two dimensional Ising model in zero field in 1944.for 1D chain:TC=0 (Ising proved)for

40、2D square lattice:KBTC/J=2/ln(1+sqrt(2)=2.26918531421 (Onsager proved)for 3D lattice:?Code example 8.5.Ising.cpp 8.6.Flip.javaHomework Write a MC code for a two-dimensional Heisenberg spin lattice with perioidc boundary conditions.izizijjiijijjiShSSDSSJHThe vector Dij points from i to j.SkyrmionA magnetic vortex with a topological monopole