Recursive piecewise integral relation with piecewise base case?

How is this recursive formula

$$f_{n+1}(z) = \int_0^1 f_{n}(z-y)\,{\rm d}y$$

implemented in Mathematica? The base case is

$$f_1(z) = \begin{cases} 1 & 0\leq z\leq 1 \\ 0 & \text{otherwise} \end{cases}$$

I have tried setting up a recursive integral relation, but my syntax must not be right.

• Is there a base case, perhaps for $f_0$ or $f_1$? – Pillsy Nov 3 '15 at 21:33
• @Pillsy no idea how to define piecewise base cases. – martin Nov 3 '15 at 21:35
• go ahead and write it out in LaTeX, then. – Pillsy Nov 3 '15 at 21:38
• @Pillsy $$f_{X+Y}(z) = \int_{0}^1 f_X(z-y)\,{\rm d}y = \cases{z&0 \le z\le 1\cr 2-z&1 \le z\leq 2\cr 0& otherwise}$$ $$f_{X+Y+Z}(z) = \int_{0}^1 f_{X+Y}(z-y)\,{\rm d}y = \cases{ \frac{z^2}{2}&0 \le z\le 1\cr -\frac 3 2- {z}^{2}+3\,z&1 \le z\le 2\cr \frac 92-3\,z+\frac{z^2}{2}&2 \le z\leq 3\cr 0& otherwise}$$ – martin Nov 3 '15 at 21:42
• f[1] = Integrate[ PDF[UniformDistribution[{0, 1}], z - y], {y, 0, 1}] /. z -> y; f[n_] := Integrate[f[n - 1] /. y -> z - y, {y, 0, 1}] /. z -> y; f[3] // Simplify ? – Dr. belisarius Nov 3 '15 at 21:50

f[1] = Integrate[PDF[UniformDistribution[{0, 1}], z - y], {y, 0, 1}] /. z -> y;
f[n_] := f[n] = Integrate[f[n - 1] /. y -> z - y, {y, 0, 1}] /. z -> y // Simplify;
f[3]


Adapting once again Leonid's solution from here,

f[1, z_] := UnitBox[z - 1/2];
f[n_Integer, z_] := Module[{zl, t},
Set @@ Hold[f[n, zl_],
Simplify[Convolve[UnitBox[t - 1/2], f[n - 1, t], t, zl]]];
f[n, z]];

f[4, z]


$\displaystyle\begin{cases} -\frac16(-4+z)^3&3\le z<4\\ \frac{z^3}{6}&0<z\le1\\ \frac23-2z+2z^2-\frac{z^3}{2}&1<z\le2\\ -\frac{22}{3}+10z-4z^2+\frac{z^3}{2}&2<z<3\\ 0&\mathtt{True} \end{cases}$

Plot[f[6, t], {t, -2, 2}]


Now, here's the surprise: there's an even shorter implementation for f[]!

f[n_Integer, z_] := BSplineBasis[n - 1, z/n]

• I believe that your final observation can be derived quickly in Fourier space. – bbgodfrey Nov 4 '15 at 14:03
• That's correct; in fact, one of the ways to present the theory of B-splines is to consider them as repeated convolutions of a boxcar function. – J. M.'s technical difficulties Nov 4 '15 at 14:04

@Winther's solution is particularly fast from here

pwf[z_] := Piecewise[{z[[#]], # - 1 <= y < #} & /@ Range@Length@z]
iidf[n_] := With[{nn = n}, ffunc = Table[If[i == 1, 1, 0], {i, 1, n}];
Do[temp = ffunc; temp[[1]] = Integrate[ffunc[[1]], {z, 0, z}];
Do[temp[[k]] = Integrate[ffunc[[k - 1]], {z, z - 1, k - 1}] +
Integrate[ffunc[[k]], {z, k - 1, z}];, {k, 2, Floor[(i + 1)/2]}];
Do[temp[[k]] = temp[[i - k + 1]] /. z -> i - z;,
{k, Floor[(i + 1)/2] + 1, i}]; ffunc = temp;, {i, 2, n}];
pwf@ExpandAll[ffunc] /. z -> y]

Plot[Evaluate[iidf@6], {y, 0, 6}, PlotPoints -> 400]


Along with @J.M.'s superfast solution, and this nice little identity, where for $X_j \text{ iid},$ uniformly distributed on $[0,1],$

$$\dfrac{1}{n!} \left\langle n \atop k \right\rangle = P\left(\sum_{j=1}^{n}X_j\in[k,k+1]\right)$$

we can get eg eulplot[6, 2], eulplot[12, 5]:

eulerian[k_, n_] :=
CoefficientList[(1 - x)^(n + 1) PolyLog[-n, x]/x, x][[k + 1]]

iid[k_, n_] := eulerian[k, n]/n!

eulplot[n_, pt_] := With[{aa = Piecewise[SortBy[(BSplineBasis[n - 1,
x/(n)] // PiecewiseExpand)[[1]], Last@# &]]},
Show[Plot[Evaluate[aa], {x, 0, n}, PlotPoints -> 1000],
Plot[Evaluate[aa[[1, pt + 1]]], {x, pt, pt + 1}, Filling -> Axis,
PlotRange -> {{0, Automatic}, {0, Automatic}}], Frame -> True,
LabelStyle -> Black, PlotLabel -> StringJoin["A=", ToString[

• There is an undocumented function for the Eulerian: NumberTheoryEulerianNumber[n, k]`. Note that the argument order in your version is reversed. – J. M.'s technical difficulties Nov 7 '15 at 7:27