Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

So I am looking at modelling the response of a system that is excited by multiple pulses over a period of time. The way to find the response at time t is to take the convolution between the impulse response (response of the system to an infinitely short pulse) and the multiple excitation pulses over the range [0, t]:

$$Response(t) = \int_0^t \! Excitations(\tau) ImpResp(t-\tau) \, \mathrm{d}\tau$$

So if I wanted to plot this as a function of time I would need to evaluate that integral for every point that I wanted to put on the plot and that is clearly very expensive.

My code at the moment is pretty simple, I define my gaussian pulses and the impulse response and then define a function which takes the convolution of those over [0, t] and then I plot my convolution as a function of t:

Gaussian[amp_, mean_, std_, t_] := amp*Exp[-(t - mean)^2/(2*std^2)]

GaWindow[amp_, mean_, std_, t_, w_] := Piecewise[{{Gaussian[amp, mean, std, t], mean - w <= t <= mean + w}}]

imp[t_] := 100*Exp[-t/0.01]

twopulse[t_] := GaWindow[1, 0.1, 0.0008, t, 0.008] + GaWindow[1, 0.15, 0.0008, t, 0.008]

conv[t_] := Integrate[twopulse[s]*imp[t-s],{s,0,t}]

Plot[conv[t], {t, 0.09, 0.18}, PlotRange -> All, PlotStyle -> Green]

I tried the built-in convolution function but it can only perform the convolution from negative infinity to positive infinity and when combined with step functions (to make it a function of time) it offers no speed benefits over the regular integral that I perform above.

I have had a go at discretized versions of the functions and using ListConvolve, although I was unable to recreate the correct plots from the working code I have above.

Using the PlotPoints and MaxRecursion arguments for the plot, I can increase the speed whilst sacrificing the quality of the graph, but it's nicer not to have to resort to that.

There seems to be a lot of information about optimization for infinite-range convolutions, but none that I have had the ability to apply to this multiple finite-range convolutions problem.

Thanks in advance for any suggestions on how I could improve the speed/efficiency of this.

share|improve this question
    
You can do the integration in conv symbolically, which means that it needs to be integrated only once. To do this change your definition of conv to conv[t_] = Integrate[twopulse[s]*imp[t - s], {s, 0, t}, Assumptions -> (t \[Element] Reals)]. Note the use of = rather than :=. –  Stephen Luttrell Aug 15 at 8:07
    
I'm quite new to mathematica and I find that a lot of its finer points are masked by its simplicity. Does the fact that I defined my convolution using := mean that it does the whole integration multiple times whereas using the = gives an expression that is already integrated and only requires substitution of numbers? Also, do the assumptions that you put in place speed up the operation or just avoid a warning message? Thanks a lot for your answer. –  Jack Aug 17 at 2:29
    
Actually, I think I get it. The := sets something as an expression whereas the = sets it to the evaluation of the expression. –  Jack Aug 17 at 2:52

2 Answers 2

up vote 1 down vote accepted

Way 1:

As Stephen Luttrell said in comment:

conv1[t_] := 
  Evaluate@Integrate[twopulse[s]*imp[t - s], {s, 0, t}, Assumptions -> t \[Element] Reals]

now conv1 is:

enter image description here

then plot it:

Plot[conv1[t], {t, 0.09, 0.18}, PlotRange -> All, PlotStyle -> Green, Exclusions -> None]

Way 2:

conv2[t_] := NIntegrate[twopulse[s]*imp[t - s], {s, 0, t}];
data = {#, conv2@#} & /@ Range[0.09, 0.18, 0.0002];
ListLinePlot[data, PlotRange -> All, PlotStyle -> Green]
share|improve this answer

Depending on your intended use, Convolve could provide a more useable representation of the solution than Integrate.

Please note this comes at some performance expense, however the main delay was the unevaluated use inside of Plot, which has been solved by @Chenminqi.

Your issue with only being able to work with infinite intervals can be resolved when you multiply with UnitStep[].

Leading to:

Convolve[twopulse[s] UnitStep[s], imp[s] UnitStep[s], s, t]
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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