In mathematics, the walk-on-spheres method (WoS) is a numerical probabilistic algorithm, or Monte-Carlo method, used mainly in order to approximate the solutions of some specific boundary value problem for partial differential equations (PDEs).[1] [2] The WoS method was first introduced by Mervin E. Muller in 1956 to solve Laplace's equation,[1] and was since then generalized to other problems.
It relies on probabilistic interpretations of PDEs, and simulates paths of Brownian motion (or for some more general variants, diffusion processes), by sampling only the exit-points out of successive spheres, rather than simulating in detail the path of the process. This often makes it less costly than "grid-based" algorithms, and it is today one of the most widely used "grid-free" algorithms for generating Brownian paths.
Let
\Omega
Rd
\Gamma
\Gamma
x
\Omega
Consider the Dirichlet problem:
\begin{cases}\Deltau(x)=0&ifx\in\Omega\\ u(x)=h(x)&ifx\in\Gamma.\end{cases}
It can be easily shown that when the solution
u
x\in\Omega
u(x)=Ex[h(W\tau)]
where is a -dimensional Wiener process, the expected value is taken conditionally on, and is the first-exit time out of .
To compute a solution using this formula, we only have to simulate the first exit point of independent Brownian paths since with the law of large numbers:
Ex[h(W\tau)]\sim
| 1 | |
| n |
| n | |
| \sum | |
| i=1 |
| i | |
| h(W | |
| \tau) |
The WoS method provides an efficient way of sampling the first exit point of a Brownian motion from the domain, by remarking that for any -sphere centred on, the first point of exit of out of the sphere has a uniform distribution over its surface. Thus, it starts with equal to, and draws the largest sphere
l{S}0
l{S}0
According to intuition, the process will converge to the first exit point of the domain. However, this algorithm takes almost surely an infinite number of steps to end. For computational implementation, the process is usually stopped when it gets sufficiently close to the border, and returns the projection of the process on the border. This procedure is sometimes called introducing an
\varepsilon
\varepsilon
Choose
\varepsilon>0
x(0)=x
d(x(n),\Gamma)>\varepsilon
rn=d(x(n),\Gamma)
\gamman
x(n+1):=x(n)+rn\gamman
d(x(n),\Gamma)\le\varepsilon
xf:=p\Gamma(x(n))
x(n)
\Gamma
xf
The result
xf
\Omega
x
In some cases the distance to the border might be difficult to compute, and it is then preferable to replace the radius of the sphere by a lower bound of this distance. One must then ensure that the radius of the spheres stays large enough so that the process reaches the border.[1]
As it is a Monte-Carlo method, the error of the estimator can be decomposed into the sum of a bias, and a statistical error. The statistical error is reduced by increasing the number of paths sampled, or by using variance reduction methods.
The WoS theoretically provides exact (or unbiased) simulations of the paths of the Brownian motion. However, as it is formulated here, the
\varepsilon
l{O}(\varepsilon)
When it is possible to use it, the Green's function first-passage (GFFP) method is usually preferred, as it is both faster and more accurate than the classical WoS.[3]
It can be shown that the number of steps taken for the WoS process to reach the
\varepsilon
l{O}(|log(\varepsilon)|)
As it is commonly the case for Monte-Carlo methods, this algorithm performs particularly well when the dimension is higher than
3
This method has been largely studied, generalized and improved. For example, it is now extensively used for the computation of physical properties of materials (such as capacitance, electrostatic internal energy of molecules, etc.). Some notable extensions include:
The WoS method can be modified to solve more general problems. In particular, the method has been generalized to solve Dirichlet problems for equations of the form
\Deltau=cu+f
More efficient ways of solving the linearized Poisson–Boltzmann equation have also been developed, relying on Feynman−Kac representations of the solutions.[7]
Again, within a regular enough border, it possible to use the WoS method to solve the following problem :
\begin{cases}\partialtu(x,t)+
| 1 | |
| 2 |
\Deltaxu(x,t)=0&ifx\in\Omegaandt<T\\ u(x,T)=h(x,T)&ifx\in\bar{\Omega}\\ u(x,t)=h(x,t)&ifx\in\Gamma.\end{cases}
of which the solution can be represented as:[8]
u(x,t)=Et,x(h(XT,T\wedge\tau))
where the expectation is taken conditionally on
Xt=x
To use the WoS through this formula, one needs to sample the exit-time from each sphere drawn, which is an independent variable
\tau0
R
E(\exp(-s\tau0))=
| R\sqrt{2s | |
The total time of exit from the domain
\tau
T
The WoS method has been generalized to estimate the solution to elliptic partial differential equations everywhere in a domain, rather than at a single point.[10]
The WoS method has also been generalized in order to compute hitting times for processes other than Brownian motions. For example, hitting times of Bessel processes can be computed via an algorithm called "Walk on moving spheres".[11] This problem has applications in mathematical finance.
The WoS can be adapted to solve the Poisson and Poisson–Boltzmann equation with flux conditions on the boundary.[12]
Finally, WoS can be used to solve problems where coefficients vary continuously in space, via connections with the volume rendering equation.[13]