Harmonic function Monte Carlo

Approximation of a discrete harmonic function by Monte Carlo method

An harmonic function \(f\) on a domain of \(\mathbb{Z}^2\) satisfies, for each couple of integer numbers \((x,y)\) in the domain, the relation \[ f(x,y) = ( f(x+1,y) + f(x-1,y) + f(x,y+1) + f(x,y-1) )/4 . \] Therefore, in all points the function \(f\) is equal to its arithmetic mean of its 4 adiacent points.

The values token by \(f\) inside the domain are identified by the values token on the boundary. Let \((X_t,Y_t)_{t \in \mathbb{N}}\) be the symmetric random walk on \(\mathbb{Z}^2\) starting from \((X_0,Y_0)=(x,y)\) and let \(T\) be the time needed to reach the boundary, then \[f(x,y) = \mathbb{E}_{(x,y)}f(X_T,Y_T)\] the expected value of \(f\) when the random walk reaches for the first time the boundary of the domain.


Here the domain is the square of integers \([0,K]^2\) and its interior is the square \(I^2\) with \(I = [1,K-1]\) and \(K=25\).

The values token on the boundary are

  • \(f(0,y) = f(K,y) = K\) for each \(y \in I\) and
  • \(f(x,0) = f(x,K) = 0\) for each \(x \in I\).

For each iteration \(N\), from each point \((x,y) \in I^2\), we start an independent random walk, and simulate each random walk until it reaches the boundary. Let \(F_N(x,y)= f(X_T,Y_T)\) be the value token by \(f\) at that hitting time. Due to the Strong Law of Large Numbers, the arithmetic mean \[ \bar F_N(x,y) = (F_1(x,y) + \dots + F_N(x,y)) / N \] converges almonst surely to \(\mathbb{E}_{(x,y)}f(X_T,Y_T)\) and it is, then, an approximation of the harmonic function \(f(x,y)\).


This is an evolutive simulation, i.e. one click on the button [COMPUTE] executes \(N\) additional iterations to the previous execution, with \(N \in \left\{1, 10, 100\right\}\). In order to reset the simulation and set the number of simulations back to 0, set \(N\) to 0.