http://www.laureescientifiche.units.it/
La simulazione numerica in fisica:
approcci stocastici
(aggiornamento e formazione insegnanti;
corso CIRD - CP e M IDIFO3: Lab_A I15)
G. Pastore e M. Peressi
Universita’ di Trieste
a.a. 2010/11
14 marzo 2011
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
Numeri casuali
Una sequenza di numeri casuali e’ una sequenza
di numeri che sembrano impredicibili ma che
hanno ben definite proprieta’ statistiche
Generazione di numeri
casuali
metodi manuali:
- lancio di dadi, di monete...
QuickTi me™ and a
GIF decompressor
are needed to see t his pict ure.
QuickTi me™ and a
GIF decompressor
are needed to see t his pict ure.
Quic kT i me™ and a
GIF dec ompres s or
are needed t o s ee thi s pi c ture.
- 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
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
<r>
0.65
0.75
0.85
0.95
Numeri pseudocasuali
e 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)
loop
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:
(pseudo)random numbers
generation:
example I1 - “Linear congruential method”
(Lehemer, 1948)
In+1 = (a In + c) mod m
Limits of the algorithm:
Numeri pseudocasuali
e problemi...
correlazione in iperpiani (routine IBM)
More subtle limits of some algorithms...
xi, p(xi)
(xi, xi+1)
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
Artigianale e’ bello, ma non troppo:
correlazioni con LCM e M=256
Quanti numeri generati? Quante coppie?
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
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
A cosa servono
i numeri
pseudocasuali?
• per calcolare integrali
definiti...
4 x Area settore circ.
Area quadrato
PiMonteCarlo.java
Pi.f90
# punti nel settore circ.
# punti nel quadrato
in onore del
“pi-day” ;-)
A cosa servono
i numeri
pseudocasuali?
• per simulare processi
casuali:
•
flusso di particelle tra due meta‘ di una
scatola: box.f90
• ..
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)
(accade...)
else
if (r <= prob) then
(non accade...)
end if
Macroscopic systems
towards
equilibrium
A simple example: non-interacting classical particles in
a box
(gas diffusion)
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.
N left(t): number of particles present at time t in the left side
Given N left(0), what is N left(t) ?
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)=?
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
Approach to equilibrium II
initial N(left)=1000
How to reduce fluctuations?
- more particles
- average over many simulation runs
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
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 most “random”!
Equilibrium =
Maximum number of
possible microstates =
Maximum entropy
A cosa servono
i numeri
pseudocasuali?
• per simulare processi
casuali:
•
flusso di particelle tra due meta‘ di una
scatola...
• decadimento radioattivo
•
...
decay.f90
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):
Radioactive decay:
numerical simulation
Note: unbounded loop
DO
DO i = 1, nleft
! loop on time
! 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 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!
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
Radioactive decay:
N(t)
results of numerical
simulation
plot of the results of decay
simulation (N vs t)
with N=1000
t
N(t) ~ N0 exp(a t)
log(N(t))
semilog plot (log(N) vs t)
=> log(N(t)) = log N0 + a t
=> slope is a
t
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
A cosa servono
i numeri
pseudocasuali?
• per simulare processi
casuali:
•
flusso di particelle tra due meta‘ di una
scatola...
• decadimento radioattivo
•
random walk
decay.f90
segue in dettaglio
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
Random walks
unidimensionali
Cosa possiamo calcolare
PER OGNI CAMMINATORE dopo N passi ?
: distanza DI UN CAMMINATORE dalla partenza
: la sua distanza QUADRATICA dalla partenza
Random walks
Cosa possiamo calcolare
unidimensionali
PER PIU’ CAMMINATORI
dopo N passi ?
: spostamento MEDIO dopo N passi
: spostamento QUADRATICO MEDIO dopo N
passi
ATTENZIONE!!! proporzionale al numero di
passi
Random walks
unidimensionali
Piu’ generale:
Random
walks
Un’altra grandezza calcolabile SU PIU’
: probabilita’CAMMINATORI:
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:
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)
RW 1D: simulazione
rw1d.f90
Proprieta’ statistiche
dei random walks (qui 1D)
Random walks
E’ facile calcolare analiticamente la probabilita’
unidimensionali
per un random
walk di N passi di finire in un
punto x:
si tratta di un triangolo di Pascal con degli zeri
!
dei random walks (qui 1D)
Random
walks
Per lunghi cammini (N grande) la probabilita’
che sia lo spostamento finale dal punto di partenza
unidimensionali
dopo N passi ha una distribuzione simmetrica rispetto
alla media
e con “larghezza”
standard”,SD)
(“deviazione
La
tavola
di
Galton
Simulazione della distribuzione di probabilita’ della
posizione finale di un random walk unidimensionale (con
passoposizione
variabile)
Biglie che partono da una stessa
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.
poche
tante
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.
al solito, attenzione
alla statistica;
gia’ abbiamo visto per quanto riguarda
Δ2N;
vediamo anche per PN:
RW 1D: simulation
RW 1D: simulation
RW 1D: simulation
RW 1D: simulation
Random Walks 2D
Random Walks 2D
Per piu’ camminatori:
Random Walks 2D
Per piu’ camminatori:
RW 2D: simulazione
rw2d.f90
rw2zoom.f90 (stampa posizioni ogni tot passi,
per vedere la caratteristica frattale del cammino)
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
if (mod(i,10)==0)
.....
1
0.6946244
2
0.9359566
3
1.8891419
4
0.9642899
5
0.1308700
6
0.2071800
7
0.9160752
8
0.2856980
9
1.0143363
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.7193726
1.6898152
1.9922019
2.3725290
2.9251692
3.9222534
4.6275673
3.8512783
3.1663797
Random Walks 2D
then
WRITE (...) i,x,y
end if
WRITE (...) i,x,y
Random Walks 2D
self-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 color
set size square
set out '1.ps'
p [-20:5][-10:15] '1.dat' u 2:3 w l
set 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...
Diffusione: tanti tipi di diffusione....
Ad es.: di neutroni all’interno di un reattore nucleare
ma anche per modellare la diffusione dei fotoni attraverso la zona
radiativa del sole:
http://solarscience.msfc.nasa.gov/interior.shtml
moto browniano
e traiettoria
frattale