In [1]:
## In this set of demo's we are going to study
# random walks on the integers.
In [2]:
## Uniform random walks
# In this problem, we can move around the integers.
# At each step, we move -1 or +1 with prob 1/2 each.
# We start at the point 0 and stop when
# we get to +6 or - 4.
# How long does this take to stop?

# In this problem, we'll show that the answer is 24.

# A simulation is easy to code in Julia
function simlength()
  S = 0
  len = 0
  while S != -4 && S != 6
    len += 1
    S += rand(-1:2:1)
  end
  return len
end
Out[2]:
simlength (generic function with 1 method)
In [3]:
## Now we need to run this a large number of times.
# in-class: write the code to simulate this
# in-class: typo on function to show scoping rules
N = 10^6
avg = 0
for i=1:N
  avg += simlength()/N
end
avg
# this fails because of Julia's weird scoping rules with Global :(
Out[3]:
24.01852000002423
In [4]:
##
avg = 5
for i=1:10
  avg = 6
end
In [5]:
## You get the expected behavior in a function
function mysim(N)
  avg = 0
  for i=1:N
    avg += simlength()/N
  end
  return avg
end
mysim(10^6)
Out[5]:
23.98940400002502
In [6]:
## Here's a fancy Julia way.
# run simlength 10^6 times
N = 10^6
expectedlen = mapreduce(x -> simlength()/N, +, 1:N)
Out[6]:
24.00348000000004
In [7]:
## Now, let's derive the following system.
using LinearAlgebra
A = collect(SymTridiagonal(ones(11), -0.5*ones(10)))
Out[7]:
11×11 Array{Float64,2}:
  1.0  -0.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
 -0.5   1.0  -0.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
  0.0  -0.5   1.0  -0.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0  -0.5   1.0  -0.5   0.0   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0  -0.5   1.0  -0.5   0.0   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0  -0.5   1.0  -0.5   0.0   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0  -0.5   1.0  -0.5   0.0   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0  -0.5   1.0  -0.5   0.0   0.0
  0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.5   1.0  -0.5   0.0
  0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.5   1.0  -0.5
  0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.5   1.0
In [8]:
##
A[1,2] = 0
A[end,end-1] = 0
b = [0;ones(9);0]
Out[8]:
11-element Array{Float64,1}:
 0.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0
 0.0
In [9]:
##
x = A\b
Out[9]:
11-element Array{Float64,1}:
  0.0              
  9.0              
 16.0              
 21.0              
 24.000000000000004
 25.000000000000004
 24.000000000000004
 21.0              
 16.0              
  9.0              
  0.0
In [10]:
## Convert A as a sparse matrix.
using SparseArrays
B = sparse(A)
Out[10]:
11×11 SparseMatrixCSC{Float64,Int64} with 29 stored entries:
  [1 ,  1]  =  1.0
  [2 ,  1]  =  -0.5
  [2 ,  2]  =  1.0
  [3 ,  2]  =  -0.5
  [2 ,  3]  =  -0.5
  [3 ,  3]  =  1.0
  [4 ,  3]  =  -0.5
  [3 ,  4]  =  -0.5
  [4 ,  4]  =  1.0
  [5 ,  4]  =  -0.5
  [4 ,  5]  =  -0.5
  [5 ,  5]  =  1.0
  ⋮
  [7 ,  7]  =  1.0
  [8 ,  7]  =  -0.5
  [7 ,  8]  =  -0.5
  [8 ,  8]  =  1.0
  [9 ,  8]  =  -0.5
  [8 ,  9]  =  -0.5
  [9 ,  9]  =  1.0
  [10,  9]  =  -0.5
  [9 , 10]  =  -0.5
  [10, 10]  =  1.0
  [10, 11]  =  -0.5
  [11, 11]  =  1.0
In [11]:
## Create A as a sparse matrix.
# we list all the non-zero entries in A.
function build_A()
  nz = 9*3 + 2
  I = zeros(Int, nz)
  J = zeros(Int, nz)
  V = zeros(nz)
  # the arrays I,J,V satisfy:
  # A[I[i],J[i]] = V[i]
  index = 1
  I[index] = 1
  J[index] = 1
  V[index] = 1.0
  index += 1
  for i=1:9
    I[index] = i+1
    J[index] = i
    V[index] = -0.5
    index += 1

    I[index] = i+1
    J[index] = i+1
    V[index] = 1.0
    index += 1

    I[index] = i+1
    J[index] = i+2
    V[index] = -0.5
    index += 1
  end

  I[index] = 11
  J[index] = 11
  V[index] = 1.0
  index += 1

  return sparse(I,J,V,11,11)
end
C = build_A()
Out[11]:
11×11 SparseMatrixCSC{Float64,Int64} with 29 stored entries:
  [1 ,  1]  =  1.0
  [2 ,  1]  =  -0.5
  [2 ,  2]  =  1.0
  [3 ,  2]  =  -0.5
  [2 ,  3]  =  -0.5
  [3 ,  3]  =  1.0
  [4 ,  3]  =  -0.5
  [3 ,  4]  =  -0.5
  [4 ,  4]  =  1.0
  [5 ,  4]  =  -0.5
  [4 ,  5]  =  -0.5
  [5 ,  5]  =  1.0
  ⋮
  [7 ,  7]  =  1.0
  [8 ,  7]  =  -0.5
  [7 ,  8]  =  -0.5
  [8 ,  8]  =  1.0
  [9 ,  8]  =  -0.5
  [8 ,  9]  =  -0.5
  [9 ,  9]  =  1.0
  [10,  9]  =  -0.5
  [9 , 10]  =  -0.5
  [10, 10]  =  1.0
  [10, 11]  =  -0.5
  [11, 11]  =  1.0