I just finished my upper divisional linear algebra course (I took a lower divisional one back in high school) and wanted more practice with applications that could be solved using matrices.
This problem is just a random walks problem. I created a $900\times 900$ Markov matrix $A$, one row/column for each square in the grid. Each $a_{i,j}$ entry in this matrix represent the probability a flea will jump from location $i$ to location $j$. Therefore, the transition matrix for this problem is $A^n$ where $n$ is the number of times the fleas move positions.
I struggled a bit to figure out what to do with this matrix once I got it. I read several articles about random walks and though it helped, none of them discussed situations where there were multiple walks going on at once.
To better understand the problem on my own, I then implemented a brute force method. This did two things. It helped me see the relationship between matrix $A$ and the number of empty squares. It also helped me write several tests.
From here I was able to take these tests back to my matrix implementation and get the final solution. For each square on the grid, I needed to find the probability that it will not have a flea on it. Then, add these probabilities together.
I used NumPy and it was blazingly fast. The speed in which it found $A$ raised to any power was nearly constant. For that reason, I was even able to calculate the expected number of empty squares after $1{,}000{,}000$ rings of the bell. After this many iterations, the system begins to find a steady state. A cool problem would have been to find the expected number of empty squares after an infinite number of bell rings. That would have been fun because then I could have gotten practice with eigenvectors and eigenvalues!
Runs in about 0.6 seconds with Python 3.9.6.
import numpy
def markov_matrix(grid_size, bells):
def neighbors_of(m):
row, col = divmod(m, grid_size)
neighbors = []
if row != 0:
neighbors.append((row - 1, col))
if row != grid_size - 1:
neighbors.append((row + 1, col))
if col != 0:
neighbors.append((row, col - 1))
if col != grid_size - 1:
neighbors.append((row, col + 1))
return neighbors
matrix = [[0 for _ in range(grid_size**2)] for _ in range(grid_size**2)]
for from_n in range(grid_size**2):
neighbors = neighbors_of(from_n)
prob = 1 / len(neighbors)
for neighbor_i, neighbor_j in neighbors:
to_n = grid_size*neighbor_i + neighbor_j
matrix[to_n][from_n] = prob
return numpy.linalg.matrix_power(matrix, bells)
def solve(grid_size, bells):
markov, zeros = markov_matrix(grid_size, bells), 0
for row in markov:
prob = 1
for col in row:
prob -= prob * col
zeros += prob
return round(zeros, 6)
if __name__ == '__main__':
import sys
n = eval(sys.argv[1])
m = eval(sys.argv[2])
print(solve(n, m))
import pytest
from problem import solve, markov_matrix
_test_markov_matrix = (
(2, 0, [[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]]),
(2, 1, [[0, 1/2, 1/2, 0],
[1/2, 0, 0, 1/2],
[1/2, 0, 0, 1/2],
[0, 1/2, 1/2, 0]]),
(2, 2, [[1/2, 0, 0, 1/2],
[0, 1/2, 1/2, 0],
[0, 1/2, 1/2, 0],
[1/2, 0, 0, 1/2]]),
(2, 3, [[0, 1/2, 1/2, 0],
[1/2, 0, 0, 1/2],
[1/2, 0, 0, 1/2],
[0, 1/2, 1/2, 0]]),
)
@pytest.mark.parametrize('n,m,expect', _test_markov_matrix)
def test_markov_matrix(n, m, expect):
assert expect == markov_matrix(n, m).tolist()
_test_solve = (
(2, 0, 0),
(2, 1, 1),
(2, 2, 1),
(2, 3, 1),
(2, 4, 1),
(2, 5, 1),
(3, 0, 0),
(3, 1, 2.725309),
(4, 0, 0),
(4, 1, 4.777778),
)
@pytest.mark.parametrize('n,m,expect', _test_solve)
def test_solve(n, m, expect):
assert expect == solve(n, m)