1
$\begingroup$

I am trying to approximate a Dirac Comb to use in an NDSolve calculation but when I Plot it over a large interval, it misses peaks out and they are not all of the same height (it works fine for a small interval). I have a feeling this is a sampling error but I don't know how to change this.

ϵ = 0.0001;
θz[t_] := \!\(
\*UnderoverscriptBox[\(∑\), \(i = 1\), \(120\)]
\*FractionBox[
SuperscriptBox[\(E\), 
FractionBox[\(-\ 
\*SuperscriptBox[\((t - i)\), \(2\)]\), \(4  ϵ\)]], \(2 
\*SqrtBox[\(ϵ*π\)]\)]\);
Plot[θz[t], {t, 0, 120}, PlotRange -> All]

enter image description here

Improve this question
$\endgroup$
9
  • 3
    $\begingroup$ As you say, this is Plot sampling too coarsely when plotting. If you plot with the command Plot[\[Theta]z[t], {t, 0, 120}, PlotRange -> All, MaxRecursion -> 10], you will see more spikes and they are the same height, but still a lot will be missing. Increasing further will take more time and hopefully show more spikes. Adding PlotPoints -> 100 or something should also help. But there is nothing wrong with your function, so don't worry :) $\endgroup$ – Marius Ladegård Meyer Jul 30 '18 at 11:38
  • $\begingroup$ Up to en.wikipedia.org/wiki/Dirac_delta_function and en.wikipedia.org/wiki/Distribution_(mathematics) , this plot has no sense in traditional math. Maybe, you make use of some nontraditional math. In this case, please give us references. $\endgroup$ – user64494 Jul 30 '18 at 13:51
  • $\begingroup$ Just to check -- you are aware that there is a built in function DiracComb -- so there is no need to approximate it for use in your NDSolve. $\endgroup$ – bill s Jul 30 '18 at 22:42
  • 1
    $\begingroup$ @user64494 OP clearly said that they are trying to plot an approximate Delta comb, which is a perfectly well-defined problem. An approximate delta is a smooth function, and there is nothing non-traditional about trying to plot it. $\endgroup$ – AccidentalFourierTransform Jul 30 '18 at 23:01
  • $\begingroup$ @bills I am aware, although I am modelling an MRI Scanner which is sending RF pulses so the finite width is slightly more realistic. Thanks everyone for their input, the problem is now sorted :) $\endgroup$ – Dave Bassito Jul 31 '18 at 8:33
1
$\begingroup$

You could do the following, which splits the desired horizontal range into smaller interval chunks and plots each of them separately before combining them with Show.

With[{interval = 10, range = 120},
 Show[Table[
   Plot[θz[t], {t, i, i + interval}, 
    PlotRangePadding -> None], {i, 0., range, interval}], 
  PlotRange -> {{0, range}, {0, 30}}]]

Mathematica graphics

Interestingly, this doesn't work without the option PlotRangePadding -> None. That's probably one of the many graphics bugs.

Improve this answer
$\endgroup$
0
$\begingroup$

Generally to be sure to get an accurate graph you should aim to have the initial sampling have at least one point in each pit and on each peak (or more precisely in each concave up segment and each concave down segment). Strictly speaking it is not necessary to hit every one because of the way Plot works, but it is sufficient. Then the recursive subdivision algorithm of Plot will reveal the features of the graph.

Your function has a spike at t = 1, 2,..., 120. So 241 points at roughly t = 0.0, 0.5, 1.0, 1.5,..., 120.0 should work. I say roughly, because Plot does an asymmetric subdivision. The distance between successive points will be approximately 0.5, but since the spikes are narrow, it is conceivable that some might be missed. (But it turns out that they are not. In fact, 121 or even 61 plot points works in this case, but not 62 and not 31.)

Plot[\[Theta]z[t], {t, 0, 120}, PlotPoints -> 121, PlotRange -> All, 
  MaxRecursion -> 9, ImageSize -> Medium] // AbsoluteTiming

Mathematica graphics

You'll notice some warning/error messages. That's what I get pasting your definition from the OP into Mathematica. If I convert it to InputForm, the messages go away, but it's quite a bit slower!

Mathematica graphics

A faster way to compute your function on floating-point input is the following:

θz[t_?NumericQ] := Total[E^(- (t - Range[120])^2/(4 ϵ))/(2 Sqrt[ϵ*π])];

Plot[θz[t], {t, 0, 120}, PlotPoints -> 121, PlotRange -> All, 
  MaxRecursion -> 9, ImageSize -> Medium] // AbsoluteTiming

Mathematica graphics

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.