Calling Correct Function for Plotting DiracDelta

I am wondering what is the correct function in Mathematica to plot the true impulse function, better known as the DiracDelta[] function. When using this inside of a function or just the function itself when plotting, it renders output = zero. Quick example:

Plot[DiracDelta[x], {x,-1,1}]

I am wondering, is this the correct delta function which is infinite in height at zero and zero everywhere else. I have seen other functions such as KroneckerDelta[], but this seems to do the same exact thing. The code is below:

eqn1wb := (2 ((I Pi f) Exp[(-I) 2 Pi f]) ((1/2) DiracDelta[f-2] +
DiracDelta[f+2]))/(1 + I 2 Pi f);
wb1 = Plot[Re[eqn1wb], {f, -5, 5},
PlotRange :> All, PlotStyle :> {Thick, Red},
AspectRatio :> 1,
GridLines :> {{-4, -2, 2, 4}, {1.0, 0.5, -0.5, -1.0}},
PlotLabel :> "Re(w) - Frequency Response",
AxesLabel :> {"f", "W(f)"}
]

To create the plot you could replace any occurrence of DiracDelta[a] with something like 10000 UnitStep[1/10000 - a^2]], so for example to plot

f[x_] := DiracDelta[x - 2] + DiracDelta[x + 2]

you could do something like

Plot[Evaluate[f[x] /. DiracDelta[a_] :> 10000 UnitStep[1/10000 - a^2]],
{x, -4, 4}, Exclusions -> None, PlotPoints -> 800] Note that for Mathematica to see the discontinuities you need to increase the number of plot points. The number of points needed will depend on the plot range, so you might have to tweak that.

• A large number of PlotPoints isn't necessary. There is a nice possibility to specify additional points with PlotPoints->{Automatic, {-2, 2}} (ref). Then recursion algorithm will plot it fine. Oct 28 '15 at 19:13

As others already have written, the Dirac delta is not a real function and it can't be plotted. Other programs that claim to plot it just fake it.

Having said that, you can roll a diracDelta of your own, that more or less mimics the Dirac Delta's behavior but is still continuous. Advantage with respect to Heike's solution is that we don't need to increase the number of PlotPoints. I'll use the NormalDistrution with a very small standard deviation:

diracDelta[x_] = PDF[NormalDistribution[0, 1/100], x];

eqn1wb := (2 ((I Pi f) Exp[(-I) 2 Pi f]) ((1/2) DiracDelta[f - 2] +
DiracDelta[f + 2]))/(1 + I 2 Pi f);

wb1 = Plot[
Re[eqn1wb] /. DiracDelta -> diracDelta // Evaluate, {f, -5, 5},
PlotStyle :> {Thick, Red}, AspectRatio :> 1,
GridLines :> {{-4, -2, 2, 4}, {1.0, 0.5, -0.5, -1.0}},
PlotLabel :> "Re(w) - Frequency Response",
AxesLabel :> {"f", "W(f)"}, PlotRange -> {Automatic, 1.5}] DiracDelta's are sometimes plotted by using arrows. Of course we can do that as well. In this case we can set the function at zero everywhere and add arraows in the Epilog part of the plot.

Plot[0, {f, -5, 5}, PlotStyle :> {Thick, Red}, AspectRatio :> 1,
GridLines :> {{-4, -2, 2, 4}, {1.0, 0.5, -0.5, -1.0}},
PlotLabel :> "Re(w) - Frequency Response",
AxesLabel :> {"f", "W(f)"}, PlotRange -> {Automatic, 1.},
Epilog -> {Red, Thick,Arrow[{{-2, 0}, {-2, 1}}], Arrow[{{2, 0}, {2, 1}}]}] I believe the convention is to scale the lengths of the arrows in proportion to the factors in front of the DiracDelta. The above can be easily extended to do that.

Factors can be found using Coefficient:

Coefficient[eqn1wb, DiracDelta[-2 + f]]

(*
==> (I E^(-2 I f \[Pi]) f \[Pi])/(1 + 2 I f \[Pi])
*)

Ratio's of coefficients:

Coefficient[eqn1wb, DiracDelta[2 + f]]/Coefficient[eqn1wb, DiracDelta[2 - f]]

(*
==> 2
*)

Locations of the Dirac delta's:

Cases[eqn1wb, DiracDelta[a__] :> Solve[a == 0, f], Infinity]

(*
==> {{{f -> 2}}, {{f -> -2}}}
*)

I'm not going to automate it all, but the idea is clear. In the end you'll get something like: The Dirac delta, $\delta(x)$ is zero everywhere except at zero, and has an integral of 1 over $\mathbb{R}$. It is not really a function in the true sense and equating $\delta(0)=\infty$ is a rather loose definition; it should technically be considered as a distribution or a delta measure. Mathematica's implementation of DiracDelta remains unevaluated at 0 and has all the other properties.

DiracDelta /@ {-1, 0, 1}
Out= {0, DiracDelta, 0}

Integrate[DiracDelta[x] , {x, -Infinity, Infinity}]
Out= 1

This explains why you weren't getting any result with Plot.

In the off-chance that you were trying to plot the impulse response of a discrete-time system, then you might be interested in the DiscreteDelta function:

ListPlot[DiscreteDelta /@ Range[-5, 5], Filling -> Axis, PlotMarkers -> {Automatic, 10}] You can also reproduce your MATLAB "plot" of the Dirac delta function in Mathematica. Let me first note that you're not using their dirac() function to produce the output you've shown below. In fact, they define it "loosely":

function Y = dirac(X)
%DIRAC  Delta function.
%    DIRAC(X) is zero for all X, except X == 0 where it is infinite.

