1D Heat Equation with Insulated Boundaries with Dirac Delta

Question

This question is very similar to the question I previously asked here; however, there is a subtlety which I did not appreciate at the time. The previous question which I asked was correctly answered, and I accepted the answer, and so I am asking a new question instead of revising the previous one.

I would like to determine the solution to the 1D heat equation where the initial condition is a Delta function arbitrarily close to the boundary. Previously, I asked for the solution where the Delta function was at the boundary;

• $$ \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2}$$
• $$x \in [0, L], t \in [0, \infty]$$
• $$u(0,x) = \beta \delta(x)$$
• $$ u_{x} (0,t) = u_{x}(L,t) =0$$

Now I would like to know the solution for: $$u(0,x) = \beta \delta(x-\epsilon)$$ And then consider the limit $\epsilon \to 0.$

Physically, a finite amount of mass $\beta$ is concentrated at one end of an insulated bar of length $L$ at time $t=0$ and I want to know the evolution of this mass distribution throughout time. The total amount of mass in the bar should always remain $\beta,$ and the steady solution should be $u(x,\infty) = \frac{\beta}{L}.$

It turns out the limit $\epsilon \to 0$ is not the same as the solution $\epsilon=0,$ because roughly speaking only half of the mass at the boundary will diffuse into the region. This has a practical effect in that a numerical simulation which starts out with all of the mass very close to the boundary will disagree with the analytic solution.

Attempt at a solution: I believe that using the method of images should suffice here. That is, to get the boundary conditions, simply add two fundamental solutions: $$\beta u_0(x-\epsilon,t) + \beta u_0(x+\epsilon, t)$$ and then sum these centered over the boundaries at $2L.$ It would be nice to get confirmation that this is correct.

Answer 1

Using the method of separation of variables, one obtains the following general solution to the PDE satisfying the given boundary conditions: $$ u(x,t)=a_0+\sum_{n=1}^{\infty}a_n e^{-\frac{n^2\pi^2Dt}{L^2}}\cos\left(\frac{n\pi x}{L}\right). \tag{1} $$ The coefficients $a_n$ are determined by the initial condition $u(x,0)=\beta\delta(x-\epsilon)$, and are given by (assuming $0<\epsilon<L$) $$ a_0=\frac{1}{L}\int_0^L\beta\delta(x-\epsilon)\,dx=\frac{\beta}{L}, \tag{2a} $$ $$ a_n=\frac{2}{L}\int_0^L\beta\delta(x-\epsilon)\cos\left(\frac{n\pi x}{L}\right)dx= \frac{2\beta}{L}\cos\left(\frac{n\pi \epsilon}{L}\right)\qquad(n=1,2,\ldots). \tag{2b} $$ Plugging $(2)$ into $(1)$, and taking the limit $\epsilon\to 0^+$, one finally obtains $$ u(x,t)=\frac{\beta}{L}\left[1+2\sum_{n=1}^{\infty}e^{-\frac{n^2\pi^2Dt}{L^2}}\cos\left(\frac{n\pi x}{L}\right)\right]. \tag{3} $$

Answer 2

Your idea is pretty much right. You have to reflect the singularities in both boundaries, so you will have delta functions at $2nL \pm \epsilon$, for $n\in \mathbb Z$: $$ u(x,t) = \frac{\beta}{2\sqrt{\pi D t}}\sum_{n=-\infty}^\infty \exp\left(-\frac{(x+2nL+\epsilon)^2}{4Dt}\right) + \exp\left(-\frac{(x+2nL-\epsilon)^2}{4Dt}\right). $$ As $\epsilon \to 0$, the two Gaussians will become identical. It approaches an initial condition with delta function of strength $2\beta$ on the boundary, hence will have the mass you expect.