# The curious case of Convolve function

Bug Introduced in Version 11 or earlier and persisting through 12.0

For integer values of $$n$$, I was trying to calculate the convolution of $$t^{-n}$$ with itself. So I wrote this:

Table[Convolve[t^-n UnitStep[t], t^-n UnitStep[t], t, x], {n,5}]


which resulted in a sequence of functions whose general form was like this: $$\frac{b_n+c_n\ln{x}}{x^{2n-1}}\text{u}(x)\qquad n=1,2,3,...$$ where $$b_n\le0$$ and $$c_n>0$$ are constant values. I couldn't find a pattern in these constants. So I tried to evaluate the general term of convolution in Mathematica. I defined f as:

f[t_,n_]:= Piecewise[{{t^-n, Element[n,Integers] && n>0 && t>0}, {0,True}}]


and

g[x_,n_]:= Evaluate[Convolve[f[t,n],f[t,n],t,x]]


to get: $$\frac{\text{u}(t)}{t^n}*\frac{\text{u}(t)}{t^n}=\frac{\Gamma(1-n)\sqrt{\pi}}{\Gamma(3/2-n)}\left(\frac 2t\right)^{2n-1}\text{u}(t)$$ But this doesn't add up. For example, if you write

h[t_,n_]:= (2/t)^(2n-1) Sqrt[Pi] Gamma[1-n] / Gamma[3/2-n];
Limit[h[t, n], n -> 1]


It gives the result as $$\pm\infty$$, which is surely not equal to the convolution of $$\frac 1t\text{u}(t)$$ with itself. (It was $$\frac {2\text{u}(t)}t\ln{t}$$ by the way).

So the question is, what am I missing here? Why Mathematica gives two completely different results for two (seemingly) same expressions?

• The form of the result suggests that the Gauss hypergeometric theorem is being applied improperly here. Let me think about it... – J. M. is away Jan 3 '17 at 8:34
• FWIW, FullSimplify[g[x,n]] evaluates to ComplexInfinity symbolically (no need to take the limit of $n\to1$). The reason is that Gamma[1-n]==ComplexInfinity for all integer n. – AccidentalFourierTransform Jan 3 '17 at 12:13
• @AccidentalFourierTransform This has been answered using the Pareto distribution defined such that the area under the curve is not infinite, see – Carl Mar 11 '18 at 4:50
• While I agree this looks like a bug based on your cited question, math.stackexchange.com/questions/2082959, this community's convention is to wait for a case-report from WRI confirming something's status as a bug before we attach the tag. If you haven't already reported it to them you can do so, attach the case number, and then put the tag back on. – b3m2a1 Mar 11 '18 at 16:46
• @b3m2a1 I am not familiar with the process. In fact, I have no idea what you meant by case-report from WRI. The post is made over a year ago, but the bug is still there even in the new version. So it looks like nobody cares, and since I am quite busy with other stuff, I think my best course of action is to just leave it be. – polfosol Mar 11 '18 at 16:58

I don't know why Mathematica returns the incorrect result, but it may be helpful to have the correct result so as to find the origin of the bug.

To find the correct result, one may generate the first few convolutions for $n=1,2,\cdots,20$ and then use FindSequenceFunction to get the general formula for arbitrary $n$. The result is

-((4^(-1 + n) x^(1 - 2 n) Gamma[-(1/2) + n] (Log[4] - 2 Log[x] + PolyGamma[0, -(1/2) + n] - PolyGamma[0, n]))/(Sqrt[\[Pi]] Gamma[n]))


which, for $n=1,2,3,\cdots$ evaluates to

(2 Log[x])/x, (-2 + 4 Log[x])/x^3, (-7 + 12 Log[x])/x^5, ...


The Latex formula for arbitrary $n$ reads $$\frac{4^{n-1} x^{1-2 n} \Gamma \left(n-\frac{1}{2}\right) \left(-\psi ^{(0)}(n)+\psi ^{(0)}\left(n-\frac{1}{2}\right)-2 \log \left(\tfrac12x\right)\right)}{\sqrt{\pi }\ \Gamma (n)}$$

• That formula of yours can be compacted further: 2 x^(1 - 2 n) Binomial[2 n - 2, n - 1] (HarmonicNumber[n - 1] - HarmonicNumber[2 n - 2] + Log[x]) – J. M. is away Jan 3 '17 at 13:45
• I made a post about this on MathSE. It turned out the results that I thought were wrong are actually RIGHT! @J.M. – polfosol Jan 4 '17 at 6:25
• @polfosol indeed, if you evaluate Assuming[x>0,Integrate[t^-1 HeavisideTheta[t](x-t)^-1 HeavisideTheta[x-t],{t,-Infinity, Infinity}]] you get Integrate::idiv: Integral does not converge on {-Infinity,Infinity}. >> which confirms your result that the integral diverges! I agree with you that this is clearly a bug of MM. – AccidentalFourierTransform Jan 4 '17 at 13:39

