Re: Question about DiracDelta
- To: mathgroup at smc.vnet.net
- Subject: [mg68662] Re: Question about DiracDelta
- From: ab_def at prontomail.com
- Date: Mon, 14 Aug 2006 06:44:46 -0400 (EDT)
- References: <ebmth1$6g3$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
vasil michev wrote:
> Hi group, I need to plot this:
>
> Plot[NIntegrate[DiracDelta[3 - Cos[x] - Cos[y] - Cos[z] - e], {x, -Pi,
> Pi}, {y, -Pi, Pi}, {z, -Pi, Pi}, MaxPoints -> 1000000]/8/Pi^3, {e, 0,
> 6}, PlotPoints -> 30]
>
> the problem is that the function always returns 0, no matter what the
> value of 'e' is. I saw that other people have reported incorect
> evaluation/simplification related to DiracDelta, but I was wondering
> can anyone tell me if there's some workaround?
>
> Thanks in advance
NIntegrate doesn't attempt to process DiracDelta symbolically, it will
just evaluate it at the sampling points. The workaround is to use
PiecewiseIntegrate (
http://library.wolfram.com/infocenter/MathSource/5117/ ) and some
tricks. Integrate once to get rid of the delta function:
In[2]:= PiecewiseIntegrate[
1/(8*Pi^3)*DiracDelta[3 - Cos[x] - Cos[y] - Cos[z] - e],
{z, -Pi, Pi}]
Out[2]= If[e + Cos[x] + Cos[y] == 4, 1/(4*Pi^3*Sqrt[Abs[1 - (3 - e -
Cos[x] - Cos[y])^2]]), 0] + If[2 <= e + Cos[x] + Cos[y] <= 4,
1/(4*Pi^3*Sqrt[1 - (3 - e - Cos[x] - Cos[y])^2]), 0]
The result can be integrated numerically but it will be quite slow. A
better way is to perform some other symbolic transformations first.
Replace the integrand with an undefined function and integrate over y
and z:
In[3]:= PiecewiseIntegrate[
PiecewiseIntegrate[% /.
1/Sqrt[1 - (3 - e - Cos[x] - Cos[y])^2] :> f[e, x, y],
{y, -Pi, Pi}],
{x, -Pi, Pi}]
This way we don't attempt to integrate 1/Sqrt[1 - (3 - e - Cos[x] -
Cos[y])^2] symbolically but for each value of e we get a smooth
integrand.
In[4]:= % /.
f[e_, x_, y_] :> 1/Sqrt[1 - (3 - e - Cos[x] - Cos[y])^2] /.
Integrate[e_ + e2_, iter_] :>
Integrate[e, iter] + Integrate[e2, iter] /.
Integrate[Integrate[e_, iter_], iter2_] :>
NIntegrate[e, iter2, iter,
MaxRecursion -> 15, SingularityDepth -> 10^6, PrecisionGoal -> 2]
In[5]:= Plot[%, {e, 0, 6}, PlotDivision -> 3] // Timing
Out[5]= {3.078*Second, -Graphics-}
Large value of SingularityDepth prevents Mathematica from trying a
change of variables.
Maxim Rytin
m.r at inbox.ru