add code for chapter 21

Author: vtjnash

demos/book/21/airyeigmax.jl

@@ -0,0 +1,51 @@
+#airyeigmax.jl
+#Algorithm <21.1> of Random Eigenvalues by Alan Edelman
+
+#Experiment:    Largets eigenvalue of a Stochastic Airy Operator
+#Plot:          Histogram of the largest eigenvalues
+#Theory:        The Tracy-Widom Law
+
+## Parameters   # Include the most interesting configurable values at the top, with brief descriptions
+t = 10000       # trials
+n = 1e9         # level of discretization
+beta = 2
+h = n^(-1/3)    # h serves as dx
+
+function airyeigmax_experiment(t,n,beta,h) # wrap the experimental logic in a function to enable faster JIT
+    ## Experiment
+    v = zeros(t)    # samples
+    x = 0:h:10      # discretization of x
+    N = length(x)
+
+    #Parallel experiment
+    b = (1 / h^2) * ones(N-1)
+    v = pmap(1:t) do i
+        ## discretize stochastic airy operator
+        # discretize airy operator
+        a = -(2 / h ^ 2) * ones(N);  # differential operator: d^2 / dx^2
+        a -= x                       # d^2 / dx^2 - x
+        # add the stochastic part
+        dW = randn(N) * sqrt(h)
+        a += (2 / sqrt(beta)) * dW / h
+        ## calculate the largest eigenvalue of tridiagonal matrix T
+        # diagonal of T: a
+        # subdiagonal of T: b
+        return eigmax(SymTridiagonal(a,b))
+    end
+
+    binsize = 1/6
+    return hist(v, -6:binsize:6)
+end
+grid, count = airyeigmax_experiment(t,n,beta,h) #run the experiment, making the global variables local for speed
+
+## Plot
+include("../1/tracywidom.jl")
+h = Histogram(count/(sum(count)*step(grid)), step(grid))    # add the histogram of data
+h.x0 = first(grid)
+add(p, h)
+
+if isinteractive()
+    Winston.display(p)
+else
+    file(p, "<filename>.png")
+end

demos/book/4/mpexperiment.jl

@@ -11,7 +11,7 @@ r = 1           # aspect ratio
 n = 100         # matrix column size
 dx = 0.1        # binsize
 
-function mp_experiment()
+function mp_experiment(n,r,t,dx)
     m = iround(n/r)
     ## Experiment
     v = Float64[]          # eigenvalue samples
@@ -28,7 +28,7 @@ function mp_experiment()
     y = real(sqrt((x.^2-a^2).*(2^2-x.^2) + 0im)./(pi*x*r))
     return (hist(v, x), (x, y))
 end
-((grid, count), (x,y)) = mp_experiment()
+((grid, count), (x,y)) = mp_experiment(n,r,t,dx)
 
 ## Plot
 using Winston

demos/book/README

@@ -2,7 +2,7 @@ The codes in this folder are translations of MATLAB code by Alan Edelman in his
 
 It is intended that these codes can be run in batch mode to generate the sample figures as output files, but that the codes are most useful when used interactively in the exploration of random matrix theory. For this reason, when run from the Julia REPL, the plots will be output to the screen.
 
-Additional codes should use the following template (modified from book/4/mpexperiment.jl). While some experiments will need to be single-threaded, many codes in Random Matrix Theory are Monte Carlo simulations which are trivial to parallelize using the pmap function. To make use of the parallelism this provides, either start Julia with the `-p N` flag or run `addprocs(N-1)` at the REPL (where N>=3 and equals CPU_CORES+1 optimally). When doing Julia-level parallelism, it is typically best to disable OpenBLAS parallelism by calling `Base.openblas_set_num_threads(1)`.
+Additional codes should use the following template (modified from book/4/mpexperiment.jl). While some experiments will need to be single-threaded, many codes in Random Matrix Theory are Monte Carlo simulations which are trivial to parallelize using the pmap function. To make use of the parallelism this provides, either start Julia with the `-p N` flag or run `addprocs(N-1)` at the REPL (where N>=3 and equals CPU_CORES+1 optimally). When doing Julia-level parallelism, it is typically best to disable OpenBLAS parallelism by calling `Base.openblas_set_num_threads(1)`. Run code at the REPL using `include("<filename.jl>")` so that it is run only on the local processor (`require` is run on all processors).
 
 ```julia
 #<filename>.jl
@@ -17,7 +17,7 @@ t = 10000       # trials
 n = 100         # matrix column size
 dx = 0.1        # binsize
 
-function mp_experiment(t,n,dx) # wrap the experimental logic in a function to enable faster JIT
+function an_experiment(t,n,dx) # wrap the experimental logic in a function to enable faster JIT
     ## Experiment
 
     #Single threaded experiment
@@ -35,7 +35,7 @@ function mp_experiment(t,n,dx) # wrap the experimental logic in a function to en
     y = x .^ 2
     return (hist(v, x), (x, y))
 end
-((grid, count), (x,y)) = mp_experiment(t,n,dx) #run the experiment, making the global variables local for speed
+((grid, count), (x,y)) = an_experiment(t,n,dx) #run the experiment, making the global variables local for speed
 
 ## Plot
 using Winston