Symmetric Random Walk on a 2D grid

Question

I am doing a simulation to get the average temperature in a 2D grid heated only at the top at 100 degrees Celsius like below:

This is the equivalent of solving the diffusion equation:

\begin{cases} \Delta u(x,y) =0, \quad \forall x,y \in ]0,9[ \\ u(x,0)=u(0,y)=u(9,y)=0, \\ u(x,9)=100 \end{cases}

The probability to move from one node to another in the grid is p = 1/4

When I launch the process at a random position (x,y) in the above grid at each iterations over N = 10000 iterations, I get an average temperature of 25 as seen in the below code and image results. Why is it the case?

####Importing Libraries
import numpy as np  
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
rng = np.random.default_rng()

##Grid
Rmax = 9  ##rows
Cmax = 9   ##columns
T = np.zeros((Rmax,Cmax))  
T0 = 100
T[:,-1] = T0 # BOUNDARY CONDITIONS 

##2d Mesh
X1  = np.arange(0,Rmax,1) 
Y2 = np.arange(0,Cmax,1)
X,Y = np.meshgrid(X1,Y2)

Direction = np.array([(1,0),(-1,0),(0,1),(0,-1)]).reshape(4,2) ## 4 possible directions

N_iter = 10000 ##Iterations

##looping through all grid points except those at boundaries
for i in range(1,Rmax-1):
    for j in range(1,Cmax-1):
        T_boundary = []  ##store the temperature of terminated point
        for k in range(N_iter):  ##random walk
            x,y = np.random.randint(1,Rmax-1),np.random.randint(1,Cmax-1) ## random inital position
            while True:  ##running infinite loop
                R_n = np.random.randint(0,4)  ##selecting one of four directons
                next_direc = Direction[R_n,:]
                ##updating x and y coordinate
                x+= next_direc[0]
                y+= next_direc[1]

                if(x%(Rmax-1)==0 or y%(Cmax-1)==0 ):  ##terminating random walks at boundary 

                    T_boundary.append(T[x,y]) ##appending array
                    break;
        T[i,j] = np.mean(np.array(T_boundary)) ##channging grid temperature by an average

#plot
def f():
    z = T[X,Y]  ##giving value to each point of mesh 
    return z;

Z = f() 
# 3dplot
print(Z)

fig = plt.figure(figsize = (18,10))  
ax = plt.axes( projection='3d')
ax.plot_wireframe(X,Y,Z)
ax.set_title('Laplace Equation using Random Walk')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel(r'$T(x,y)$')
plt.show()

Answer 1

the simulation is not capturing the right problem. to get the temp. at a point $(x,y)$, start many random walks at that point, then average the results. Run a loop over the starting points $(x,y)$. There is no reason to initialize at a random point, you are just adding noise to the solution by doing that. As a sanity check, you can also solve the system of linear equations that govern the temperature, and compare to the simulation.

After doing that, one can average the temperatures obtained at the different points.

That the answer will be 25 is completely predictable. Call the desired average $A_1$. Define $A_2$ similarly, with boundary value 100 on the right and $0$ on other edges. Also, let $A_3$ correspond to the bottom edge only having boundary value 100, and $A_4$ to the left edge. Clearly $A_i$ are all equal, but also, $A_1+A_2+A_3+A_4=100$ because of linearity of the solution in the boundary conditions.

Answer 2

This is the sort of result you should have got

      0%      0%      0%      0%      0%      0%      0%    
0%   1.7%    3.2%    4.2%    4.5%    4.2%    3.2%    1.7%   0%
0%   3.8%    6.9%    8.9%    9.6%    8.9%    6.9%    3.8%   0%
0%   6.4%   11.7%   15.0%   16.2%   15.0%   11.7%    6.4%   0%
0%  10.3%   18.4%   23.3%   25.0%   23.3%   18.4%   10.3%   0%
0%  16.4%   28.2%   35.0%   37.1%   35.0%   28.2%   16.4%   0%
0%  26.9%   43.1%   51.2%   53.6%   51.2%   43.1%   26.9%   0%
0%  48.3%   66.1%   73.1%   74.9%   73.1%   66.1%   48.3%   0%
    100%    100%    100%    100%    100%    100%    100%    

You can read each values as the probability of a random walk from any particular point you eventually leave at the $100^\circ$ side. The probability at each point is going to be about the average of probabilities of the four points above, below, right and left.

Fortunately, because the boundaries are $100^\circ$ and $0^\circ$, you can also read the $\%$ as ${\,}^\circ$ to give a suggested temperature at each point.

The overall average is the probability of a random walk from a point chosen uniformly at random eventually leaves at the $100^\circ$ side. Given the symmetry of the square boundary, this is obviously $\frac14=25\%$.