Solution to diffusion equation of a random walk

Question

In my class of statistical physics, we studied the classic problem of random walk for the discrete case. In the end, we made the changes necessary for the master equation to be in the continuous format, in which the diffusion equation appeared:

$\frac{\partial}{\partial t} p(x,t) = D\frac{\partial^2 }{\partial x^2}p(x,t)$

So, in a homework assignment, the teacher asked to solve this equation. So I just used the standard method of separation of variables and, jumping to the end, I derived this result:

$p(x,t) = e^{-k^2t}\left(A\sin{\frac{K}{\sqrt{D}}x}+B\cos{\frac{K}{\sqrt{D}}x}\right)$

OK, assuming that the particle starts at the origin, $p(0,0) = 1$ so $B=1$. The next thing I did was to calculate $p(x,0)$:

$p(x,0) = \left(A\sin{\frac{K}{\sqrt{D}}x}+\cos{\frac{K}{\sqrt{D}}x}\right)$

And argued that the derivative at $x=0$ must be $0$, because this point should be a maximum since it's the only point that we know for sure what is the location of the particle, so all the other points should have $p(x,t) < 0 \ \forall x,t \neq0$. So I concluded that $A = 0$ and ended up with my 'final solution':

$p(x,t) = e^{-k^2t}\cos{\frac{K}{\sqrt{D}}x}$

Now to my question, this equation generates an inconsistency. First, the cosine ranges from $-1$ to $1$, in which $p$ can take negative values, which does not make sense since it measures a probability. Besides that, I have no clue how to normalize it, if I try to calculate

$\int_{0}^{\infty} e^{-k^2t} dt\int_{-\infty}^{\infty} \cos{\frac{K}{\sqrt{D}}x} dx$

The second integral gives me what? $0$? $\infty$?

I researched in the books I own about the diffusion equation, but all of them solve it in the context of a heat diffusion, so the problem has boundary condition at the ends of the material (a metal bar, a metal plate, etc.) that make appear infinite eigenfunctions and the final solution becomes an infinite sum of sines and cosines, rendering the heat function (that cannot take negative values too) always positive. But I don't have such boundary conditions in my random walk, it's just like a free particle. So it reminded me of the Quantum Mechanical formalism that says I have to square the wave function to calculate the probability.. it surely would make my solution always positive.

Answer 1

Separation implies summation

When we perform a separation of variables of this form,$$\rho(x, t)= X(x) T(t), \\ X\dot T = DX''T,\\ \dot T/T = D X''/X = -Dk^2 \text{ for some constant } k\\ \rho(x,t) = A \cos(kx)~e^{-Dk^2 t} + B \sin(kx)~e^{-Dk^2 t} \text{ for arbitrary } A,B,$$it is important to remember that we started from a really narrow starting place with $\rho$ and that this is not a fully general solution. The fully general solution requires, as you have said, taking a sum of these. In this case you can use the linearity of the differential equation $\dot \rho = D \rho''$ to argue that if you can determine that $$\rho(x,0) = \int_{-\infty}^{\infty}dk~\cos(kx)~f(k)$$ for some $f(k)$ then that is a sum of coefficients times the functions you got above, and the solution is then $\int dk~e^{-Dk^2t}~\cos(kx)~f(k).$

In other words: the boundary condition is the full function $\rho(x, 0)$. It's a 2D space consisting of half of $\mathbb R^2$ and your boundary condition fixes the function on that boundary $t=0.$

Fourier transforms to the rescue

