Please answer the following questions in complete sentences in a clearly prepared manuscript and submit the solution by the due date on Blackboard (around Sunday, September 9nd, 2018.)
Remember that this is a graduate class. There may be elements of the problem statements that require you to fill in appropriate assumptions. You are also responsible for determining what evidence to include. An answer alone is rarely sufficient, but neither is an overly verbose description required. Use your judgement to focus your discussion on the most interesting pieces. The answer to "should I include 'something' in my solution?" will almost always be: Yes, if you think it helps support your answer.
Please identify anyone, whether or not they are in the class, with whom you discussed your homework. This problem is worth 1 point, but on a multiplicative scale.
Make sure you have included your source-code and prepared your solution according to the most recent Piazza note on homework submissions.
In class, we showed that we could form a linear system of equations in order to determine the expected time that a random walk on the integers spends before it reaches the endpoints and . We can do the same thing to determine the expected length of a game of Candyland!
Recall that the data for our Candyland game is available from the website urls:
The starting cell is given by cell 140 and the ending cell is given by cell 134. Cells 135 and 136 are are bridge cells. Cells 137, 138, and 139 are used to model the sticky cells.
Let be the probability of moving from Candyland cell to Candyland cell given a draw from the desk of cards.
(Warm-up) Show how to identify the game-cell associated with sticky cells 137, 138, and 139 based on the transition matrix . Your answer must include both the answer in the given data and and procedure that would apply more generally.
Use the same type of implicit formulation where is the expected length of a Candyland game starting from cell to derive a linear system of equations. Build and solve this linear system of equations in Julia. For this problem, assume that you use one turn for any cell other than the ending cell . (That is, we are doing no modeling associated with the bridges or sticky states.) What is the expected length of the game starting from the start? (Hint, I get 33.66804940823773). You should provide the linear system in terms of the matrix , but you should not provide explicit entries. Also, presumably, the creators would have made the start state the place of maximum game length. Is this the case or is there another cell or set of cells that result in a longer game?
For this problem, I found it helpful to visualize a solution to make sure it made sense. I used the following line of code to visualize a solution :
using Plots
pyplot(size=(600,600))
scatter(xc,yc,markersize=x, zcolor=x)
gui()
where xc
and yc
are the coordinates from the coords file.
Recall that in class we showed that gave the probability of being in each state after steps. By definition, the expected length of the game is: Develop a computer program to approximately compute this sum. Is there a good way to determine when to stop adding terms? What value do you get?
Now, ideally, the answer to parts 3 and 4 are the same. (Hint: they are.) Use the Neumann series of a matrix to prove that this is the case.
Going back to part 2, we didn't get our model of Candyland quite right because we assumed that each of the bridge and sticky states should be associated with one turn, including a transition to the virtual state . Show how to modify your linear system of equations so that these virtual transitions are no longer included. (Note, there is more than one way to solve this.) Just to be precise, we want to assume the following game rules:
How does this change the expected length of the game? Also, are there still a cell or set of cells that results in a longer game than the actual start?
In this problem, we'll meet one of the most common matrices
studied in numerical linear algebra: the -Laplacian. We
arrive at this matrix by discretizing a partial differential
equation. Poisson's equation is:
where is a continuous function defined over the unit-plane
(i.e. ),
is a continuous function defined over the same region,
and is the Laplacian operator:
Given a function , we want to find a function that satifies this equation. There are many approaches to solve this problem theoeretically and numerically. We'll take a numerical approach here.
Suppose we discretize the function at regular points ,
and where so that we have:
For this discretization, note that
What we've done here is use the approximation:
for both partial terms.
We need this equation to hold at each point .
But note that there are some issues with this equation at the boundary
values (where x=0 or 1, or where y=0 or 1).
For this problem, we'll make it very simple and set:
Now, we'll do what we always do! Turn this into some type of matrix equation!
Let be an matrix that we'll index from zero instead of one: where . At this point, we are nearly done. What we are going to do is turn Poisson's equation into a linear system. This will be somewhat like how we turned image resampling into a matrix vector equation in the first homework.
In order to write as a vector, we'll keep the convention from last time: Let be the vector of elements here. Note that our approximation to , just involved a linear combination of the elements of . This means we have a linear system: where the rows of and correspond to equations of the form:
Let . Write down the linear equation for including all the boundary conditions. Note that you can encode the boundary conditions by adding a row of where: .
Write a Julia or Matlab/Python code to construct a sparse matrix and vector when and . Here's some pseudo-code to help out:
function laplacian(n::Integer, f::Function)
N = (n+1)^2
nz = <fill-in>
I = zeros(Int,nz)
J = zeros(Int,nz)
V = zeros(nz)
fvec = zeros(N)
# the transpose mirrors the row indexing we had before.
G = reshape(1:N, n+1, n+1)' # index map, like we saw before;
h = 1.0/(n)
index = 1
for i=0:n
for j=0:n
row = G[i+1,j+1]
if i==0 || j == 0 || i == n || j == n
# we are on a boudnary
fvec[row] = 0.0
# fill in entries in I,J,V and update index
else
fvec[row] = f(i*h, j*h)*h^2
# fill in entries in I,J,V and update index
end
end
end
A = sparse(I,J,V,N,N)
return A, fvec
end
Solve for using Julia's or Matlab's backslash solver, and show the
result using the mesh
function (Matlab) or surface
function (Plots.jl in Julia).
We can use the Neumann series to turn an inverse into an infinite summation. What happens if we try and use that approach to solve this linear system of equations? Does it work or not?
Write working code for the following operations for a matrix
given by Compressed Sparse Column arrays: colptr
, rowval
, nzval
along with dimensions m
and n
as in the Julia sparse matrix
storage.
Sparse matrix-transpose multiplication by a vector
""" Returns y = A'*x where A is given by the CSC arrays
colptr, rowval, nzval, m, n and x is the vector. """
function csc_transpose_matvec(colptr, rowval, nzval, m, n, x)
end
Row-inner-product
""" Returns = A[i,:]*x where A is given by the CSC arrays
colptr, rowval, nzval, m, n and x is the vector. """
function csc_row_projection(colptr, rowval, nzval, m, n, i, x)
end
Column inner-product
""" Returns = A[:,i]'*x where A is given by the CSC arrays
colptr, rowval, nzval, m, n and x is the vector. """
function csc_column_projection(colptr, rowval, nzval, m, n, i, x)
end