Relation of Brownian Motion to Helmholtz Equation

one can obtain solutions to the Laplace equation $$\Delta\psi(x) = 0$$

or even for the Poisson equation $\Delta\psi(x)=\varphi(x)$ in a Dirichlet boundary value problem using a random-walk approach, see e.g. Introduction to Brownian Motion.

Now, already this fascinating concept is not really clear to me. Nevertheless, it would be worth to get into details if there also exists such a connection to the Helmholtz equation $$\Delta\psi(x) +k^2\psi(x)= 0$$

Hence my question:

Can we use some random walk to calculate solutions of the Helmholtz equation numerically?

It would also be interesting if this is still true for open systems where the boundary conditions are different than in the Dirichlet case and for which $k$ is now a domainwise constant function.
Thank you in advance

Robert

Solution 1:

The general form for the infinitesimal generator of a continuous diffusion in $\mathbb{R}^n$ is $$ Af(x) = \frac12\sum_{ij}a_{ij}\frac{\partial^2 f(x)}{\partial x_i\partial x_j}+\sum_ib_i\frac{\partial f(x)}{\partial x_i}-cf(x).\qquad{\rm(1)} $$ Here, $a_{ij}$ is a positive-definite and symmetric nxn matrix, $b_i$ is a vector and $c$ is a non-negative scalar, with the coefficients $a,b,c$ functions of position $x$. Such operators are said to be semi-elliptic second order differential operators. The case with $c=0$ is the most common - being the generator of a Markov process (or semigroup). However, the $c > 0$ case does occur, and is then a generator of a submarkovian semigroup.

The coefficients $a,b,c$ can be understood as follows: $a$ gives the covariance matrix of the motion over small time intervals (i.e., the level of noise). $b$ gives the mean over small time intervals (the drift) and $c$ is the rate at which the process is "killed". To be precise, a processes $X_t$ can be modeled formally by adding an additional state $\it\Delta$ called the cemetary. So, we represent the state space for the killed diffusion as $\mathbb{R}^n\cup\{{\it\Delta}\}$. In a small time interval $\delta t$, the process has probability $c\delta t$ of being killed, in which case it jumps to the cemetary, and stays there. So, $X_t={\it\Delta}$ for all $t\ge\tau$ with $\tau$ being the (random) time when the process is killed. The terminology I am using here is taken from Revuz & Yor (Continuous Martingales and Brownian Motion), and will vary between different authors.

