• 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?
>

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

```

• Prev by Date: Hadamard -Sylvester Matrix Self-Similarity by substitution and reparatitioning
• Next by Date: RE: RE: [heur] Scientific Notation on Y-axis of ListPlots