G. Pastore e M. Peressi Universita’ di Trieste a.a. 2010/11

1. La simulazione numerica in fisica:approcci stocastici(aggiornamento e formazione insegnanti;corso CIRD - CP e M IDIFO3: Lab_A I15) http://www.laureescientifiche.units.it/ • G. Pastore e M. Peressi • Universita’ di Trieste • a.a. 2010/11 14 marzo 2011

2. Approcci stocastici Preferibili per: - fenomeni intrinsecamente casuali - sistemi con gran numero di gradi di liberta’ in cui l’interesse non e’ focalizzato sul singolo valore ma su proprieta’ medie dell’insieme Utilizzano: - numeri casuali

3. Numeri casuali Una sequenza di numeri casuali e’ una sequenza di numeri che sembrano impredicibili ma che hanno ben definite proprieta’ statistiche

4. Generazione di numeri casuali metodi manuali: - lancio di dadi, di monete... - estrazione da urne, da tabelle... Ottime proprieta’ di casualita’, ma... lunghi !!! praticamente inutilizzabili!!! - dispositivi elettronici ad hoc... NO! altri problemi... generazione con il computer: Generazione mediante metodi aritmetici

5. Numeri pseudocasuali e proprieta’ statistiche Anche se si lavora con sequenze di numeri pseudocasuali (cosi’ sono quelli generati al computer, perche’ deterministicamente generati), va bene a patto che abbiano ben definite proprieta’ statistiche: 1) non correlazione tra i numeri 2) determinata distribuzione e relative proprieta’ (es. : uniforme) 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 <r>

6. Numeri pseudocasuali e algoritmi per la generazione un po’ di storia in brevissimo...

7. (pseudo)Random numbers generation

8. (pseudo)random numbers generation: example 1 - “Middle square” algorithm (Von Neumann, 1946) To generate a 10-digit integer sequence: - take a first one - square it - take the middle 10 digits of the result loop Limits of the algorithm: depending on the initial choice, you can be trapped into short loops:

9. (pseudo)random numbers generation: example I1 - “Linear congruential method” (Lehemer, 1948) In+1 = (a In + c) mod m Limits of the algorithm:

10. Numeri pseudocasuali e problemi... correlazione in iperpiani (routine IBM)

11. More subtle limits of some algorithms... xi, p(xi) (xi, xi+1)

13. Numeri pseudocasuali e algoritmi per la generazione il Linear Congruential Method implementato in modo semplice in: random_lc.f90 La funzione intrinseca del Fortran90: modulo ((a*old+c), m) restituisce il resto della divisione tra i due argomenti

14. Artigianale e’ bello, ma non troppo: correlazioni con LCM e M=256 Quanti numeri generati? Quante coppie?

15. Numeri pseudocasuali e procedure intrinseche Tipicamente ci sono procedure intrinseche un po’ piu’ sofisticate di questo nostro LCM che generano numeri casuali tra [0,1] con distribuzione uniforme e migliori proprieta’ statistiche Es. : in Fortran90: la subroutine random_number() in Pascal: la funzione Random

16. Numeri pseudocasuali e procedure intrinseche Es. in Fortran90: la subroutine random_number() Si puo’ dichiarare un ‘seme’ (seed) di partenza rantest_intrinsic.f90 e rantest_intrinsic_con_seed.f90

17. A cosa servono i numeri pseudocasuali? • per calcolare integrali definiti... 4 x Area settore circ. Area quadrato # punti nel settore circ. # punti nel quadrato PiMonteCarlo.java Pi.f90 in onore del “pi-day” ;-)

18. A cosa servono i numeri pseudocasuali? • flusso di particelle tra due meta‘ di una scatola: box.f90 • .. • per simulare processi casuali:

19. Come si usano i numeri pseudorandom per simulare eventi che avvengono con una certa probabilita’ (prob)? 1) genero un numero pseudorandom r tra 0 e 1 2) lo confronto con prob 3) se r e’ < o = prob , l’evento accade, altrimenti no; genero un altro numero random e itero... call random_number(r) if (r <= prob) then (accade...) else (non accade...) end if

20. Macroscopic systems towards equilibrium A simple example: non-interacting classical particles in a box (gas diffusion) Another version: particles black/white in both sides (interdiffusion of two gases): per unit time, one from each side is picked at random and put in the other side: Nleftwhite(t)+Nleftblack(t)=constant; Nleftwhite(t)=?

21. Approach to equilibrium - I macrostate: specified by the number of particles n on the left side; microstate: specified by the specific list of the n particles on the left side

22. Approach to equilibrium - II initial N(left)=1000 How to reduce fluctuations? - more particles - average over many simulation runs

23. Equilibrium and entropy The number of microstates for the “particle in a box” model with N=10. The macrostate is specified by the number of particles on the left side, n. The total number of microstates for N=10 is 210=1024 the most “random”! Equilibrium = Maximum number of possible microstates = Maximum entropy

24. A cosa servono i numeri pseudocasuali? • flusso di particelle tra due meta‘ di una scatola... • decadimento radioattivodecay.f90 • ... • per simulare processi casuali:

25. Random processes: radioactive decay Atoms present at time Probability for each atom to decay in Atoms which decay between and we use the probability of decay of each atom to simulate the behavior of the number of atoms left; we should be able to obtain (on average):

