La simulazione numerica in fisica:

approcci stocastici

(aggiornamento e formazione insegnanti;corso CIRD - CP e M IDIFO3: Lab_A I15)

G. Pastore e M. PeressiUniversita’ di Trieste

a.a. 2010/11

http://www.laureescientifiche.units.it/

14 marzo 2011

Approcci stocasticiPreferibili 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

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

Numeri casuali

Generazione di numeri casuali

- 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...

metodi manuali:

generazione con il computer:

Generazione mediante metodi aritmetici

QuickTime™ and aGIF decompressorare needed to see this picture.QuickTime™ and aGIF decompressorare needed to see this picture.

QuickTime™ and aGIF decompressorare needed to see this picture.

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:

0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95

Numeri pseudocasualie proprieta’ statistiche

1) non correlazione tra i numeri 2) determinata distribuzione e relative proprieta’ (es. : uniforme)

<r>

Numeri pseudocasualie algoritmi per la generazione

un po’ di storia in brevissimo...

(pseudo)Random numbers generation

(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

Limits of the algorithm: depending on the initial choice, you can be trapped into short loops:

loop

(pseudo)random numbers generation:

example I1 - “Linear congruential method”(Lehemer, 1948)

In+1 = (a In + c) mod m

Limits of the algorithm:

Numeri pseudocasualie problemi...

correlazione in iperpiani (routine IBM)

More subtle limits of some algorithms...

xi, p(xi)

(xi, xi+1)

Numeri pseudocasualie algoritmi per la generazione

il Linear Congruential Method implementatoin 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

Artigianale e’ bello, ma non troppo:correlazioni con LCM e M=256

Quanti numeri generati? Quante coppie?

Tipicamente ci sono procedure intrinseche un po’ piu’ sofisticate di

questo nostro LCMche generano numeri casuali tra [0,1]con distribuzione uniforme e migliori

proprieta’ statisticheEs. : in Fortran90:la subroutine random_number()

in Pascal: la funzione Random

Numeri pseudocasualie procedure intrinseche

Es.

in Fortran90:la subroutine random_number()Si puo’ dichiarare un ‘seme’ (seed) di partenza

Numeri pseudocasualie procedure intrinseche

rantest_intrinsic.f90e

rantest_intrinsic_con_seed.f90

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” ;-)

A cosa servono i numeri

pseudocasuali?

•flusso di particelle tra due meta‘ di una scatola: box.f90

•..

•per simulare processi casuali:

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 12) lo confronto con prob3) 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

Macroscopic systems towards

equilibrium

A box is divided into two parts communicating through a small hole. One particle randomly can pass through the hole per unit time, from the left to the right or viceversa.

Nleft(t): number of particles present at time t in the left sideGiven Nleft(0), what is Nleft(t) ?

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: Nleft

white(t)+Nleftblack(t)=constant; Nleft

white(t)=?

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

initial N(left)=1000

Approach to equilibrium - II

How to reduce fluctuations?- more particles- average over many simulation runs

Equilibrium and entropy

n # of microstates

