Message Boards Message Boards

[✓] Plot this gaussian function?

GROUPS:

Consider the following code:

    gauss[x_] := E^(-x*x)

    sinc[x_] := Sinc[x]

    Plot[gauss[x], {x, -5, 5}, PlotRange -> All]

    Plot[sinc[x], {x, -5, 5}, PlotRange -> All]

    Plot[Convolve[gauss[x], sinc[x], x, y], {y, -5, 5}, PlotRange -> All]

The first two plots look good, but it's been evaluating the final plot for an hour now, and I don't remember it working when I tried a few months ago. Incidentally, this doesn't evaluate to any close form:

Convolve[gauss[x], sinc[x], x, y]

What strategies can I try to get a plot of it?

POSTED BY: Joe Donaldson
Answer
3 months ago

Will a quick rough approximation help? (Tweak the constants if needed)

ListPlot[Table[{y, NIntegrate[gauss[x] sinc[y-x],{x,-40,40}]},{y,-40,40,1/4}],PlotRange->All,Joined->True]

enter image description here

POSTED BY: Bill Simpson
Answer
3 months ago

Yes, that definitely helps. Thanks.

POSTED BY: Joe Donaldson
Answer
3 months ago

Joe,

While Bill's solution is excellent, another alternative is to do the convolution digitally:

gd = gauss@Range[-40, 40, .25];

sd = sinc@Range[-40, 40, .25];

ListLinePlot[ListConvolve[gd, sd, (Length[gd] - 1)/2], PlotRange -> All]

getting the same plot.

Regards

POSTED BY: Neil Singer
Answer
3 months ago

Thanks. I'm grateful for the solutions, but I'm still curious what was wrong with my attempt?

I understand that the convolution entails integration which Mathematica can't put into closed form, so it's a tough problem, but I'd thought in that case Mathematica would recognize the difficulty and switch to numerical estimation? In other words, I'd thought Mathematica would automatically do something like what you've both shown me how to do manually.

Admittedly, I haven't studied the Plot options lately, so I might have overlooked something in there.

POSTED BY: Joe Donaldson
Answer
3 months ago

Edited version -- I grabbed the wrong plot and forgot the Fourier Constant

Joe,

Mathematica can Solve it exactly by using the FourierTransform:

In[81]:= ga = FourierTransform[E^(-x*x), x, w]

Out[81]= E^(-(w^2/4))/Sqrt[2]

In[82]:= sa = FourierTransform[Sinc[x], x, w]

Out[82]= 1/2 Sqrt[\[Pi]/2] (Sign[1 - w] + Sign[1 + w])

In[84]:= ans = Sqrt[2*Pi]*InverseFourierTransform[ga*sa, w, x]

Out[84]= 1/2 E^-x^2 \[Pi] (Erf[1/2 - I x] + Erf[1/2 + I x])

In[86]:= Plot[ans, {x, -40, 40}, PlotRange -> All]

To get:

enter image description here

This relies on the relationship that the Convolution of two functions is the InverseFourierTransform of the multiplication of their FourierTransforms (with the appropriate constant which depends on the definition of Fourier Transform that you use -- in Mathematica it is Sqrt[2*Pi] unless you change options).

POSTED BY: Neil Singer
Answer
3 months ago

Oops, I forgot the scale factor. I'll fix it later when I get back to my machine.

--Fixed above.

POSTED BY: Neil Singer
Answer
3 months ago

Please don't worry about the scale-factor on my account. I understand the relation between fourier-xform and convolution, and am fairly familiar with the Gaussian and Sinc functions.

Why doesn't this direct method yield the results you obtained indirectly by using FourierTransform?

gauss[x_] := E^(-x*x)

sinc[x_] := Sinc[x]

Convolve[gauss[x], sinc[x], x, y]

Why did my original attempt fail? I want to understand why the workarounds are necessary instead. (Mostly I wanted to get a plot, but now that I have a good plot, I'm curious why the direct method fails.)

POSTED BY: Joe Donaldson
Answer
3 months ago

Joe,

  1. You raise a good point that it was fairly easy for me to get the closed form solution to the Convolution using the Fourier Transform -- I am surprised that the Convolve function does not try to do this as part of its internal algorithm. Maybe someone at Wolfram can shed light on this.

  2. I fixed the errors I made so this thread will be more useful to people searching in the future.

  3. One other point: -- even if Convolve would have been able to solve for the closed form solution, I would recommend against your Plot command:

    Plot[Convolve[gauss[x], sinc[x], x, y], {y, -5, 5}, PlotRange -> All]

Doing the plot this way forces Mathematica to recompute the closed form solution each time a point is plotted and could be very slow. You are better off doing:

ans= Convolve[...]
Plot[ans, {y, -5, 5}, PlotRange -> All]

This solves the (difficult) calculation once and then substitutes in values.

Regards

POSTED BY: Neil Singer
Answer
3 months ago

Thanks, that is helpful to me.

POSTED BY: Joe Donaldson
Answer
3 months ago

Group Abstract Group Abstract