Improving the speed/efficiency of multiple (finite range) convolutions for a causal response plot

Question

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.

Answer 1

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:

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]

Answer 2

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]