Correct numerical modeling of the 3D Dirac Delta function

[Telemachus]

Hi. I was trying to test a code for the diffusion equation, using the analytical solution for infinite media, with a Dirac delta source term:

$q(\mathbf{r},t)=\delta (\mathbf{r}) \delta (t)$.

The code is not giving the analytical solution, and there might be several reasons why this is so. One reason could be that my code has something wrong, a bug, something I did wrong. Another reason might be that my modelling of the analytic situation is incorrect, for example, the boundary conditions might be affecting my solution (my numerical model is of course not infinite).

I am solving in a cube, and looking at the solution at the center of the domain in order to avoid border effects.

The equation I am solving is: $\displaystyle \frac{1}{c}\frac{\partial \phi(\mathbf{r},t)}{\partial t}-D\nabla^2 \phi(\mathbf{r},t)=q(\mathbf{r},t)$.

I have tried with a cube of side two, and of side one. And I plot the solution at a distance of 0.2 units from the source.

So, to model the Dirac delta source I was using:

      sigg=0.1d0
      siggan=0.01d0
      Power=20.d0
      rxs=0.5d0 ! x position of source
      rys=0.d0 ! y position of source
      rzs=0.5d0 ! z position of source
      rts=0.3d0 ! at this time the pulse is emitted

        do k=1,zrows
         do j=1,yrows
          do i=1,xrows
           rq(i,j,k)= Power*dexp(-((x-rxs)**2 +
     +       (y-rys)**2+(z-rzs)**2)/sigg**2)*dexp(-(t-rts)**2/siggan**2)
          enddo
         enddo
        enddo

There is also, of course, a loop in time.

However, I have noted that this is very sensitive to the sigma values I use. I am also using a few points in the spatial coordinates, but I don't think that is an issue (perhaps it could be affecting the distance I am considering to the source, but not more than that, the solution looks converged, and if I use many more points I don't get considerable changes).

I am a bit worried about this, because I thought the code was working fine, the solutions looks physically as I was specting, and I get curves similar to the one reported in the papers, but when I've tried to see in detail how it behaves with respect to an analytic solution, I've found with this issue.

So first, I would like to know if my modelling of the Dirac delta function is appropriate. Second, how big should be my domain in order to avoid any effects from the boundaries? I think that, due to the small distance to the source, and considering the speed of propagation of the wave, the borders shouldn't be bothering. Should I look my problem somewhere else in the code? or should I focus in the proper modelling of the analytic situation, with an infinite homogeneous media, and the Dirac delta source term?

A possible way to model a one dimensional Dirac delta functions is to use that:

$\delta (x)=\displaystyle \lim_{\sigma \to 0} \frac{e^{-x^2/\sigma^2}}{\sigma\sqrt{\pi}}$

For a 3D dirac delta, with an impulse in time, I thought of using:

$\delta (\mathbf{r})\delta (t)=\displaystyle \lim_{\sigma \to 0} \frac{e^{-[(x-x̣_0)^2+(y-y_0)^2+(z-z_0)^2+(t-t_0)^2]/\sigma^2}}{\sigma^4 \pi^2}$

However, dividing for a small sigma seemed to cause trouble. I am not totally sure about this. But I think that the modelling of the Dirac delta function causes some problem, and that might be the reason why I am not being able to get the analytical situation. When I do this numerically, I take a small $\sigma$, and I've noted that small variations on it produces big changes in the solution.

 

[Orodruin]

A possible way to model a one dimensional Dirac delta functions is to use that:

$\delta (x)=\displaystyle \lim_{\sigma \to 0} \frac{e^{-x^2/\sigma^2}}{\sigma\sqrt{\pi}}$

For a 3D dirac delta, with an impulse in time, I thought of using:

$\delta (\mathbf{r})\delta (t)=\displaystyle \lim_{\sigma \to 0} \frac{e^{-[(x-x̣_0)^2+(y-y_0)^2+(z-z_0)^2+(t-t_0)^2]/\sigma^2}}{\sigma^4 \pi^2}$

