Convolve[Sinc[x], Exp[-x^2], x, X]
(* E^-X^2 π *)

is obviously false, but why? Any suggestions ?

  • 1
    My Mathematica returns this unevaluated. 10.0.2 on OS X 10.10.3 – happy fish Mar 15 '15 at 15:01
  • Mine too (10.0.2.0 on Win 8.1 x64). – Jinxed Mar 15 '15 at 15:06
  • 1
    Consider using different variable names as this brings confusion. Convolve[Sinc[v], Exp[-w^2], v, w] evaluates to E^-w^2 [Pi] @ Mathematica 9.0.0.1 and Mathematica 10 quits the kernel – Sektor Mar 15 '15 at 15:09
  • Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – user9660 Mar 15 '15 at 15:10
  • 1
    @Sektor Are you sureConvolve[Sinc[v], Exp[-w^2], v, w] and Convolve[Sinc[v], Exp[-v^2], v, w] are the same thing? – happy fish Mar 15 '15 at 15:19
up vote 16 down vote accepted
$Version

(*
"8.0 for Microsoft Windows (64-bit) (October 7, 2011)"
*)

Direct attack fails:

Timing[Convolve[Sinc[x], Exp[-x^2], x, y]]

(*
Out[218]= {59.296, Convolve[Sinc[x], E^-x^2, x, y]}
*)

or, equivalently,

Timing[Integrate[Sinc[x] Exp[-(x - y)^2], {x, -∞, ∞}] ]

$\left\{49.92,\int_{-\infty }^{\infty } e^{-(x-y)^2} \text{Sinc}[x] \, dx\right\}$

Let us now solve the problem, using Fourier transform.

We have

FourierTransform[Sinc[x], x, t]

(*
Out[206]= 1/2 Sqrt[π/2] (Sign[1 - t] + Sign[1 + t])
*)

Therefore we can write

InverseFourierTransform[1/2 Sqrt[π/2] (Sign[1 - t] + Sign[1 + t]), t, x]

(*
Out[207]= Sin[x]/x
*)

or, explicitly

Integrate[Exp[-I t x] (1/
    4 (Sign[1 - t] + Sign[1 + t])), {t, -∞, ∞}]

(*
Out[211]= Sin[x]/x
*)

Therefore doing the x integration first we have

Integrate[Exp[-I t x] (1/
    4 (Sign[1 - t] + Sign[1 + t])) Exp[-(x - 
      y)^2], {x, -∞, ∞}]

(*
Out[212]= 1/4 E^(-(1/4) t (t + 4 I y)) Sqrt[π] (Sign[1 - t] + Sign[1 + t])
*)

and the t integration finalizes the calculation giving for the convolution the following result:

Integrate[1/4 E^(-(1/4) t (t + 4 I y))
   Sqrt[π] (Sign[1 - t] + Sign[1 + t]), {t, -∞, ∞}]

(*
Out[215]= 1/2 E^-y^2 π (Erf[1/2 - I y] + Erf[1/2 + I y])
*)

EDIT #1

17.03.15 Comparision and analysis

Let us first compare the answers of Wolfgang and Jens

fWolfgang[y_] := 1/2 E^-y^2 π (Erf[1/2 - I y] + Erf[1/2 + I y])

fJens[X_] := -(1/2) E^-X^2 Pi Erfc[1/2 - I X] - 1/2 E^-X^2 Pi Erfc[1/2 + I X]

Because $erfc(z) = 1 - erf(z)$ this can be written

f1Jens[y_] := -(1/2) E^-y^2 π (1 - Erf[1/2 - I y] + 1 - Erf[1/2 + I y])

The difference is just

fWolfgang[x] - f1Jens[x] // Simplify

(*
Out[9]= E^-x^2 π
*)

Now, in order to see the "true" result, let's define the numeric integral

fNum[y_] := 
 NIntegrate[Sin[x]/x Exp[-(y - x)^2], {x, -∞, ∞}]

Comparing the results graphically (ignoring the error messages of the integration) gives