Here is the email I just received from Wolfram product support regarding this issue:

RE: [CASE:4209568] Your Wolfram Product Feedback

Hello *******,

Thanks for contacting Wolfram Technical Support.

It does appear that it is not behaving properly. I have forwarded an issue report to our developers with the information you provided. I also included your contact information in my report.

Furthermore, I would request you to report such issues to the Wolfram Technical Support (http://support.wolfram.com) so that it can be forwarded to the developer team.

Please do not hesitate to contact us if you have any further queries.

So it does look like the bug will hopefully be resolved in the next version.

In the meantime, it's still unclear for me that why would someone remove the [bug] tag from this question?

I don't know if this is related to the same underlying problem or if it is a different problem with Convolve. I also realize this is not an answer per se, but it's too large to fit the limits of a comment. I'm posting it here in case it may be related somehow to the strange behavior the OP encountered. When I was working with Convolve a while ago, I identified the following strange behavior. It's my understanding that the convolution operation is supposed to be both commutative and associative. Below is some demonstration code that produces wildly different results when convolving three functions together. If Convolve was commutative and associative, it should've given the same results for all 6 variants.

s[f_, t_, delay_, dispersion_] := Apply[f, {t - delay}]*UnitStep[t];
u[t_, delay_, dispersion_] := Exp[-t/(delay + dispersion)]*UnitStep[t];
v[t_, dispersion_] := Exp[-t/dispersion]/dispersion*UnitStep[t];

delayedAndDispersedConcentrationVariant1[f_, t_, delay_, dispersion_] :=
Module[{x, y},
Evaluate[Convolve[Evaluate[Convolve[s[f, x, delay, dispersion], u[x, delay, dispersion],
x, y, Assumptions -> x >= 0]], v[y, dispersion], y, t, Assumptions -> y >= 0]]];
delayedAndDispersedConcentrationVariant2[f_, t_, delay_, dispersion_] :=
Module[{x, y},
Evaluate[Convolve[Evaluate[Convolve[s[f, x, delay, dispersion], v[x, dispersion],
x, y, Assumptions -> x >= 0]], u[y, delay, dispersion], y, t, Assumptions -> y >= 0]]];
delayedAndDispersedConcentrationVariant3[f_, t_, delay_, dispersion_] :=
Module[{x, y},
Evaluate[Convolve[Evaluate[Convolve[u[x, delay, dispersion], s[f, x, delay, dispersion],
x, y, Assumptions -> x >= 0]], v[y, dispersion], y, t, Assumptions -> y >= 0]]];
delayedAndDispersedConcentrationVariant4[f_, t_, delay_, dispersion_] :=
Module[{x, y},
Evaluate[Convolve[Evaluate[Convolve[v[x, dispersion], s[f, x, delay, dispersion],
x, y, Assumptions -> x >= 0]], u[y, delay, dispersion], y, t, Assumptions -> y >= 0]]];
delayedAndDispersedConcentrationVariant5[f_, t_, delay_, dispersion_] :=
Module[{x, y},
Evaluate[Convolve[Evaluate[Convolve[u[x, delay, dispersion], v[x, dispersion],
x, y, Assumptions -> x >= 0]], s[f, y, delay, dispersion], y, t, Assumptions -> y >= 0]]];
delayedAndDispersedConcentrationVariant6[f_, t_, delay_, dispersion_] :=
Module[{x, y},
Evaluate[Convolve[Evaluate[Convolve[v[x, dispersion], u[x, delay, dispersion],
x, y, Assumptions -> x >= 0]], s[f, y, delay, dispersion], y, t, Assumptions -> y >= 0]]];


If we then compute the convolution for an arbitrary function (q[t] in this example) using the different variants (which differ only by commutative and associative rearrangements), the resulting plots are very different.

q[t_] := Piecewise[{{0.0, t < 1.0}, {3.0, 1.0 <= t <= 2.0}, {0.0,
2.0 < t < 3.5}, {5.0, 3.5 <= t <= 5.0}, {0.0, t > 5.0}}];

Plot[q[t], {t, 0, 20}]

Plot[delayedAndDispersedConcentrationVariant1[q, t, 3.0, 2.0], {t, 0, 20}]
Plot[delayedAndDispersedConcentrationVariant2[q, t, 3.0, 2.0], {t, 0, 20}]
Plot[delayedAndDispersedConcentrationVariant3[q, t, 3.0, 2.0], {t, 0, 20}]
Plot[delayedAndDispersedConcentrationVariant4[q, t, 3.0, 2.0], {t, 0, 20}]
Plot[delayedAndDispersedConcentrationVariant5[q, t, 3.0, 2.0], {t, 0, 20}]
Plot[delayedAndDispersedConcentrationVariant6[q, t, 3.0, 2.0], {t, 0, 20}]


Variants 1 and 2 give the expected result, a delayed and smoothed version of the input function q[t].

• I'll look into it and forward this issue as well – polfosol Jan 24 at 10:24