Integrating over multivariate dirac deltafunctions

[Matthew Rocklin]

Hi Everyone, 

I'd like to compute multivariate integrals that contain Dirac Deltafunctions. I.e. expressions like

integrate(exp(-(x**2+y**2))/pi * delta(2*x+3*y), (x,-oo, oo), (y,-oo, oo))

The deltaintegrate function inside sympy fails to compute these correctly, see issue 2630.

Wikipedia says that you can compute general expressions of this form as follows: 

(hopefully the above image makes it through, if not go here)

http://upload.wikimedia.org/wikipedia/en/math/3/3/f/33fbbed28ec715257d268faefc9e0e9f.png

How hard would it be to compute the right hand side in sympy? In particular I'm confused by how to express the domain and what they mean by \delta\sigma(x)

Or, if there is a better way of going about this I'm happy to hear it. 

-Matt

[Tom Bachmann]

I think there are two basic cases to consider:

1) delta functions with *linear* arguments.

This is (I believe) by a large margin the most common case. For this
deltaintegrate should simply be fixed; one only needs to correctly
implement integral(f(x) delta(a*x+b), (x, c, d)) and this shouldn't be
hard. (Note that this will also yield the multi-dimensional cases, by
repeated integration.)

2) delta functions with more complicated arguments

Getting this to work is (I believe) much more work. It is non-trivial to
even define this.

Tom

On 20.12.2011 16:21, Matthew Rocklin wrote:
> Hi Everyone,
>
> I'd like to compute multivariate integrals that contain Dirac
> Deltafunctions. I.e. expressions like
>
> integrate(exp(-(x**2+y**2))/pi * delta(2*x+3*y), (x,-oo, oo), (y,-oo, oo))
>
> The deltaintegrate function inside sympy fails to compute these
> correctly, see issue 2630.

> <http://code.google.com/p/sympy/issues/detail?id=2630>
>
> Wikipedia says
> <http://en.wikipedia.org/wiki/Dirac_delta_functions#Properties_in_n_dimensions>

> that you can compute general expressions of this form as follows:
>

> \int_{\mathbb{R}^n} f(\mathbf{x}) \, \delta(g(\mathbf{x})) \,
> d\mathbf{x} =
> \int_{g^{-1}(0)}\frac{f(\mathbf{x})}{|\mathbf{\nabla}g|}\,d\sigma(\mathbf{x})
>
> (hopefully the above image makes it through, if not go here)
> http://upload.wikimedia.org/wikipedia/en/math/3/3/f/33fbbed28ec715257d268faefc9e0e9f.png
>
> How hard would it be to compute the right hand side in sympy? In
> particular I'm confused by how to express the domain and what they mean
> by \delta\sigma(x)
>
> Or, if there is a better way of going about this I'm happy to hear it.
>
> -Matt
>

> --
> You received this message because you are subscribed to the Google
> Groups "sympy" group.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to
> sympy+un...@googlegroups.com.
> For more options, visit this group at
> http://groups.google.com/group/sympy?hl=en.

[krastano...@gmail.com]

@Tom, the operation is well defined for all sufficiently nice functions.
@Matthew, The notation that is used there is a bit unclear (for me). The domain is all zeros of g and the sigma is just some measure theory notation meaning that you should actually take the sum.

To solve this in sympy you can try the following thing (it depends on how good 'solve' is in sympy):

0) to solve integrate(delta(f(x,y,...)*g(x,y,...)), (x,...), (y,...), ...)
1) take out the g(x,y,...) function from the delta function
2) find the zeros of g(...) (presumably using 'solve')
3) evaluate the derivatives of g(...) at those zeros
4) your integral is equal to the sum of f(...) / |grad(g(...))| evaluated at those zeros.

I think that a nice explanation can be found in Appel. I can send a pdf with the book to you if you are interested.

sympy+unsubscribe@googlegroups.com.

For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To post to this group, send email to sy...@googlegroups.com.

To unsubscribe from this group, send email to sympy+unsubscribe@googlegroups.com.

[Matthew Rocklin]

@Tom - I agree that it makes sense to separate out the common easy case from the very rare very difficult case. 

@Stephan - In the multivariate case the zeros are likely to consist of continuous lines/surfaces/manifolds of co-dimension 1. We'll need to do an integral over complex surfaces. This could be challenging. 

From my experience deltaintegrate has some issues, even in the linear case. The code also seems like it might be more complex than necessary, allowing odd bugs to filter through. Wolfram handles deltas more robustly than SymPy does so at least we know that there exists a nicer solution. I wonder if a system something like what Stephan outlined would be cleaner/more robust. 

To unsubscribe from this group, send email to sympy+un...@googlegroups.com.

[Aaron Meurer]

Does someone know of an example integral containing a delta function
with a non-linear argument, with the correct answer? That would help
with trying to find the algorithm.

Perhaps it can be done by a change of variable, to make the argument
linear. You probably have to be careful about some things, though,
since delta functions don't follow all the "rules" that normal
functions do.

Aaron Meurer

[krastano...@gmail.com]

I vaguely remember the trajectory of a relativistic particle presented in jackson's electrodynamics. I suppose it's convoluted enough for use in tests.

A simpler one will be some surface density over a strange surface. For example the total charge of a unitary surface charge distribution on a spherical surface (of radius = 1) will be:

<delta(sqrt(x**2+y**2+z**2)-1) | 1> = 4 pi

Well this was just the surface of the sphere so it was not that interesting, but one can calculate for example a field using such integral.

I'll try to dig some nice examples.

PS Mathematica can calculate those (at least in 1D and 2D)

In[6]:= Integrate[DiracDelta[-1 + Sqrt[x^2]], {x, -Infinity, Infinity}]
Out[6]= 2

In[7]:= Integrate[DiracDelta[-1 + Sqrt[x^2 + y^2]],  {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]
Out[7]= 2 Pi

[krastano...@gmail.com]

and also

In[8]:= Integrate[DiracDelta[-1 + Sqrt[x^2 + y^2 + z^2]],  {x, -Infinity, Infinity}, {y, -Infinity, Infinity}, {z, -Infinity, Infinity}]

Out[8]= 4 Pi

[Ondřej Čertík]

On Fri, Dec 30, 2011 at 2:39 AM, krastano...@gmail.com
<krastano...@gmail.com> wrote:
> I vaguely remember the trajectory of a relativistic particle presented in
> jackson's electrodynamics. I suppose it's convoluted enough for use in
> tests.

I have it at my table here, but I couldn't find it in the special
relativity and the subsequent chapter.

Ondrej