26. Radioactive decay:numerical simulation Note: unbounded loop DO ! loop on time DO i = 1, nleft ! loop on all the nuclei left call random_number(r) IF (r <= lambda) THEN ! BASIC ALGORITHM nleft = nleft -1 ! update the nuclei left (*) ENDIF END DO WRITE (unit=7,fmt=*) t , nleft if (nleft == 0) exit t = t + 1 END DO Note: “exit” =/= “cycle” (*) Notice that the upper bound of the inner loop (nleft) is changed within theexecution of the loop; but with most compilers, in the execution the loop goes on up to thenleftset at the beginning of the loop; this ensures that the implementation of thealgorithm is correct. The programprova.f90in the same directory is a test for the behavior of the loop. Look also at decay-new.f90. If nleft would be changed (decreased) during the execution, the effect would be an overestimate of teh decay rate. CHECK with your compiler!

27. Details on Fortran: unbounded loops [name:] DO exit [name] or [name:] DO END DO [name] possible forms of "do while": DO if (condition)exit END DO or: DO WHILE (.not. condition) ... END DO NOTE: first is better (if () exit can be placed everywhere in the loop, DO WHILE must execute the loop up to the end - Other note: Difference between EXIT and CYCLE

28. Radioactive decay:results of numerical simulation plot of the results of decay simulation (N vs t) with N=1000 N(t) ~ N0 exp(a t) semilog plot (log(N) vs t) => log(N(t)) = log N0 + a t => slope is a N(t) t log(N(t)) t

29. Radioactive decay:results of numerical simulation Semilog plot of the results of decay simulation for the same decay rate and different initial number of atoms: almost a straight line, but with important deviations (stochastic) for small N numerical simulations: OK on average and for large numbers

30. A cosa servono i numeri pseudocasuali? • flusso di particelle tra due meta‘ di una scatola... • decadimento radioattivodecay.f90 • random walk segue in dettaglio • per simulare processi casuali:

31. Random walks • come simulare? • quali proprieta’ interessano?

32. Random walks unidimensionali Un camminatore puo’ andare a destra o a sinistra: : numero di passi : lunghezza del singolo passo (direzione casuale) ( relativo spostamento al passo -simo) : spostamento dal punto di partenza dopo passi ( ) , : probabilita’ di spostamento a destra o a sinistra

33. Random walks unidimensionali Cosa possiamo calcolare PER OGNI CAMMINATORE dopo N passi ? : distanza DI UN CAMMINATORE dalla partenza : la sua distanza QUADRATICA dalla partenza

34. Random walks unidimensionali Cosa possiamo calcolare PER PIU’ CAMMINATORI dopo N passi ? : spostamento MEDIO dopo N passi : spostamento QUADRATICO MEDIO dopo N passi ATTENZIONE!!! proporzionale al numero di passi

35. Random walks unidimensionali Piu’ generale:

36. Random walks Un’altra grandezza calcolabile SU PIU’ CAMMINATORI: : probabilita’ che IN MEDIA sia lo spostamento finale dal punto di partenza dopo N passi (ma lasciamo da parte per ora ...) In piu’ dimensioni??? per il singolo camminatore per piu’ camminatori si definiscono anche in due o piu’ dimensioni:

37. RW 1D: simulazione The basic algorithm: ix = position of the walker x_N, x2_N = cumulative quantities rnd(N) = sequence of N random numbers (1 run= 1 particle)

38. RW 1D: simulazione rw1d.f90

39. Proprieta’ statistiche

40. Un’altra proprieta’ statistica dei random walks (qui 1D) Random walks unidimensionali E’ facile calcolare analiticamente la probabilita’ per un random walk di N passi di finire in un punto x: si tratta di un triangolo di Pascal con degli zeri !

41. Un’altra proprieta’ statistica dei random walks (qui 1D) Random walks unidimensionali Per lunghi cammini (N grande) la probabilita’ che sia lo spostamento finale dal punto di partenza dopo N passi ha una distribuzione simmetrica rispetto alla media e con “larghezza” (“deviazione standard”,SD)

42. La tavola di Galton Simulazione della distribuzione di probabilita’ della posizione finale di un random walk unidimensionale (con passo variabile) Biglie che partono da una stessa posizione centrale e cadono deviate da chiodi disposti secondo uno schema triangolare. Ogni volta che una sferetta colpisce un chiodo ha il 50% di probabilità di cadere a sinistra e il 50% a destra. Le palline si accumulano nelle scanalature collocate alla base della struttura, formando delle pile con una certa distribuzione. In questa simulazione le pile nelle scanalature sono la somma di N=8 variabili casuali, ma se si aggiungono più file di chiodi (cioè se si aumenta N, il numero di variabili casuali), l'andamento della distribuzione approssima quello della distribuzione normale. poche tante

43. al solito, attenzione alla statistica;gia’ abbiamo visto per quanto riguarda Δ2N;vediamo anche per PN:

44. RW 1D: simulation

45. RW 1D: simulation

46. RW 1D: simulation

47. RW 1D: simulation

48. Random Walks 2D

49. Random Walks 2D Per piu’ camminatori:

50. Random Walks 2D Per piu’ camminatori: