I would like to compute the integral

$$q_n(x)=\int_{B}^A q_{n-1} (\omega)f(x-\omega)\mathrm{d} \omega,\quad q_1=f,\quad n\geq 1.$$

for a given density function $f$. Basically it is the convolution of a pdf with itself for $n$ times in the range $A$ and $B$. when $n=1$, I simply plot the original pdf $f$, if $n=2$, then I convolve it with itself only once and plot etc..

Eventually I would like to plot the densities $q_n$ for $n=1,2,3..,7$ in the same figure.

My code with the help of some people from this forum is

f0[x_] := PDF[NormalDistribution[-2, 2], x]
f1[x_] := PDF[NormalDistribution[2, 2], x]

q[1, B_, A_, f_] := f

q[n_ /; n > 1, B_, A_, f_][x_] := Module[{\[Omega]}, 
tempIntegrate[Evaluate[q[n - 1, B,A, f]][\[Omega]] f[x - \[Omega]], {\[Omega],B, A}]]

p[n_, B_, A_, f_] := 
tempIntegrate[q[n, B, A, f][x], {x, -\[Infinity], B}]+tempIntegrate[q[n, B, A, f][x], {x,A, \[Infinity]}] //. 
{s_ tempIntegrate[b_, a__] :>tempIntegrate[s b, a], 
tempIntegrate[tempIntegrate[b_, a__], c__] :>tempIntegrate[b, a, c]} /. 
tempIntegrate -> NIntegrate

I can not

Plot[q[2, -2, 2, f0][x],{x,-10,10}] 

correctly and I cannot see the analytical formulation when I type

q[2, B, A, f0][x]

Could you please help me? Thank you very much.


My eventual intention was to calculate

$$P[n]=\int_{(-\infty,B)\cup (A,\infty)} q_n (x) \mathrm{d}x$$

and p[2, -2, 2, f0] provides me a numerical value but I am also interested in some analytical values.

