Phys690: Homework Assignment 2

2/19/07

==> Due on Fr 3/2 by 4pm in my mailbox <==

(1). Monte Carlo sampling of an exponential distribution.

The exponential probability density function has frequent applications in radiation transport, queuing theory, etc. As we learned in class, it is straightforward to sample such a distribution, p(x)= exp(-x/lambda)/lambda, (x>0), i.e., to produce x's that are distributed according to the probability density p(x). Here we first write an exponential random number generator that does this sampling, and then measure properties of the sampled x's and compare them with analytic results from p(x).

Below are some Maple codes. Some of what the codes do can be done simply by Maple library packages and made much more fancy. Exploring this possibility is good and encouraged, after you have understood the provided codes. These codes can be run through main, which reads in the other two. It looks like:

#------------------------------- 0 -----------------------------------
 # read in the other 2 programs, run & mkhist:
read run:
read mkhist:
with(plots):
#------------------------------- 1 -----------------------------------
 # set parameters to execute programs:
n_mea := 10;             # number of measurements
 #       ^^ <-- MODIFY to change accuracy of final results
n_smpls_per_mea := 100;  # number of samples in each measurement
 #                 ^^^ <-- MODIFY to change statistics of measurements
lambda := 1.5;           # parameter lambda in exponential PDF
 #        ^^^ <-- MODIFY to change PDF being sampled
#------------------------------- 2 ------------------------------------
 # execute programs:
run(lambda, n_mea, n_smpls_per_mea):
mkhist(lambda, n_mea, n_smpls_per_mea):
#------------------------------- 3 ------------------------------------
 # view results:
 #    plot of histogram of sampled distribution vs analytical PDF:
display({pex,phi},title=`Sample (pts) vs ext (line) PDF`);
 #    computed <x> from samples and its statistical error:
x_avg := sum(xacc[m], m=1..n_mea)/n_mea;
x_err := sqrt((sum(xacc[m]^2, m=1..n_mea)/n_mea-x_avg^2)/(n_mea-1));

Download run and mkhist. Please note that there are two blanks you need to fill to make the codes work correctly. One is in the expression for the exponential random number generator ranexp in run; the other is the histogram normalization in mkhist.

Fill in the two blanks. Include the answer and derivations in your write-up.
Run the code with the the current choice of parameters in block 1 in main. How does the running time depend on these parameters? Use this as a guideline when changing parameters in the next step.
Now change the three parameters, run the code (i.e., repeat blocks 1, 2, and 3) and justify your results. Do this and refer back to the three source files until you understand the structure. With a fixed value of lambda=1.5, do a few runs with different choices of the other 2 parameters to demonstrate the sqrt dependence of the error of <x> on both parameters.
Evaluate the integral INT_{0, infty} [ x^3 exp(-2*x) dx ] by Monte Carlo. As always, include error estimate in your result. Compare your result with the exact answer.

run:

 # ranxi() gives a random # uniform on (0,1)
 # We're trusting maple's generic RNG here, which may not be safe!   
ranxi := () -> (10e-13)*rand();            

 # ranexp(lambda) should give a random # distributed according to the 
 #    exponential probability density function exp[-x/lambda]/lambda
 #    where x is on ( 0, +infinity ) and lambda ( >0 ) is a parameter.   
ranexp := (lambda) -> ranxi()     ;
 #                    ^^^^^^^^^^^^ <-- MODIFY -- fill in correct form
####
 # program that calls ranexp(lambda) to generate a histogram of 
 #    distribution and also measures the first moment 
run := proc(lambda, n_mea, n_smpls_per_mea)

global xacc, hist, bin_range, bin_size, n_bins, bin_size_inv, midpt, bin;
local bi, st, m, i, x;
 # -- define parameters and procedures for histogram:
bin_range := 6. *lambda;                #NOTE: prefactor can be changed
bin_size := 0.1 *lambda;                #NOTE: prefactor can be changed
n_bins := trunc(bin_range/bin_size)+1;  #last bin to store x > bin_range
bin_size_inv := 1./bin_size;
hist := array(1..n_bins);
for bi to n_bins do hist[bi] := 0 od:
bin  := proc(x)    
    global n_bins, bin_size_inv; 
    min(trunc(x*bin_size_inv)+1,n_bins) 
end:                                    # bin(x) gives which bin x is in.  
midpt := proc(i) global bin_size; (i-0.5)*bin_size end:
                                        # midpt(i) gives mid-pt of bin i
####
 # -- xacc[m] stores the m^th measurement of 
xacc := array(1..n_mea);
 # -- n_mea measurements [first do loop (with m)], 
 #    EACH with n_smpls_per_mea samples [second do loop (with i)] :
for m from 1 to n_mea do
    xacc[m] := 0:
    for i to n_smpls_per_mea do
        x := ranexp(lambda):
        xacc[m] := xacc[m] + x:
        bi := bin(x):
        hist[bi] := hist[bi] +1
    od;
    xacc[m] := xacc[m]/n_smpls_per_mea:
od;
 #
end:

mkhist:

 # program to fix up histogram (normalization, etc) and 
 #  prepare plots with both histogram and the analytical PDF
 # called AFTER a call to 'run'
mkhist := proc(lambda, n_mea, n_smpls_per_mea)
 #
global xacc, hist, bin_range, bin_size, n_bins, midpt, phi, pex;
local histnorm_inv, bi, hi, exact_pdf;
 #
 # normalize histogram (histnorm_inv is INVERSE of normalization):
histnorm_inv := 1.                        :
 #              ^^^^^^^^^^^^^^^^^^^^^^^^^^ <-- MODIFY -- 
 #                                                  fill in correct form 
 # Note histnorm_inv appears in the following line:
for bi to n_bins do hist[bi] := hist[bi]*histnorm_inv od:
 # plot histogram:    
hi := [[midpt(n),hist[n]] $n=1..n_bins]:
phi := plot(hi,x=0..bin_range+bin_size,style=point,symbol=circle,color=red):
 # plot exact probability density function
exact_pdf := proc(x)
    global lambda;
    exp(-x/lambda)/lambda:
end:
pex := plot(exact_pdf(x),x=0..bin_range+bin_size,style=line,color=blue):
 # end:

(2). Random walks with traps, and Green's function.

We consider again the one-dimensional version of the drunkard problem discussed in class. This time he starts out at position x (e.g., x=10), and takes unit steps left or right with equal probability. Suppose at 0 and a (a>x>0, e.g., a=30) are two bars. When he reaches either one, he is trapped forever.

We want to determine the average time t (the mean number of steps) before the drunkard is trapped. This is a diffusion problem and the 'time' t is called the first passage time. Do a simulation to compute t. (Again note that only the average of t has a meaning and it is that average that we want to compute.) Then run your simulation for several different x and try to relate t (as a function of x) to x*(a-x). Does the result make sense? Use your favorite programming language. The structure of the code for problem (1) may be helpful.
We now use your new program to study the distribution of the drunkard (on (0,a), excluding end points) as a function of time (i.e., number of steps). (This is related to the so-called Green's function which is an important concept in diffusion, E&M, wave propagation, quantum mechanics, etc.) Given a drunkard that starts out at x=10, what is his distribution after 20 steps. Simulate the random walk enough times and accumulate the distribution, until you have a satisfactory histogram. Note that we are only interested in where the drunkard is after exactly 20 steps.
Before he is trapped, how many times on the average does the drunkard visit each position? Note that here we are interested in the entire path history of the drunkard (as opposed to an instantaneous position in part 2). Again, simulate the random walk enough times until you have a satisfactory histogram.