Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I would like to calculate the convolution of a function with itself for $n$ times. Let $n$ be $10$. For $n=1$, I wont make any calculation and $f0$ is my function, when $n=2$, I will convolve $f0$ with itself and get $g=f0*f0$ when $n=3$ I will get $gg=g*f0$ etc..

As it can be seen the total number of convolutions that I need to make is $9$. At each iteration I will compute some other functions. In Matlab I do it as follows

for n=1:9
some computations related to "f0"
some computations related to "h"

I would like to do the same in Mathematica and below is the code that I typed:

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_] :=Block[{\[Omega]},q[n, B, A, f] =Function[x,Evaluate[Integrate[q[n - 1, B, A, f][\[Omega]] f[x - \[Omega]],{\[Omega], B, A}]]]]

p[n_, B_, A_, f_] :=Integrate[q[n, B, A, f][x], {x, -\[Infinity], B}] +Integrate[q[n, B, A, f][x], {x, A, \[Infinity]}]

In short $q$ is calculating the convolution and $p$ some integrals over a given $q$.

If I want to calculate for example

Sum[p[n, B, A, f0]*(k3*n + k2*p[n, B, A, f1]), {n, 1, 10}]

where $k2$ and $k3$ are some constants. Then Mathematica will make 108 convolutions and the following calculations. Instead I can make the same thing in Matlab with 9 convolutions and the following calculations.

As a result of this I decided to change my code parts as follows

q0[1] := PDF[NormalDistribution[-2, 2], x]
For[n = 2, n < 10, n++,q0[n][B_, A_] :=Convolve[q0[n - 1][B, A], f0(UnitStep[x - B] - UnitStep[x - A]), x,y]]

this code shows what I want to do but It is incorrect. I want to calculate q0[1][x], q0[2][x], q0[3][x],.. which are simply the one time, 2 times, and 3 times convolutions of $f0$ with itself (in the range from $B$ to $A$ ).

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

I would like to do some other operations for example

Sum[p[n, B, A, f0]*(k3*n + k2*p[n, B, A, f1]), {n, 1, 10}]

Thank you very much. Please let me know if something is unclear in the same for loop. How I can manage it?

share|improve this question
Look into using NestList and Nest – RunnyKine Oct 14 '13 at 13:13
@RunnyKine thank you very much for your comment, just looking now. – Seyhmus Güngören Oct 14 '13 at 23:51
@RunnyKine isnt it possible to assign the results of nestlist authomatically to some indexed functions such as $g_j$ for $j=1:10$; for the results of $9$ nested convolutions as well as the initial function $f_0$ – Seyhmus Güngören Oct 14 '13 at 23:54
up vote 1 down vote accepted

This appears to be a follow up of this question.

As suggested altering using list will provide the intermediate results. Using:

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

To achieve you goal and assuming k3,k2 are or will be defined your goal can be achieved:

Total@MapThread[#1(k3 #3+k2 #2)&,{list1,list2,Range[10]}]

18 convolutions (Matlab wil have to do 9 for each initial each function). NestList will produce 10 elements (starting condition and 9 convolutions).

I am certain there may be more efficient methods but I believe this is consistent with your aims.

share|improve this answer
I typed your code. Instead of $k2$ and $k3$, I put the numbers $2$ and $3$. I got a result as follows: – Seyhmus Güngören Oct 14 '13 at 23:48
{((Erf[(2 + a)/(2*Sqrt[2])] - Erf[(2 + b)/(2*Sqrt[2])])^7* (-Erf[(2 + a)/(2*Sqrt[2])] + Erf[(2 + b)/(2*Sqrt[2])])* (Erfc[a/2 - y/4] - Erfc[b/2 - y/4]))/(2048*E^((4 + y)^2/16)* Sqrt[Pi]) + (#1*(3*#3 + 2*#2) & ),((Erf[(-2 + a)/(2*Sqrt[2])] - Erf[(-2 + b)/(2*Sqrt[2])])^7* (-Erf[(-2 + a)/(2*Sqrt[2])] + Erf[(-2 + b)/(2*Sqrt[2])])* (Erfc[a/2 - y/4] - Erfc[b/2 - y/4]))/(2048*E^((-4 + y)^2/16)* Sqrt[Pi]) + (#1*(3*#3 + 2*#2) & ), {1 + (#1*(3*#3 + 2*#2) & ), – Seyhmus Güngören Oct 14 '13 at 23:48
2 + (#1*(3*#3 + 2*#2) & ), 3 + (#1*(3*#3 + 2*#2) & ), 4 + (#1*(3*#3 + 2*#2) & ), 5 + (#1*(3*#3 + 2*#2) & ), 6 + (#1*(3*#3 + 2*#2) & ), 7 + (#1*(3*#3 + 2*#2) & ), 8 + (#1*(3*#3 + 2*#2) & ), 9 + (#1*(3*#3 + 2*#2) & ), 10 + (#1*(3*#3 + 2*#2) & )}} – Seyhmus Güngören Oct 14 '13 at 23:49
I am sorry you find the code unhelpful. I note:(i) your intended function as described will be a function of x (indefinite as starting/initial function is PDF) + numerical values (assuming you define your limits of integratiion) from the definite integrals (convolutions).(ii) it would be helpful to know if the expected results for low numbers match between Mathematica and Matlab. I look forward to other respondents and hope they achieve your desired goal. – ubpdqn Oct 15 '13 at 8:03
@ubpdqn: As Convolve yields a function of y, your code will maybe work and produce some result, but probably not as intended. You might need to add /.y->x at the end of the Convolve functional. – Jinxed Feb 22 '15 at 23:15

Does this do what you want:

cf[f_, n_] := NestList[Convolve[#1, f, x, y] /.y -> x &, f, n]


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


cf[f0[x], 3]


 {E^(-(1/8) (2 + x)^2)/(2 Sqrt[2 \[Pi]]), E^(-(1/16) (4 + x)^2)/(
 4 Sqrt[\[Pi]]), E^(-(1/24) (6 + x)^2)/(2 Sqrt[6 \[Pi]]), E^(-(1/32) (8 + x)^2)/(4 Sqrt[2 \[Pi]])}

If you want cf[f0[x],1] + cf[f0[x],2] + ... + cf[f0[x],10] i.e. the sum of each elememt, then just do:

Total[cf[f0[x], 10]]
share|improve this answer
if I want to have $cf[f0[x], 1]+cf[f0[x], 2]+cf[f0[x], 3][f0[x], 10]$ how many convolutions will I do altogether? it is $54$ I want to do it in only $9$ operations as in MATLAB. Whenever I do one time convolution, I will make some calculations over the result that I have obtained. Then when I convolve again I will do again some calculations etc.. this will happen only 9 times. With this code I cannot do it. – Seyhmus Güngören Oct 14 '13 at 23:30
@SeyhmusGüngören, see my edit. – RunnyKine Oct 15 '13 at 1:06
Thank you very much. Is it possible to assign each output of Nestlist as a function authomatically? For example: {E^(-(1/8) (2 + x)^2)/(2 Sqrt[2 [Pi]]),(E^(-(1/16) (4 + x)^2) (-Erfc[A/2 - x/4] + Erfc[B/2 - x/4]))/(8 Sqrt[[Pi]])} is the result of cf[f0[x], 1]. I want fx1=E^(-(1/8) (2 + x)^2)/(2 Sqrt[2 [Pi]]) and fx2=(E^(-(1/16) (4 + x)^2) (-Erfc[A/2 - x/4] + Erfc[B/2 - x/4]))/(8 Sqrt[[Pi]]) , or even better fx[1]=.. and fx[2]=.. – Seyhmus Güngören Oct 15 '13 at 9:50
meanwhile cf[f0[x], 9] is still running)) – Seyhmus Güngören Oct 15 '13 at 12:38
@SeyhmusGüngören, Note that I used NestList in my definition so that my output for cf[f0[x], n] is a list that contains iterations from 1 to n. If you want just individual values use Nest. Second, I updated my definition sometime yesterday to include /. y -> x. On my machine cf[f0[x], 10] takes 1 second to compute. – RunnyKine Oct 15 '13 at 13:05

Your Answer


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.