Convolve does not get the correct answer

Question

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

is obviously false, but why? Any suggestions ?

Answer 1

$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 *)

"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.

Answer 2

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.