Plot[{fWolfgang[y] + 0.1, fNum[y]}, {y, -6, 6}]
(* 150317_Plot _fW _fN *)

enter image description here

"proves" that fWolfgang is correct.

Do you want still another incorrect result from correct input? Here we go:

Writing Sin[x] = 1/(2 I) (Exp[I x] - Exp[-I x]) our integral becomes

fSplit[y_] = 
 Integrate[(Exp[I x] - Exp[-I x])/(2 I x)
    Exp[-(y - x)^2], {x, -∞, ∞}]

(*
Out[16]= 1/2 I E^-y^2 (Log[-I - 2 y] - Log[I - 2 y] + Log[-I + 2 y] - Log[I + 2 y])
*)

Here there even is no error function. And the result is obviously wrong:

fSplit[0.]

(*
Out[19]= 3.14159 + 0. I
*)

fWolfgang[0.]

(*
Out[20]= 1.6352 + 0. I
*)

The same (wrong) result is obtained usind the option PrincipalValue->True in order to tell Mathematica how to deal with the false pole at x = 0.

But let's look at the ostensible pole in more detail. This integral is obviously divergent at x = 0:

Integrate[Exp[I x]/(2 I x) Exp[-(y - x)^2], {x, -∞, ∞}]

During evaluation of In[21]:= Integrate::idiv: Integral of E^(I x-(-x+y)^2)/x does not converge on {-∞,∞}. >>

$\int_{-\infty }^{\infty } -\frac{i e^{i x-(-x+y)^2}}{2 x} \, dx$

But taking the pricipal value the result is finite

Integrate[Exp[I x]/(2 I x) Exp[-(y - x)^2], {x, -∞, ∞}, 
 PrincipalValue -> True]

(*
Out[22]= 1/2 I E^-y^2 (Log[-I - 2 y] - Log[I + 2 y])
*)

and it is part of the wrong result fSplit[] above.

The integral can also be written as a fourier transform

Sqrt[2 π] FourierTransform[1/(2 I x) Exp[-(y - x)^2], x, t] /. t -> 1

(*
Out[27]= -(1/2) I E^-y^2 (-Log[-I - 2 y] + Log[I + 2 y])
*)

But it still leads to the same wrong result.

Summarizing we find that splitting the Sin[] into a sum of complex exponentials the resulting integral leads rather consistently to a wrong result.

Concluding

(i) it is not only Convolve which produces wrong results but also related integrals do.

(ii) I realize that I was just lucky having found the correct result by attacking the problem using Fourier transformation.

EDIT #2

In order to mitigate the pessimistic outlook here's a positive message:

We can replace the lengthy Fouriertransform approach by this one.

A pole 1/x can be produced by an auxiliary integration.

Indeed, we can write

Integrate[Cos[t x], {t, 0, 1}]

(* Out[45]= Sin[x]/x *)

Changing the order of integration, doing the x-integral first, we get

Integrate[Cos[t x] Exp[-(x - y)^2], {x, -∞, ∞}]

(* 
Out[46]= 1/2 E^(-(1/4) t (t + 4 I y)) (1 + E^(2 I t y)) Sqrt[π] 
*)

and doing the t-integral subsequently we have

Integrate[%, {t, 0, 1}]

(*
Out[47]= 1/2 E^-y^2 π (Erf[1/2 - I y] + Erf[1/2 + I y])
*)

which is the correct result.

Check:

% /. y -> 1.

(* Out[44]= 1.39248 + 0. I *)

We can also "save" convole.

We have to consider (before the t-integration)

Convolve[Cos[t x], Exp[-x^2], x, y]

(*
Out[57]= 1/2 E^(-(1/4) t (t + 4 I y)) (1 + E^(2 I t y)) Sqrt[π]
*)

and the t-integral

Integrate[%, {t, 0, 1}]

(*
Out[55]= 1/2 E^-y^2 π (Erf[1/2 - I y] + Erf[1/2 + I y])
*)

% /. y -> 1.

(*
Out[56]= 1.39248 + 0. I
*)

