Diffusion with Laplacian Operator (NumPy)

diffusion.py

import numpy as np
import matplotlib.pyplot as plt

def laplacian(grid):
    return np.roll(grid, 1, 0) +\
           np.roll(grid, -1, 0) +\
           np.roll(grid, 1, 1) +\
           np.roll(grid, -1, 1) - 4 * grid

def evolve(grid, df, D=1):
    return grid + df * D * laplacian(grid)

grid_shape = (10, 10)
grid = np.zeros(grid_shape)
block_low = int(grid_shape[0] * 0.4)
block_high = int(grid_shape[1] * 0.5)
grid[block_low:block_high, block_low:block_high] = 0.005

for i in range(100):
    grid = evolve(grid, 0.1)

plt.imshow(grid)
plt.show()