However, dividing for a small sigma seemed to cause trouble. I am not totally sure about this. But I think that the modelling of the Dirac delta function causes some problem, and that might be the reason why I am not being able to get the analytical situation. When I do this numerically, I take a small $\sigma$, and I've noted that small variations on it produces big changes in the solution.

You cannot use these definitions when you have discretised your spatial dimensions. You will end up with one point where the source term is $\sigma$-dependent and for all other grid points you will get essentially zero. You must adapt your approximation of the delta function to your grid. Essentially, what you need to do is to put a source term in a single grid point with the appropriate strength, that is the best approximation of the delta function you can hope to do on a grid. You can never do better than the resolution provided by the grid.

 

[Telemachus]

Ok. So my delta function should be something like:

$\delta(\mathbf{r})\delta (t) = \left\{ \begin{array}{@{}l@{\thinspace}l} A & \text{if} & \mathbf{r}=\mathbf{r}_0 & \text{and} & t=0 \\ 0 & \text{in any other case} & \\ \end{array} \right.$

Is that right? now, the impulse magnitude is 1 for the analytical situation. How do I determine the magnitude of A?

 

[Orodruin]

That would depend on your step size. What is the defining characteristic of the delta function? Compare to the corresponding statement in the discretised case.

 

[Telemachus]

I think that the defining characteristic is that $\int_{-\infty}^{+\infty}\delta(x)dx=1$, and that $\delta(0)=\infty$. How can I relate this to the step size?

However, I've been trying with the piece wise definition for the delta function, and the solution gives zero everywhere. I have tried with different values for A.

 

[Orodruin]

I think that the defining characteristic is that $\int_{-\infty}^{+\infty}\delta(x)dx=1$, and that $\delta(0)=\infty$. How can I relate this to the step size?

Almost, but not quite. The delta function is not really a function, it is a distribution, it makes no sense to talk about its value at zero. The defining feature is that
$$
\int \delta(x) f(x) dx = f(0)
$$
as long as $f(x)$ is a sufficiently nice function. What would be the discretised version of this?

 

[Telemachus]

Ok, right.

The discretised version would be, I think: $\sum_i^N \Delta x f(x_i) \delta(x_i)=f(0)$, where $\Delta x$ is the grid spacing, and the sum is over every point in the grid. Is that right?

So, from this, I think I can define: $\delta(x_i)=0$ for all $x_i\neq 0$ and $\delta (x_i)=1$ for $x_i=0$.

So, maybe under this setting I could think of defining A as $f(0)/\Delta x$, but the thing is that in my problem I don't have the test function $f(x)$ anywhere. Am I right?

Thanks for your guidance.

 

[Orodruin]

So, from this, I think I can define: $\delta(x_i)=0$ for all $x_i\neq 0$ and $\delta (x_i)=1$ for $x_i=0$.

Almost. What would the sum give if you put $\delta(x_i) = 1$ for $x_i = 0$?

So, maybe under this setting I could think of defining A as $f(0)/\Delta x$, but the thing is that in my problem I don't have the test function $f(x)$ anywhere. Am I right?

It does not really matter whether you actually have a test function or not, the point is in how you define $\delta$. The $f(0)$ should come from the $f(x)$ in the sum/integral.

 

[Telemachus]

Almost. What would the sum give if you put $\delta(x_i) = 1$ for $x_i = 0$?

I would have $f(0) \delta(0) \Delta x=f(0)$, I didn't realize of this before, does this mean that $A=1/\Delta x$?

Now, the issue I have using this definition, is that the program gives zero. I am defining the source at a single point in space and time, and giving it the magnitude $A=\frac{1}{\Delta x \Delta y \Delta z \Delta t}$. So basically (in the sake of clarity) my discrete Dirac delta function is (with $x_i=(i-1)\Delta x,y_j=(i-j)\Delta y, z_k=(k-1)\Delta z, t_m=(m-1)\Delta t$:

$\delta_{i,j,k,m}= \left\{ \begin{array}{} \frac{1}{\Delta x \Delta y \Delta z \Delta t} & \text{if} & i,j,k,m=1. \\ 0 & \text{in other case} & \\ \end{array} \right.$

Edit: it doesn't give zero, I had an error in the program. However, I still don't have the analytic solution. It looks like if there is a difference of a factor. I will check the program again and see if find more mistakes.

 

[Orodruin]

Not sure if you really want to put the delta at $x_1$, $y_1$, $z_1$. Make sure you have put it in the correct position. Otherwise, yes, it is the closest discretised approximation you can make to the delta distribution.

 

[Telemachus]

Yes, it was only illustrative, I have actually put it at the middle of the domain to avoid boundary effects. However, I have tried by considering a bigger domain, and the solution got closer to the analytical, so there are clearly other things I have to look at.

Thank you very much!

 

[Orodruin]

You have a typical time-scale in your problem which is $t_0 = L^2/Dc$, where $L$ is the distance to your boundary. For times $t \ll t_0$, your solution should be reminiscent of the solution in all of space. Once $t \sim t_0$ you will invariably suffer from boundary effects.

 

[Telemachus]

Thanks a lot!

 

[Telemachus]

Ok, actually the equation I am solving is this one (I've introduced the thread with diffusion because it was easier, and it is asymptotic to it for large $\mu_s$:

$\frac{1}{c}\frac{\partial I(\mathbf{r},\hat \Omega, t)}{\partial t}+ \hat \Omega \cdot \nabla I(\mathbf{r},\hat \Omega, t)+\mu_s I(\mathbf{r},\hat \Omega, t)=\frac{\mu_s}{4\pi}\int_{4\pi} I(\mathbf{r},\hat \Omega, t) d\hat \Omega+\delta(\mathbf{r})\delta(t)$

This is a linearized Boltzmann equation.

What I am plotting is the photon scalar flux, which is defined as: $\phi(\mathbf{r},t)=\frac{1}{4\pi}\int_{4\pi} I(\mathbf{r},\hat \Omega, t) d\hat \Omega$.

The analytic solution for infinite (3D) homogeneous media also involves a Dirac delta function, which physically represents the unscattered particles:

$\phi(\mathbf{r},t)_{an}=\frac{e^{-ct\mu_s}}{4\pi r^2}\delta(r-ct)+\frac{(1-r^2/c^2t^2)}{(4\pi ct/3\mu_s)} e^{-ct\mu_s}\times e^{ct\mu_s[1-r^2/c^2 t^2]}\sqrt{1+\frac{2.026}{ct\mu_s[1-r^2/c^2 t^2]}} \Theta (ct-r)$

Here $\Theta$ is the heaviside step function.

This is actually an approximation to the analytic solution, which is not exactly known in closed form, but it is proven in a paper that the error is of the order of 2% to the truth solution.

The thing is that this solution involves a Dirac delta function, so I was thinking that perhaps I am not being able to numerically reproduce this behavior because it is not so well behaved. The domain I am working with has 10 units per side, I am using c=1, and in the picture I have used $\mu_s=1$.

I have also tried with $\mu_s=0.1$, and the propagating wave looks like a gaussian, it doesn't resemble the Dirac delta that is clearly seen in the analytic solution.

The r for the analytic solution is $r=\sqrt{0.12}$.

Also, how should I work the dimensional analysis for this equation in order to obtain a relation for $t_0$?

I am thinking that perhaps the smoothing of the Dirac delta impulse is due to the spatial grid. The analytical solution is obtained from Fourier space, and then a Fourier transformation is used. It is exact in that sense, the approximation comes from the fact that the inverse has no closed solution in 3D, but it does in 2D and 4D, so an interpolation formula is used.

So, I was thinking that perhaps I should look far away from the source, so the delta function, within a scattering medium, will be hardly attenuated far away from the source, and I would expect that the numerical situation looks more like the analytical one. The stepsize I am using is $\Delta x=\Delta y=\Delta z=0.1$. I was trying to avoid using a huge domain, because the computational work involved scales very rapidly.

The numerical solution I've obtained in this regime (the one in the picture) is probably wrong due to another fact: I should have used a much finer discretization in the $\hat Omega$ than the one I was using for the highly scattering regime.