log(# of micr.)

0 1 0.00

1 10 2.30

2 45 3.81

3 120 4.79

4 210 5.35

5 252 5.53

6 210 5.35

7 120 4.79

8 45 3.81

9 10 2.30

10 1 0.00

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

A cosa servono i numeri

pseudocasuali?

•flusso di particelle tra due meta‘ di una scatola...

•decadimento radioattivo decay.f90

• ...

•per simulare processi casuali:

Random processes: radioactive decay

Atoms present at time Probability for each atom to decay inAtoms which decay between and

we use the probability of decay of each atomto simulate the behavior of the number of atoms left;we should be able to obtain (on average):

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

Radioactive decay:numerical simulation

Note: “exit” =/= “cycle”

Note: unbounded loop

(*) Notice that the upper bound of the inner loop (nleft) is changed within the execution of the

loop; but with most compilers, in the execution the loop goes on up to the nleft set at the beginning

of the loop; this ensures that the implementation of the algorithm is correct. The program prova.f90 in 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!

[name:] DO exit [name]

or [name:] DO END DO [name]

possible forms of "do while": DO if (condition)exit END DOor: 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

Details on Fortran: unbounded loops

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

Radioactive decay:results of numerical

simulation

log(N(t))

N(t)

t

t

Semilog plot of the results of decay simulation for the same decay rate and different initial number of atoms:almost a straight line, but withimportant deviations (stochastic) for small N

Radioactive decay:results of numerical

simulation

numerical simulations: OK on average and for large numbers

A cosa servono i numeri

pseudocasuali?

•flusso di particelle tra due meta‘ di una scatola...

•decadimento radioattivo decay.f90

• random walk segue in dettaglio

•per simulare processi casuali:

Random walks

•come simulare?

•quali proprieta’ interessano?

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

Cosa possiamo calcolare PER OGNI CAMMINATORE dopo N passi ?

: distanza DI UN CAMMINATORE dalla partenza : la sua distanza QUADRATICA dalla partenza

Random walks unidimensionali

Cosa possiamo calcolare PER PIU’ CAMMINATORI dopo N passi ?

: spostamento MEDIO dopo N passi : spostamento QUADRATICO MEDIO dopo N passi

Random walks unidimensionali

ATTENZIONE!!! proporzionale al numero di passi

Piu’ generale:

Random walks unidimensionali

(ma lasciamo da parte per ora ...)

per il singolo camminatore

per piu’ camminatorisi definiscono anche in due o piu’

dimensioni:

: probabilita’ che IN MEDIA sia lo spostamento finale dal punto di partenza dopo N passi

Random walksUn’altra grandezza calcolabile SU PIU’

CAMMINATORI:

In piu’ dimensioni???

RW 1D: simulazioneThe basic algorithm:ix = position of the walkerx_N, x2_N = cumulative quantitiesrnd(N) = sequence of N random numbers

(1 run= 1 particle)

RW 1D: simulazione

rw1d.f90

Proprieta’ statistiche

Un’altra proprieta’ statisticadei random walks (qui 1D)

Random walks unidimensionaliE’ 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 !

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)

Un’altra proprieta’ statisticadei random walks (qui 1D)

Random walks unidimensionali

La tavola di Galton

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.

Simulazione della distribuzione di probabilita’ della posizione finale di un random walk unidimensionale

(con passo variabile)

poche tante

al solito, attenzione alla statistica;

gia’ abbiamo visto per quanto riguarda Δ2

N;vediamo anche per PN:

RW 1D: simulation

RW 1D: simulation

RW 1D: simulation

RW 1D: simulation

Random Walks 2D

Random Walks 2DPer piu’ camminatori:

Random Walks 2DPer piu’ camminatori:

RW 2D: simulazione

rw2d.f90rw2zoom.f90 (stampa posizioni ogni tot passi,

per vedere la caratteristica frattale del cammino)

Random Walks 2D 0 0.0000000 0.0000000 10 0.2242774 3.7794106 20 -1.7333623 1.3218992 30 -1.4481916 -3.1119978 40 -2.2553353 -3.5246484 50 -3.8911035 -6.6665235 60 -3.6508965 -8.0110636.....

0 0.0000000 0.0000000 1 0.6946244 0.7193726 2 0.9359566 1.6898152 3 1.8891419 1.9922019 4 0.9642899 2.3725290 5 0.1308700 2.9251692 6 0.2071800 3.9222534 7 0.9160752 4.6275673 8 0.2856980 3.8512783 9 1.0143363 3.1663797 10 0.2242774 3.7794106 11 -0.7752404 3.8104627 12 -1.7280728 3.5069659 13 -2.0930278 4.4379911 14 -3.0587580 4.1784425 15 -2.0729706 4.0104446 16 -1.8304152 3.0403070 17 -2.2890768 2.1516960 18 -1.7717266 1.2959222 19 -1.1920205 0.4810965 20 -1.7333623 1.3218992 21 -1.5798329 0.3337551.....

if (mod(i,10)==0) then WRITE (...) i,x,yend if

WRITE (...) i,x,y

Random Walks 2Dself-similarity!

plot every 10 steps plot every step

‘pl’: macro for gnuplot for plotting trajectories(suppose column 1 is ‘time’, 2 is x, 3 is y)and check self-similarity:

set term postscript colorset size squareset out '1.ps'p [-20:5][-10:15] '1.dat' u 2:3 w lset out '10.ps'p [-40:20][-10:50] '10.dat' u 2:3 w l, 'contour' u 1:2 w l

Use:gnuplot$ load ‘pl’

Random Walks 2D

A cosa servono?un lungo elenco...

- finanza- polimeri (SAW)- diffusione...

http://solarscience.msfc.nasa.gov/interior.shtml

ma anche per modellare la diffusione dei fotoni attraverso la zona radiativa del sole:

Diffusione: tanti tipi di diffusione....Ad es.: di neutroni all’interno di un reattore nucleare

moto browniano e traiettoria

frattale