蒙卡方法初步蒙卡方法介绍：粒子衰变盒中粒子蒙卡积分：中心极限定理，ImportanceSampling，VegasMC模拟：随机行走Markov链Metropolis高能物理事例产生子：权重，Unweight，直方图，Tree，bin，error拟合，Chi-2，Likelihood
蒙卡方法初步 蒙卡方法介绍： 粒子衰变 盒中粒子 蒙卡积分： 中心极限定理, Importance Sampling , Vegas MC模拟： 随机行走 Markov链 Metropolis 高能物理事例产生子: 权重, Unweight, 直方图, Tree, bin, error 拟合，Chi-2, Likelihood
在很多物理问题中，可以直接模拟物理过程，而不引言：粒子衰变用写下比如描述系统的微分方程。唯一的要求是需要知道几率分布函数。假设在初始时刻t=0，我们有N(O)个x核子，该种核子会辐射衰变，单位时间衰变几率为の，则有dN(t) = -oN(t)dt,1解为可定义平均寿命为=N(t) = N(O)e-ot0假设X核子辐射衰变到Y，Y也会衰变，则有dNx(t)= -0xNx(t),dtdNy(t)-0yNy (t) + OxNx(t)dt可以取大量样本，每个样本包含多个粒子，模拟这些粒子的演进：产生随机数与衰变几率比较
引言: 粒子衰变 假设在初始时刻t=0，我们有N(0)个X核子， 该种核子会辐射衰变，单位时间 衰变几率为ω，则有 解为 可定义平均寿命为 假设 X核子辐射衰变到Y，Y也会衰变，则有 可以取大量样本，每个样本包含多个粒子，模拟这些粒子的演进： 产生随机数 与 衰变几率 比较 在很多物理问题中，可以直接模拟物理过程，而不 用写下比如描述系统的微分方程。唯一的要求是需 要知道几率分布函数
引言：粒子衰变voidmc sampling(intinitial n_particles,intmax time,doubledecay_probability,int *ncumulative)(inttime,np,n_unstable,particle_limit90000n unstable=initial n particles;80000// accumulatethenumberof particlespertime70000step per trial60000ncumulative[o]=initial_n_particles;50000// loopovereachtimestep40000for (time=1;time<= max_time; time++)particlelimit=nunstable30000for(np=1; np<=particle_limit; np++)(20000if(double(rand()/RAND_MAX<=10000decay_probability)(n unstable=n unstable-1;101001人Time1//endofloopoverparticlesncumulative[time]=n_unstable;1//endof loopovertimesteps//endmcsamplingfunction
引言: 粒子衰变 void mc_sampling(int initial_n_particles, int max_time, double decay_probability, int *ncumulative) { int time, np, n_unstable, particle_limit; n_unstable = initial_n_particles; // accumulate the number of particles per time step per trial ncumulative[0] = initial_n_particles; // loop over each time step for (time=1; time <= max_time; time++){ particle_limit = n_unstable; for ( np = 1; np <= particle_limit; np++) { if( double(rand())/RAND_MAX <= decay_probability) { n_unstable=n_unstable-1; } } // end of loop over particles ncumulative[time] = n_unstable; } // end of loop over time steps } // end mc_sampling function
在很多物理问题中，可以直接模拟物理过程，而不引言：盒中粒子用写下比如描述系统的微分方程。唯一的要求是需要知道几率分布函数。t=o单位时间内，左边粒子移动到右边的几率为/N.Choosethenumberof particlesN.Make a loop over time,where the maximum time (or maximum number of steps)should belarger than thenumber of particlesN· For everytime step At there is a probabilityni/N for a move to the right.Comparethisprobabilitywith arandomnumberx..If x<n/N,decrease thenumber of particles inthe lefthalfbyone, i.e.,n=n-1.Else,moveaparticlefromtheright half tothe left,i.e.,n=n+l.Increasethetimebyoneunit(theexternal loop)
引言: 盒中粒子 在很多物理问题中，可以直接模拟物理过程，而不 用写下比如描述系统的微分方程。唯一的要求是需 要知道几率分布函数。 单位时间内，左边粒子移动到 右边的几率为 / l n N
引言：盒中粒子int main(int argc, char* argvll)(srand(unsigned(time(0)));// seedtherandomizerint initial_n_particles,max_time;inttime,random_n,nleft;Particle In Boxofstreamofile("PinBox.dat")10000cout<<"initialnumber of particles="<<endl;sim9500exactcin >>initial n_particles;9000nleft=initial_n_particles8500max time=2*initial n particles8000for(time=0;time<=max_time;time++)N7500randomn=(int)(initial nparticles7000*double(rand()/RAND_MAX);6500if (random_n <= nleft)6000nleft -= 1;550015000101001000110000100000else (Timenleft += 1;1ofile<<setiosflags(ios::showpoint/ios::uppercase)ofile << setw(15)<< time;ofile<<setw(15)<<nleft;return O;1//endmainfunction
引言: 盒中粒子 int main(int argc, char* argv[]) { srand(unsigned(time(0))); // seed the randomizer int initial_n_particles, max_time; int time, random_n, nleft; ofstream ofile("PinBox.dat"); cout << "Initial number of particles = " << endl ; cin >> initial_n_particles; nleft = initial_n_particles; max_time = 2*initial_n_particles; for( time=0; time <= max_time; time++){ random_n = (int)(initial_n_particles *double(rand())/RAND_MAX); if ( random_n <= nleft){ nleft -= 1; } else { nleft += 1; } ofile << setiosflags(ios::showpoint | ios::uppercase); ofile << setw(15) << time; ofile << setw(15) << nleft; } return 0; } // end main function