gives the correct result.

EDIT #3

I have found a transparent way to generate the result of Jens. This show where the problem lies.

Let us again consider the integral

h = 1/(2 I )
   Integrate[Exp[I x]/x Exp[(x - y)^2], {x, -∞, ∞}]

and let us shift the integration variable thus

Exp[I x]/x Exp[-(x - y)^2] /. x -> u + y

(*
Out[66]= E^(-u^2 + I (u + y))/(u + y)
*)

completing the square

Expand[-(u - I/2)^2]

(*
Out[67]= 1/4 + I u - u^2
*)

we can write

$\text{Exp}[i y-1/4]\int_{-\infty }^{\infty } \frac{e^{-(u-i/2)^2}}{y+u} \, du$

Now shifting again, this time into the complex plane

E^-(u - I/2)^2/(y + u) /. u -> v + I/2

(*
Out[68]= E^-v^2/(I/2 + v + y)
*)

giving

$\text{Exp}[i y-1/4]\int_{-\infty -i/2}^{\infty -i/2} \frac{e^{-v^2}}{y+v+i/2} \, dv$

Now the trick from the good old university days when calculating the Fourier transform of Exp[-x^2]: we shift the integration path in the u-plane which lies 1/2 unitites below the real axis and paralell to it, to the real axis. This gives

$\frac{1}{2i}\text{Exp}[i y-1/4]\int_{-\infty }^{\infty } \frac{e^{-w^2}}{y+w+i/2} \, dw$

(*
Out[69]= ConditionalExpression[-(1/2) I E^(-(1/4) + I y - 
   1/4 (I + 2 y)^2) (I π Erf[1/2 - I y] + Log[-I - 2 y] - Log[I + 2 y]), 
 Im[y] != -(1/2)]
*)

Taking the input format we can add the condition that y>0 (for ins

(1/(2*I))*Exp[I*y - 1/4]*
  Integrate[1/(E^w^2*(y + w + I/2)), {w, -Infinity, Infinity}, 
   Assumptions -> y ∈ Reals] // Simplify

(*
Out[71]= -(1/2) E^-y^2 π Erfc[1/2 - I y]
*)

For the complete (sinc) integral we need to add the complex conjugate h* of it, giving

fShift[y_] = -(1/2) E^-y^2 π Erfc[1/2 - I y] + -(1/2) E^-y^2 π Erfc[
    1/2 + I y]

(*
Out[72]= -(1/2) E^-y^2 π Erfc[1/2 - I y] - 1/2 E^-y^2 π Erfc[1/2 + I y]
*)
fShift[1.]

(* 
Out[73]= 0.236748 + 0. I 
*)

We can easily verify that this function is exactly the result fJens[] of Jens.

This means, however, that the clue lies in the shifting of the integration path. This shifting blurres the required exact treatment of the pole.

EDIT #4: Miscellaneous results

1) Proof by series expansion

Just to fill a small gap: in the "proof" of correctness of fWolfgang we resorted to numerical integration.

Now we shall do it by expansion into power series

fWolfgang[y]

(*
Out[142]= 1/2 E^-y^2 π (Erf[1/2 - I y] + Erf[1/2 + I y])
*)

Series[fWolfgang[y], {y, 0, 6}] // Normal

(*
Out[154]= π Erf[1/2] + y^2 (Sqrt[π]/E^(1/4) - π Erf[1/2]) + 
 y^6 ((71 Sqrt[π])/(360 E^(1/4)) - 1/6 π Erf[1/2]) + 
 y^4 (-((7 Sqrt[π])/(12 E^(1/4))) + 1/2 π Erf[1/2])
*)

Expanding the expression Exp[-(x-y)^2] in the integrand with respect to y, and integrating term by term gives up to the order y^6:

Collect[Integrate[
   Sin[x]/x Series[Exp[-(x - y)^2], {y, 0, 6}] // 
    Normal, {x, -∞, ∞}] // Expand, y]

(*
Out[163]= π Erf[1/2] + y^2 (Sqrt[π]/E^(1/4) - π Erf[1/2]) + 
 y^6 ((71 Sqrt[π])/(360 E^(1/4)) - 1/6 π Erf[1/2]) + 
 y^4 (-((7 Sqrt[π])/(12 E^(1/4))) + 1/2 π Erf[1/2])
*)

which agrees with the expansion of fWolfgang.

This is not a strict proof, of course, as we have considered only a finite number of terms. But I promise to the first one who finds a term which does not agree a bottle of fine German beer.

2) The innocent "pole"

The "pole" at x = 0 alone is not the cause of trouble.

Look at this example where I have replaced the Gaussian by a Cauchy weight

Convolve[(Sin[x]/x), 1/(1 + x^2), x, y, Assumptions -> y > 0]

(*
Out[197]= (π (E - Cos[y] + y Sin[y]))/(E (1 + y^2))
*)

or, in explicit form,

Integrate[(Sin[x]/x) 1/(1 + (x - y)^2), {x, -∞, ∞}, 
 Assumptions -> y > 0]

(*
Out[200]= (π (E - Cos[y] + y Sin[y]))/(E (1 + y^2))
*)

Both operations are performed by Mathematica without problems.

I conclude that it is the combination of the "pole" and the esssential singularity of Exp[-x^2] at infinity which gives rise to the observed difficulties.

  • You are right , another answer could be : In[84]:= InverseFourierTransform[ Sqrt[Pi/2] UnitBox[t/2], t, x] Out[84]= Sinc[x] In[85]:= InverseFourierTransform[ 1/Sqrt[2] Exp[-(t/2)^2], t, x] Out[85]= E^-x^2 In[86]:= FourierTransform[UnitBox[t/2] Sqrt[Pi/2] 1/Sqrt[2] Exp[-(t/2)^2], t, x] Out[86]= 1/2 E^-x^2 Sqrt[[Pi]/2] (Erf[1/2-I x]+Erf[1/2+I x]) But, is there an option in Convolve to get the right answer directly ? – Hergé Mar 15 '15 at 17:34
  • 2
    I always enjoy your detailed analyses and answers. +1 – ciao Mar 15 '15 at 20:38
  • Your answer is probably the only one that works because Convolve in my answer is buggy. In that case, I think the question is actually a duplicate of Convolving/integrating problems. – Jens Mar 16 '15 at 18:56
  • 1
    @rasher: thank you very much for your kind remark which made my day :-) As a reward I have taken the liberty to continue the analysis a bit (see my EDIT #1) – Dr. Wolfgang Hintze Mar 17 '15 at 11:43
  • 1
    I like that change-of-order (edit#2). Cannot give a second upvote though. – Daniel Lichtblau Mar 17 '15 at 14:14

Here is another way that allows you to directly use Convolve:

Convolve[TrigToExp@FunctionExpand[Sinc[x]], Exp[-x^2], x, X]

(*
==> -(1/2) E^-X^2 Pi Erfc[1/2 - I X] - 
 1/2 E^-X^2 Pi Erfc[1/2 + I X]
*)

In order to get a successful evaluation, I just had to break up the Sinc function into its complex exponential terms.

  • Nice answer ! Thank you. – Hergé Mar 15 '15 at 18:21
  • Well. Humm I am sorry for my enthusiast answer yesterday, but Convolve doesnt work again. If i plot the answer (Jens) and the answer (Dr Hintze),Jens result is wrong. Convolve is the problem, because TrigToExp@FunctionExpand[Sinc[x]] works. – Hergé Mar 16 '15 at 18:12
  • Right, this seems to be a huge bug in Convolve. – Jens Mar 16 '15 at 18:43
  • I did some testing and sent along in house what I found. Agreed, it's a bug. – Daniel Lichtblau Mar 16 '15 at 21:07
  • @DanielLichtblau Thanks for taking a look. It's sneaky because the c in Erfc is easy to overlook, and with Erf it would almost be right (up to a sign). – Jens Mar 16 '15 at 23:02

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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