Monte Carlo Methods

Monte Carlo Methods Mark A. Novotny Dept. of Physics and Astronomy Mississippi State U. man40@ra.msstate.edu Random topics about random things for random people Supported in part by NSF DMR0120310 and ITR/AP(MPS)0113049

Where is Monte Carlo? • Europe • Principality of Monaco • Monte Carlo is 1 of 5 regions of Monaco • Monte Carlo founded in 1866 by Prince Charles III • Renowned casino, luxurious hotels, beaches, …

What is a Monte Carlo simulation? In a Monte Carlo simulation we attempt to follow the `time dependence’ of a model for which change, or growth, does not proceed in some rigorously predefined fashion (e.g. according to Newton’s equations of motion) but rather in a stochastic manner which depends on a sequence of random numbers which is generated during the simulation. Landau and Binder, A Guide to Monte Carlo Simulations in Statistical Physics (Cambridge U. Press, Cambridge U.K., 2000), p. 1.

Los Alamos National Laboratory 1953 Physical Review Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller Conferences/MonteCarloMethods/ Birth of the Monte Carlo method “The only good Monte Carlo is a dead Monte Carlo.” Trotter and Turkey, 1954

Random numbers Integrals via Monte Carlo Importance Sampling Monte Carlo Dynamic Monte Carlo Decay of nuclei (birth-death processes) Random motion in 1 dimension Long-time simulations (Monte Carlo with Absorbing Markov Chains, MCAMC) Parallel dynamic Monte Carlo simulations OUTLINE

Uniformly distributed numbers in [0,1] How good is `good enough’? `religious question’ Linear congruential, R250, Marsaglia, 4-tap, system supplied (not-really, but almost) Random Numbers Anyone who considers arithmetical methods of producing random digits is, of course, in a state of sin. John von Neumann (1951)

Xn = XOR( Xn-p, Xn-q ) XOR(,) is exclusive OR operator With p2+q2+1=prime (p and q are Mersine primes) R98: p=98, q=27 R250: p=250, q=103 R1279: p=1279, q=216 or 418 R9689: p=9689, q=84, 471, 1836, 2444, or 4187 (not-really, but almost) Random Numbers from R250

Set Nin=0 Do N times Calculate 3 random numbers, r1, r2, r3 Let x=r1 Let y=r2 Use r3 to choose quadrant (change signs of x and y) If x2+y2.le.1 set Nin=Nin+1 Estimate for p=pi = 4Nin / N Calculate p=pi using Monte Carlo

Calculate p using Monte Carlo N=104 N=105 N=103 Live or die by the Law of Large Numbers

Importance Sampling Monte Carlo • Old way: choose points with equal probability from rectangle with a<xi<b and 0<y<y0; and then use estimate yestimate=Sif(xi) (b-a) y0 • Rather choose points with importance of the value of the function at that point to the integral, p(x) • Estimate of integral yestimate=Sip-1(xi) f(xi) y0 0 a b

Dynamic Monte Carlo • For Static Monte Carlo the order of generation of points does not matter (like finding the integrals) • For Dynamic Monte Carlo the order does matter • This gives a Markov chain method, governed by the Master Equation

Decay (death) of nuclei (organisms) rate l=0.001 /msec N0=100 N(t) = N0exp(-l t) t1/2 = ln(2)/l = 693.15 Birth and Death processes

Start with N=N0individuals Start time t=0 Assume that decay (death) constant l<1 For each nucleus which has not yet decayed Generate a random number 0<r<1 Nucleus decays if r<l, and set N=N-1, else the nucleus remains When done sampling all remaining nuclei, set t=t+1 Repeat until all nuclei have decayed Birth and Death processes

Start with N=N0individuals Start time t=0 Assume that decay (death) constant l<1 For each nucleus which has not yet decayed Generate a random number 0<r<1 Nucleus decays if r<l, and set N=N-1, else the nucleus remains When done sampling all remaining nuclei, set t=t+1 Repeat until all nuclei have decayed Birth and Death processesEvent Driven --- n-fold way May be VERY SLOW if decay rate l is small!!! Does the nucleus decay? No no no no no no no no no no no no no no no ...

Start with N=N0individuals Start time t=0 Assume that decay (death) constant l<1 Until all the nuclei have decayed do: Generate two random number 0<r1 ,r2<1 Calculate time to leave current state: Dt=-ln(r1)/(N l) Set t = t + Dt Set N = N - 1 Use r2to pick which of remaining nuclei decayed Birth and Death processesEvent Driven -- or -- n-fold way One decay every algorithmic step, no matter how small l is

Birth and death of organisms lgrow=0.001 /msec ldie =0.001 /msec N0=100 Birth and Death processes

Random Walker coagulation Albena, Bulgaria model (2002)

Random Walker coagulation

Random Walker coagulation Event driven (n-fold way) One algorithmic step step from current configuration

Random Walker coagulation Monte Carlo with Absorbing Markov Chains (MCAMC) (s=2) In one algorithmic step step out of low valleys

Random Walker coagulation Average time until coagulation All algorithms statistically the same

Random Walker coagulation Different algorithms require different amounts of computer time

Random motion --- ion channel flow Random motion (of a random walker) through a channel

Random motion --- ion channel flow

Random motion --- ion channel flow Average lifetimes the same for all algorithms Computer times different for different algorithms

Random motion --- ion channel flow Square lattice for random walker, lifetime is time to get from one end to another

Monte Carlo for Ising models

Monte Carlo for Ising models Use spin classes for Ising model --- square lattice has 10 spin classes Every spin belongs to one of 10 classes Probability of flipping a spin in each class is the same If chosen, probability of flipping spin in class i is p(i) Number of spins in class i is ci

Monte Carlo for Ising models A GAME

Monte Carlo for Ising models MCAMC n-fold way

Monte Carlo for Ising models n-fold way with needed Bookkeeping

Monte Carlo for Ising models MCAMC

Monte Carlo for Ising models Square lattice Ising model

Monte Carlo for Ising models MCAMC results square Ising

Monte Carlo for Ising models Age of universe in femtoseconds? 10-15 seconds

Dynamic Monte CarloParallel Discrete Event Simulations • Dynamics for Materials • Dynamics for biological and ecological models • Dynamics of Magnets • Cell-phone switching • Spread of infectious diseases • Resource allocation following terrorist attacks • War-game scenarios

Best is trivial parallelization Each processor runs same program with different random number sequence Parallel Discrete Event Simulations computers & algorithms “Nature” ?

Two approaches to parallelization d=1 d=2 i (site index)i • Optimistic (or speculative) • PEs assume no causality violations • Rollbacks to previous states once causality violation is found (extensive state saving or reverse simulation) • Rollbacks can cascade (“avalanches”) • Conservative • PE “idles” if causality is not guaranteed • utilization, u: fraction of non-idling PEs

Monte Carlo Methods