Anyway, getting back to PDEs. Suppose we want to solve the PDE $A\psi(x)=0$ on an open domain $U\subseteq\mathbb{R}^n$ with boundary condition $\psi(x)=\psi_0(x)$ for $x$ on the boundary of $U$ ($\partial U$, say). You can do the following. Simulate the process $X_t$ with initial condition $X_0=x$. Wait until the first time $T$ at which it hits the boundary and, when this occurs (if the process doesn't get killed first. i.e., $T < \tau$), take the expected value. $$ \psi(x)=\mathbb{E}_x\left[\psi_0(X_T)1_{\{T < \tau\}}\right].\qquad\qquad{\rm(2)} $$ Then $\psi$ satisfies the PDE $A\psi=0$.

This is all very general. Getting back to the Helmholtz equation, we can let $a_{ij}$ be the identity matrix and $b_i=0$, and $c$ is a constant. In that case our generator becomes $A\psi=\frac12\Delta\psi-c\psi$. [Edit: This is not quite the same as the Helmholtz equation, which has $c=-k^2$, because here we have $c > 0$. There is a sign difference which changes the behaviour of the solutions. See below] The process then is the following: run a Brownian motion starting from $x$ until the first time $T$ it hits the boundary. Decide if it has not been killed, which has probability $e^{-cT}$ conditional on $T$. If it hasn't, take the value $\psi_0(X_T)$. Finally, take the average of this process (e.g., using Monte Carlo). There is one practical issue here though. Throwing away all the paths on which the process gets killed is a bit wasteful, so you would simply multiply by the probability of not being killed on each path, rather than actually discarding them. i.e., you simulate a regular Brownian motion, and then calculate $$ \psi(x)=\mathbb{E}_x\left[\psi_0(X_T)e^{-cT}\right].\qquad\qquad{\rm(3)} $$

We can even go the whole hog and solve $A\psi(x)=\varphi(x)$ for a general $x$-dependent coefficients and source term $\varphi$, $$ \begin{align} \psi(x)&=\mathbb{E}_x\left[1_{\{T<\tau\}}\psi_0(X_T)-\int_0^{T\wedge\tau}\varphi(X_t)\,dt\right]\\ &=\mathbb{E}_x\left[e^{-\int_0^Tc(\hat X_s)\,ds}\psi(\hat X_T)-\int_0^Te^{-\int_0^tc(\hat X_s)\,ds}\varphi(\hat X_t)\,dt\right]. \end{align}\qquad{\rm(4)} $$ Here, $X$ is the process killed at (state-dependent) rate $c$, and I'm using $\hat X$ for the process without killing, which requires multiplying by the survival probabilities $e^{-\int c(\hat X_s)\,ds}$ instead.

One other area in which you have a '$-cf$' term in the PDE governing diffusions is in finance, and it occurs in two different, but closely related ways. Prices of financial assets are frequently modeled as diffusions (even, jump-diffusions), and the value of a financial derivative would be expressed as the expected value of its future value - under a so-called "risk-neutral measure" or "martingale measure" (which are just a special probability measures). However, you need to take interest rates into account. If the rate is $r$, then you would multiply the future (time $t$) value by $e^{-rt}$ before taking the expected value, which is effectively the same as adding a $-rf(x)$ term to the generator. And, as in the general case above, $r$ can be a function of the market state. The second main way (which occurs to me) in which such terms appear in finance is due to credit risk. If a counterparty has probability $rdt$ of defaulting in any time $t$, then you would have a $-rf(x)$ term occurring in the diffusion. This is more in line with the "killing" idea discussed above, but behaves in much the same way as interest rates.

Finally, I'll mention that the PDE in the more general time-dependent situation is of the form $\partial f/\partial t + Af =0$, where $f$ is the expected value of some function of the process at a future time (and not necessarily the first time it hits the boundary). As mentioned in the other answer this is sometimes known as the Feynman-Kac formula, generally by physicists, and also as the Kolmogorov backward equation by mathematicians. Actually, the backward equation in the Wikipedia link doesn't have the $-cf$ term, but it would in the more general case of diffusions with killing. The adjoint PDE applying to probability densities is known as the Fokker-Planck equation by physicists and the Kolmogorov forward equation to mathematicians.

Edit: As mentioned above, what we have here does not quite correspond to the Helmholtz equation, because of the sign of $c$, and the behaviour of the solutions does change depending on whether $c$ is positive or negative. In the probabilistic interpretation, $c > 0$ is the rate at which the process is killed. Looking at (3) and (4), we can see that solutions to the PDE will decay exponentially as we move further from the boundary. Furthermore, if the values on the boundary are non-negative, then $\psi$ has to be non-negative everywhere. The probabilistic method naturally leads to $\psi(x)$ being a positive linear combination of its boundary values (i.e., an integral with respect to a measure on the boundary). On the other hand, the Helmholtz equation has oscillating wavelike solutions. The values of $\psi(x)$ can exceed its values on the boundary and, even if $\psi\vert_{\partial U}$ is positive, it is possible for $\psi(x)$ to go negative inside the domain. So, it is not a positive linear combination of its boundary values. We could just try using a negative $c$ in (3) and (4) but, for the reasons just mentioned, this cannot work in general. What happens is that $e^{\vert c\vert T}$ is not integrable. To get around this, it is possible to transform the Helmholtz equation so that the zeroth order coefficient $-c$ is positive. We can make a substitution such as $\psi(x)=\tilde\psi(x)e^{ikS(x)}$, where $S$ is any solution to $\Vert\nabla S\Vert=1$. Then, the Helmholtz equation becomes, $$ \frac12\Delta\tilde\psi + ik\nabla S\cdot\nabla\tilde\psi + \frac{ik}{2}(\Delta S)\tilde\psi=0. $$ So we have a zeroth order term of the form $-\tilde c\tilde\psi$ for $\tilde c=-ik\Delta S/2$. This is imaginary, so does not make sense as a "killing" rate any more. However, as its real component is nonnegative (zero here), equations such as (3) and (4) above give bounded and well-defined expressions for $\tilde\psi$. A google search gives the following paper which uses such a substitution: Novel solutions of the Helmholtz equation and their application to diffraction. I expect that the papers linked to in the other answer use similar techniques to transform the equation into a form which can be handled by the probabilistic method - although I have not had a chance to look at them yet, and are not free access.

Solution 2:

It does not seem that anybody has a more complete reply, so I might as well summarize what I have said in the comments. The technique you are refering to in the OP is an application of the Feynman-Kac formula.

This technique has in fact been applied to solving the Helmholtz equation, as a Google search will reveal:

• Probabilistic solutions of the Helmholtz equation

• Random walk approach to wave propagation in wedges and cones