When the boundary conditions are not either finite-range bounded or else periodic, the cosine-and-sine choice above proves somewhat awkward and we will find it easier to just complexify: $$ \rho(x, t) = \alpha e^{ikx}~e^{-D k^2 t} + \beta e^{-ikx} e^{-D k^2 t}.$$ We can then see that the two cases are really overlapping and once we take a sum over all $k$ we do not need both of them because we will get both positive and negative $k$, the separated solution can be simplified to just $\alpha e^{ikx} e^{-Dk^2 t}.$ Which is good because the equivalent integral there is known as the inverse Fourier transform, $$\rho(x, 0) = \frac{1}{2\pi}\int_{-\infty}^{\infty} e^{ikx}~dk~\int_{-\infty}^{\infty} e^{-ikx'}dx'~\rho(x', 0),$$ and therefore we have that by this prescription $$ \rho(x,t)=\frac{1}{2\pi}\int_{-\infty}^{\infty} e^{ikx} dk~e^{-Dk^2t}\int_{-\infty}^{\infty} e^{-ikx'} dx' \rho(x', 0).$$

But the boundary condition is a $\delta$-function

You have had some trouble specifying $\rho(x', 0)$ to be a probability density expressing certainty that a particle is at the origin. This is because to do that properly we need a function that is impossible (but surprisingly well-behaved!) called the Dirac $\delta$-function. It is a ‘function’ $\delta(x)$ such that $$\int_a^b dx~\delta(x) f(x) = \begin{cases}f(0)&\text{if } a < 0 < b,\\ -f(0)&\text{if } a > 0 > b,\\ 0&\text{otherwise}. \end{cases}$$ One typically thinks of it as an infinitely thin, infinitely large spike normalized so that the area under the curve is $1$, but other constructions also involve functions that oscillate wildly at $x\ne 0$ so that they cannot create nonzero integrals when multiplied by a standard smooth function. Part of this idea does mean that there is a great ambiguity about what happens if $a=0$ or $b=0$ in the integral above, too. But it can be realized as the limit of smooth functions e.g. a Gaussian with finite area 1 pushed to be thinner and thinner, so one can imagine that it itself is a smooth function too by just "backing off" from that limit by the tiniest bit.

Anyway in probability densities, the area under a probability density curve is the actual probability, and the Dirac $\delta$-function lets us define a probability density for "I know with probability $p$ that it is at point $x_0$" as the density contribution $p~\delta(x - x_0).$ In this case therefore you want to consider $\rho(x, 0) = \delta(x),$ and using the definition we immediately find that the Fourier transform of $\rho(x, 0)$ was just $e^{-ik0} = 1$, $$ \rho(x,t)=\int_{-\infty}^{\infty}\frac{dk}{2\pi}~ e^{ikx-Dk^2t}.$$

Complete the square, complete the derivation

So now we want to just complete the square, because we know the Gaussian integral that $\int_{-\infty}^\infty du~e^{-u^2} = \sqrt{\pi}.$ $$D~t~k^2 - ix~k =D~t~\left(k - \frac{ix}{2Dt}\right)^2 + \frac{x^2}{4Dt}.$$ Performing the Gaussian integral just gives a constant and the solution then comes out as $$\rho(x, t) = \frac{e^{-x^2/(4Dt)}}{\sqrt{4\pi Dt}}, $$a normal distribution with standard deviation $\sqrt{2Dt}.$

Now, since I have just shown you $\delta$-functions, I would like to stretch you one stage further, to tell you that actually linearity lets you use this result, instead of the separation-of-variables result, to see the evolution of an arbitrary system and alleviate your concerns about negative probabilities.

This solution that we came up with is special precisely because its boundary condition is a $\delta$-function. We call this result a “Green's function solution” of the equation for that reason. Every other function $\rho_0(x)$ could be expressed as $$\rho_0(x) = \int_{-\infty}^\infty dx'~\rho_0(x')~\delta(x - x'),$$ per the definition of $\delta$ above. As a consequence under the heat equation that function evolves over time $t$ to $$\rho(x, t) = \int_{-\infty}^{\infty} dx'~\rho_0(x')~\frac{e^{-(x-x')^2/(4Dt)}}{\sqrt{4\pi D t}},$$ from which it is easier to see that if $\rho_0(x)$ is nonnegative and integrates to $1$ then $\rho(x, t)$ for any fixed $t$ will also be nonnegative and integrate to $1$.