To do this I tried to change ->NIntegrate to Integrate but then the code didnt function properly.

  • $\begingroup$ Maybe you are not finished editing yet, but what is your definition of tempIntegrate? $\endgroup$ Commented Oct 10, 2013 at 12:36
  • $\begingroup$ If you replace tempIntegrate with Integrate, then Plot[q[2, 1, 3, f0][y], {y, 0, 2}] works for me. Note that Evaluate does nothing in your code. Evaluate is often used in Plot to ignore the fact that Plot does not evaluate its first argument. Integrate however does not hold any arguments. In addition, the Evaluate would only override this if it would surround the entire first argument. I suggest removing the Evaluate :). Note that in the Plot code I said that worked for me, the function gets integrated a lot of times. I think that is not necessary $\endgroup$ Commented Oct 10, 2013 at 12:43
  • $\begingroup$ @JacobAkkerboom I tried your suggestions before but the results which I got was incorrect. Please give me 5 minutes and I get some correct numerical values for comparison. $\endgroup$ Commented Oct 10, 2013 at 12:50
  • $\begingroup$ 0.5227501319481026` , 0.2973909442488435` , 0.11261520561472473` , 0.04211024363961929` , 0.01573950408880347` are the true values for $p[n, -2, 2, f0]$ for $n=1,...,5$ $\endgroup$ Commented Oct 10, 2013 at 12:55
  • $\begingroup$ @JacobAkkerboom i did your suggestion, removed evaluate and changed tempintegrate to Integrate and eventually p[3, -2, 2, f0] gave me nothing)) $\endgroup$ Commented Oct 10, 2013 at 13:08

3 Answers 3


Maybe this will do what you want

q[1, B_, A_, f_] := f
q[n_ /; n > 1, B_, A_, f_] :=
  q[n, B, A, f] =
      q[n - 1, B, A, f][\[Omega]] f[x - \[Omega]], {\[Omega], B, A}]]

Plotting q2 will now be very fast, as Mathematica will evaluate the integral for q2 to a simpler form, even before x is known, which greatly speeds up computation. That is, in this case you do not evaluate an integral for every x you wish to plot.

However, it turns out Mathematica is not able to find a nice expression for q3. We can still plot it, but at every step Mathematica tries to evaluate a difficult integral analytically, so this takes long and is suboptimal. Maybe I can improve it :).


Plot[q[2, 1, 3, f0][x], {x, 0, 2}]

Mathematica graphics

Note that the definition of q[2,1,3,f0] is now


Function[x$,(E^(-(1/16) (4+x$)^2) (-Erf[(2-x$)/4]+Erf[(6-x$)/4]))/(8 Sqrt[\[Pi]])]

so that we indeed see that here the Integral has been analytically "solved"/"simplified"/"evaluated".

However, again, this does not work for q3. The definition of that becomes

q[3, 1, 3, f0]
(*output in InputForm ;)*)

Function[x$, Integrate[(E^(-(2 + x$ - \[Omega])^2/8 - (4 + \[Omega])^2/16)*
    (-Erf[(2 - \[Omega])/4] + Erf[(6 - \[Omega])/4]))/(16*Sqrt[2]*Pi), {\[Omega], 1, 3}]]

where we see that the function Integrate is still present.

  • $\begingroup$ Note that to get values for q4 you have to use the function N, that is, write for example q[4, 1, 3, f0][2] // N. $\endgroup$ Commented Oct 10, 2013 at 13:21

There may be a better way to approach this problem. If you transform your function f[t] into the frequency domain ft[w], convolution becomes multiplication, and the iterative convolution of f[t] with itself becomes raising ft[w] to a power. Here is a simple example where

f[t_] := Exp[-t^2/2];
ft[w_] := FourierTransform[f[t], t, w]

Now you can plot the powers of ft[w]:

Plot[Table[ft[w]^n, {n, 1, 5, 1}], {w, -3, 3}]

enter image description here

which are the transforms of the time functions. You can return to the time domain by inverse-transforming, for example, at the 5th step,

InverseFourierTransform[ft[w]^5, w, t]


Of course you can replace this f[t] with your f0[t] without any problems. Below, the OP asked for how to do this when the range of the function is restricted. This can be handled by truncating the function. Here's a truncated normal distribution:

f1[t_, a_, b_] := UnitStep[t - a] UnitStep[b - t] PDF[NormalDistribution[-2, 2], t];

This can be transformed (almost) as easily as above:

f1t[w_, a_, b_] = FourierTransform[f1[t, a, b], t, w];
Plot[Table[N[Abs[f1t[w, -3, 3]]]^n, {n, 1, 5, 1}], {w, -3, 3}, PlotRange -> All]

enter image description here

  • $\begingroup$ Thanks for the answer. It is not a bad idea but my convolutions is in some prespecified range $A$ and $B$. how can I do in in a parametric way? $\endgroup$ Commented Oct 10, 2013 at 15:00
  • $\begingroup$ Define your function to be whatever you want in the region A to B, and zero outside. UnitStep is one way to do this. $\endgroup$
    – bill s
    Commented Oct 10, 2013 at 15:01
  • $\begingroup$ thank you very much for the edit. I have bigger problem but the idea of frequency domain sounds very interesting, After some progress I will get back to you. $\endgroup$ Commented Oct 10, 2013 at 15:43
  • $\begingroup$ I dont get the result of inverse Fourier transform $\endgroup$ Commented Oct 15, 2013 at 20:21
con[f_, g_, x_, l_, u_] := 
 Integrate[(f /. x -> x ) (g /. x -> (y - x)), {x, l, u}];
func[u_, ll_, ul_] := 
 Nest[con[#, PDF[NormalDistribution[], x], x, ll, ul] &, 
  PDF[NormalDistribution[], x], u];

Using integration limits of -1 and 1:

tab = Table[func[j, -1, 1], {j, 5}];
Plot[tab, {y, -3, 3}, PlotLegends -> Range[5]]

enter image description here gives:


To allow use of f0[x] and f1[x] and incorporating alternative suggested by bills):

q[f_, g_, x_, b_, a_, n_] := 
 Nest[Convolve[#, f (UnitStep[x - b] - UnitStep[x - a]), x, y] &, g, 

Applying to f0[x] (note must put as expression of x not just f0]:

qtab = Table[q[f0[x], f0[x], x, -1, 1, j], {j, 5}];

and plotting:

Plot[qtab, {y, -4, 4}, PlotLegends -> Range[5]]

enter image description here

  • $\begingroup$ It might be easier to use the Convolve command, which does more or less what your con does. $\endgroup$
    – bill s
    Commented Oct 11, 2013 at 13:01
  • $\begingroup$ @bills thank you for the feedback. Convolve uses -Infinity to Infinity as limits of integration (as far as I understand the documentation). I could not find options to alter the limits of integration, hence custom function. $\endgroup$
    – ubpdqn
    Commented Oct 11, 2013 at 13:14
  • $\begingroup$ You can limit the range of integration by changing the function (Unitstep for instance) because any given convolution has fixed endpoints. But your custom function is fine too... $\endgroup$
    – bill s
    Commented Oct 11, 2013 at 14:16
  • $\begingroup$ @bills thank you for this, another option for my future reference $\endgroup$
    – ubpdqn
    Commented Oct 11, 2013 at 14:17
  • $\begingroup$ I am not able to use it in my program. I replaced PDF[NormalDistribution[] parts with either $f_0$ or $f_0[x]$ and both gave an error. I would like to have a form q[n,B,A,f] for which I can select an $f$ too, for example $f_0$ as given in the question. How one should modify the code? $\endgroup$ Commented Oct 11, 2013 at 14:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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