Y = zeros(size(X));
Y(X == 0) = Inf;

and this would not have given you the plot you showed, because of the Inf. You have probably replaced Inf with 10000 or written a similar function. Heike, image_doctor and Sjoerd have shown you ways of plotting it. You can also roll your own, a la MATLAB style as:

dirac[x_] := Piecewise[{{10000, x == 0}, {0, True}}]
ListLinePlot[{#, dirac@#} & /@ Range[-1, 1, 1/1000], Frame -> True, Axes -> False] Now use this definition to redefine DiracDelta inside a Block and you can use your original code with some slight modifications:

wb1 = Block[{DiracDelta = Piecewise[{{10000, # == 0}, {0, True}}] &},
ListLinePlot[Re[eqn1wb] /. f -> # & /@ Range[-5, 5, 1/100],
DataRange -> {-5, 5}, PlotRange :> {All, {0, 1}},
PlotStyle :> {Thick, Red}, AspectRatio :> 1,
GridLines :> {{-4, -2, 2, 4}, {1.0, 0.5, -0.5, -1.0}},
PlotLabel :> "Re(w) - Frequency Response",
AxesLabel :> {"f", "W(f)"}]
] If I were to use a schematic to show the delta function, I would do it similar to Sjoerd's arrow plots.

• Okay, so what your saying is, is that it cannot be plotted over $\mathbb{R}$? I can accomplish it in another program, but wanted to do it in Mathematica. Mar 26 '12 at 8:26
• Can you show your output plot that you generated in your other program?
– rm -rf
Mar 26 '12 at 8:30
• You can define your own Dirac Delta function which has the behaviour that you would like. DiracD[x_] := DiracDelta[x] /. DiracDelta -> Infinity . You can swap infinity for any large value that suits your application if you wish. Mar 26 '12 at 8:36
• @R.M: Here is the output graph of the Delta function. i.stack.imgur.com/VSEr9.png Mar 26 '12 at 9:00
• @image_doctor: This seems to not work. I used Plot[DiracD[x], {x,-1,1}, PlotStyle:> Thick] to get this: i.stack.imgur.com/Yh6r8.png Mar 26 '12 at 9:03

Building on Sjoerd C. de Vries's solution, I wrote this little function to make the task easy for anyone interested in plotting Dirac's Delta as arrows in a Plot command in Mathematica. The function receives an equation and a variable, and returns a list of {Arrow[]}'s in the correct location of the Diract Deltas, that can be included in the Epilog of a Plot command.

ArrowsDeltaFunction[eqn_, x_] := Module[{xsubs, listDeltas, coefDeltas, locationDeltas},
xsubs = (x /. Cases[eqn, DiracDelta[a__] :> Solve[a == 0, x], Infinity]) /. x -> {};
listDeltas = DiracDelta[x - x0] /. x0 -> xsubs;
coefDeltas = Flatten[Coefficient[eqn, listDeltas]];
locationDeltas = Flatten[xsubs];

Arrow[Table[{{locationDeltas[[i]], 0}, {locationDeltas[[i]],
coefDeltas[[i]]}}, {i, 1, Length[locationDeltas]}]]
]

So you can use this function throughout your code by including it in the Plot Epilog, like this:

eqn1 := x + (3/4) DiracDelta[x - 0.2] + (\[Pi]/3) DiracDelta[x + 0.5]

Plot[eqn1, {x, -5, 5},
PlotStyle -> {Thick, Red},
Epilog -> {Thick, Red, ArrowsDeltaFunction[eqn1, x]}
]

You can see the arrows appear in the "gaps" of the function where the Delta's exist, and their height correspond to the coefficient of the Delta. The code works even if the function does not have Deltas (no arrows are plotted).

Edit: This can have nice applications for teachers to explain Fourier Transforms, its properties, etc.

(* Time Domain Signal (in this case complex) *)
Ft :=  (1 - I/2) Exp[-I t] +  Cos[2 t]  + 2 Sin[4 t]

(* Fourier Transform *)
Fw := FourierTransform[Ft, t, w]
ReFw := ComplexExpand[Re[Fw]]
ImFw := ComplexExpand[Im[Fw]]

(* Plot *)
Plot[{ReFw, ImFw}, {w, -5, 5},
PlotStyle -> {Thick, Red, Blue},
PlotRange -> 3,
AxesLabel -> {"w", "F(w)"},
PlotLegends -> LineLegend[{Red, Blue}, {"Real part", "Imag part"}],
Epilog -> {Thick,
Red, ArrowsDeltaFunction[ReFw, w],
Blue, ArrowsDeltaFunction[ImFw, w]
}
] • In the case of functions like eqn=DiracDelta[2x-2], the ArrowsDeltaFunction function attempts to extract the coefficient of DiracDelta[x-1], instead, even though no such term exists. This results in some arrows incorrectly having 0 height. I implemented a fix, which I'll post below. May 18 '16 at 10:11
• ArrowsDeltaFunction[eqn_, x_] := Module[{xsubs, listDeltas, coefDeltas, locationDeltas}, xsubs = (x /. Cases[eqn, DiracDelta[a__] :> Solve[a == 0, x], Infinity]) /. x -> {}; locationDeltas = Flatten[xsubs]; Arrow[Table[{{locationDeltas[[i]], 0}, {locationDeltas[[i]], (eqn/ DiracDelta[x - locationDeltas[[i]]]) /. x -> locationDeltas[[i]]}}, {i, 1, Length[locationDeltas]}]]] May 18 '16 at 10:12

At your own risk, you can redefine the DiracDelta function to have the behaviour you want:

Unprotect@DiracDelta;