The Monte Carlo Method: an Introduction Detlev Reiter

Vorlesung HHU Düsseldorf, WS 07/08 March 2008

Research Centre Jülich (FZJ) D 52425 Jülich http://www.fz-juelich.de

e-mail: [email protected]

Tel.: 02461 / 61-5841

There are two dominant methods of simulation for complex many particle systems • • • • 1) Molecular Dynamics Solve the classical equations of motion from mechanics.

Particles interact via a given interaction potential.

Deterministic behaviour (within numerical precision).

Find temporal evolution.

• • 2) Monte Carlo Simulation Find mean values (expectation values) of some system components.

Random behaviour from given probability distribution laws.

The Monte Carlo technique is a very far spread technique, because it is not limited to systems of particles.

This lecture

•Brief introduction: simulation •What is the Monte Carlo Method •Random number generation •Integration by Monte Carlo Tomorrow: one (of many) particular application: •particle transport by Monte Carlo

ASDEX-UPDRADE (IPP Garching)

4

Monte Carlo particle trajectories, ions and neutral particles

Trilateral Euregio Cluster TEC Inst it ut f ür Plasmaphysik Assoziat ion EURATOM-Forschungszent rum Jülich

Basic principle of the Monte Carlo method

• The task: calculate (estimate) a number

I

(one number only. Not an entire functional dependence).

Historic example: A dull way to calculate p – Numerically: look for an appropriate convergent series and evaluate this approximately – by Monte Carlo: look for a stochastic model (i.e.: ( W , s ,

p , X

): probability space with random variable

X

) Example: throw a needle an a sheet with equidistant parallel stripes. Distance between stripes: D, length of needle: L, L

First application of Monte Carlo Method

The needle experiment of Comte de Buffon, 1733 (french biologist, 1707-1788)

What is the probability p, that a needle (length L), which randomly falls on a sheet, crosses one of the lines (distance D)?

(N trials, n „hits“)

Y t =1, if crossing, Y t =0 else, then

Today: Using a computer to generate random events: We need to be able to generate random numbers X with any given probability function f(x), or a given cumulative distribution F(x) .

1) Uniformly distributed random numbers 2) General random numbers : can be obtained from a sequence of independent uniform random numbers

Random number generation f(x) 1/(b-a) a b

We will see next: Any continuous distribution can be generated from uniform random numbers on [0,1] Any discrete distribution can be generated from uniform random numbers on [0,1] Hence: Any given distribution can be generated from uniform random numbers on [0,1]

Strategy: try to transform F to another distribution, such that inverse of new F is explicitly known.

Example: Normal (Gaussian) distribution Cumulative distr. function Inverse cumul. distr. fct.

best format of storing distributions for Monte Carlo applications: „Inverse cumulative distribution function F -1 (x)“, x uniform [0,1]

Exercise (and most important example:) Generate random numbers from a Gaussian.

Let X, Y two independent Gaussian random numbers.

Transform to polar coordiantes (Jacobian!) R, Φ Sample Φ (trivial, it is uniform on 2π) Apply inversion method for R Transform sampled Φ, R back to X, Y.

This is a pair of Gaussians. ( Box-Muller Method )

Exponential distribution by „inversion“ Note: Z and 1-Z have same distrib.

(see tomorrow )

Cauchy: e.g.: natural Line broadening

( stepwise constant, with steps at points T)

Rejection

y=f(x) f(x): distribution density enclosing rectangle

Accept z, take x=z Reject z y uniform

sample x from f(x)

X

z, uniform

NEXT: Any Monte Carlo estimate can be regarded as a mean value, i.e. an integral (or sum) over a given probability distribution, ususally in a high dimensional space (e.g. of random walks….) Generic Monte Carlo: Integration Hence: How does Monte Carlo integration work?

Hit or Miss

I: unknown area known area

hit

f(x)

x 2 uniform

I = ∫ f(x) dx

x 1 , uniform

X

miss

Suggestion: try again with previous example from dull and crude Monte Carlo

Outlook: next lecture (tomorrow)

END