Computational Techniques in Theoretical Physics
Tutorial 3

• Generate a random number x according to the following distributions:

(a)                    2x exp(-x2) ,    x > 0
p(x) =  {
0 ,    x < 0

Solution:

F(x) = int -infinity to x p(x) dx

= int 0 to x 2x exp{-x2 } dx
= int 0 to x exp{-x2 } dx2
= 1 - exp(-x2)

x = F-1(xi) =  [ - ln( 1 - xi ) ]1/2 = [ - ln xi ]1/2

(b)                    exp(-x)/[ 1 + exp(-x) ] ,   x > 0
p(x) = {
0 ,    x < 0

Solution:

F(x) = int -infinity to x p(x) dx

= int 0 to x exp(-x) / [ 1 + exp(-x) ]  dx
= - int 0 to x 1/[ 1 + exp(-x) ] d exp(-x)
=  ln 2 - ln [1 + exp(-x) ]

x = F-1(xi) =  - ln( 2e-xi - 1)

(c)                     x exp(-x)   ,    x > 0
p(x) = {
0  ,       x < 0

Solution:

F(x) = int -infinity to x p(x) dx

= int 0 to x  x exp{-x } dx
= - int 0 to x x d exp{-x }
= - x exp{-x } +  int 0 to x  exp{-x } dx
= 1 - ( x + 1 )exp(-x)

xi = 1 - ( x + 1 ) exp(-x)

This function can not be solved analytically. It needs to
solved numerically. Write a C-code to do this.

• Write down the Monte Carlo estimator and the error for the following integrals:

(a)   int 0 to infinity  2x exp(-x2)

Solution:

Monte Carlo estimator:

G = int   g(X) p(X) dX

GN =  sum i=1 to N g(Xi) / N

Monte Carlo error:

sigma2 ~

N/(N - 1){1/N sum i=1 to N g2(Xi) -
(1/N sum i=1 to N g(Xi))2}

g(x) = 2x

GN =  sum i=1 to N 2 Xi / N

sigma2 ~
N/(N - 1){1/N sum i=1 to N 4 Xi2 -
(1/N sum i=1 to N 2Xi)2}

(b)   int 0 to infinity exp(-x)/[ 1 + exp(-x) ]

Solution:
g(x) = 1/ [ 1 + exp(-x) ]

GN =  sum i=1 to N  [ 1 + exp(-Xi) ]-1 / N

sigma2 ~
N/(N - 1){1/N sum i=1 to N [ 1+ exp(-Xi) ]-2 -
(1/N sum i=1 to N [ 1 + exp(-Xi) ]-1)2}

(c)   int 0 to infinity x exp(-x)

Solution:

Monte Carlo estimator:

g(x) = x

GN =  sum i=1 to N  Xi / N

sigma2 ~
N/(N - 1){1/N sum i=1 to N  Xi2 -
(1/N sum i=1 to N Xi)2}

• The transition probability W(i to j)  between 3 states is:
W(1 to 1) = W(1 to 3) = 0, W(1 to 2) = 1

W(2 to 1) = 1/6, W(2 to 2) = 1/2,
W(2 to 3) = 1/3

W(3 to 1)=0, W(3 to 2)=2/3, W(3 to 3) = 1/3

(a)  Show that the solution for the probabilities is:

P1 = 1/10, P2 = 6/10, P3=3/10

(b)  Use P0=(1,0,0) as the initial distribution,
how fast the system converges to the equilibrium
probability distribution?

Solution:

(a)

0       1       0

W = 1/6     1/2     1/3

0     2/3     1/3

P = ( P1, P2, P3)

P = P W

P1 =  0  + 1/6P2 + 0
P2 = P1 + 1/2 P2 + 2/3 P3
P3 = 0 + 1/3 P2 + 1/3 P3

P1 + P2 + P3 = 1

which gives:  P1 = 1/10, P2 = 6/10, P3=3/10

(b)

P0 = (1, 0, 0)

P1 = P0 W = (0, 1, 0)

P2 = P1 W = (1/6, 1/2, 1/6)

P3 = P2 W = (1/12, 19/36, 2/9) = (0.08, 0.52, 0.22)

P4 = P3 W = ....

• Consider a situation where we have only three possible states 1, 2, and 3 (e.g. discrete quantum states). The transition probability is given as
W(1 to 2) = 1, W(2 to 1) = W(2 to 3) = W(3 to 2) = W(3 to 3) = 1/2, and W(1 to 3) = W(3 to 1) = W(1 to 1) = W(2 to 2) = 0.

(a) Is the Markov chain specified by the above transition
probabilities ergodic?

(b) Calculate the equilibrium probability distribution Pi,
where i = 1,2,,3.

(c) Does the equilibrium distribution satisfy detailed balance
condition?

Solution:

(a)

0        1        0

W =   1/2     0        1/2

0       1/2     1/2

1/2      0        1/2

W2 =   0       3/4       1/4

1/4    1/4       1/2

0       3/4      1/4

W3 =    3/8    1/8      1/2

1/8    1/2       3/8

W4  = nonzero matrix.

Thus the Markov chain is ergodic.

(b)
Equilibrium distribution P=(P1, P2, P3) satisfies the equation:

P = P W

which gives:

P1 =  0 + 1/2 P2 + 0
P2 = P1 + 0 + 1/2 P3
P3 = 0 + 1/2 P2 + 1/2 P3

In addition there is the condition:

P1+P2+P3 = 1

From these equations P can be solved.

(c)
Detailed balance condition:

P(X') W(X' to X) = P(X) W(X to X')

Check whether P1, P2, P3 satisfies:

Pi W(i, j) = Pj W(j, i)

If yes, then P satisfies the detailed balance condition.