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