Correct numerical modeling of the 3D Dirac Delta function

• Telemachus
In summary, the conversation discusses the testing of a code for the diffusion equation using the analytical solution for infinite media with a Dirac delta source term. The code is not giving the analytical solution, and there are several possible reasons for this. One reason could be that there is a bug in the code, while another reason could be that the modelling of the analytic situation is incorrect. The speaker is solving in a cube and looking at the solution at the center of the domain to avoid border effects. They also mention the equation they are solving and the parameters they have used. The speaker is concerned about the sensitivity to the values of sigma they are using to model the Dirac delta function and how it may be affecting the solution. They also mention a
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:

Fortran:
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.

Last edited:
Telemachus said:
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?

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.

Telemachus said:
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?

Telemachus said:
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
Orodruin said:
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.

Last edited:
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!

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!

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.

Attachments

• num1000.png
9.5 KB · Views: 627
• num1000-png.png
9.5 KB · Views: 655
Last edited:

1. What is the 3D Dirac Delta function?

The 3D Dirac Delta function, also known as the three-dimensional Dirac delta distribution, is a mathematical concept used in physics and engineering to represent a point source or a localized concentration of mass or charge in three-dimensional space. It is characterized by having an infinite value at the origin and zero value everywhere else, and it integrates to one over all space.

2. How is the 3D Dirac Delta function used in numerical modeling?

In numerical modeling, the 3D Dirac Delta function is used to simulate point sources or localized concentrations in a three-dimensional space. This allows for more accurate and efficient calculations in simulations, as it avoids the need for discretization of the entire space. It is especially useful in solving partial differential equations, such as the wave equation, in which point sources are common.

3. What are some challenges in correctly modeling the 3D Dirac Delta function?

One challenge in correctly modeling the 3D Dirac Delta function is its infinite value at the origin, which can cause numerical instability and errors in simulations. Another challenge is accurately representing its behavior in all three dimensions, as it is a three-dimensional concept.

4. How can one ensure correct numerical modeling of the 3D Dirac Delta function?

To ensure correct numerical modeling of the 3D Dirac Delta function, one can use regularization techniques, such as replacing the delta function with a smooth function with similar properties. It is also important to carefully choose the numerical scheme and parameters to minimize errors and instabilities.

5. In what fields is the correct numerical modeling of the 3D Dirac Delta function particularly important?

The correct numerical modeling of the 3D Dirac Delta function is particularly important in fields such as fluid mechanics, electromagnetics, and quantum mechanics, where point sources or localized concentrations play a significant role in the behavior of the system. It is also crucial in simulations of physical phenomena, such as explosions, shock waves, and quantum tunneling.

• Set Theory, Logic, Probability, Statistics
Replies
12
Views
3K
• Calculus and Beyond Homework Help
Replies
3
Views
1K
• Calculus
Replies
2
Views
1K
Replies
5
Views
2K
• General Math
Replies
2
Views
2K
• Calculus and Beyond Homework Help
Replies
31
Views
2K
• Programming and Computer Science
Replies
26
Views
4K
Replies
3
Views
2K
• Calculus
Replies
2
Views
881
• Programming and Computer Science
Replies
1
Views
1K