Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Let's say I want to convolve two functions (f and g), a gaussian with a breit-wigner:

f[x_] := 1/(Sqrt[2 π] σ)Exp[-(1/2) ((x - μ)/σ)^2];
g[x_] := 1/π (γ/((x - μ)^2 + γ^2));

One way is to use Convolve like:

Convolve[f[x],g[x],x,y];

But that gives:

(γ Convolve[E^(-((x - μ)^2/(2 σ^2))),1/(γ^2 + (x - μ)^2), x, y])/(Sqrt[2] π^(3/2) σ)

,which means it couldn't do the convolution.

I then tried the integration (the definition of the convolution):

Integrate[f[x]*g[y - x], {x, 0, y}, Assumptions->{x > 0, y > 0}]

But again, it couldn't integrate. I know that there are functions that can't be integrated analytically, but it seems to me that whenever I go into convolution, I find another function that can't be integrated.

Is the numerical integration the only way to do convolution in Mathematica (besides those simple functions in the examples), or am I doing something wrong?

My target is to convolute a crystal-ball with a breit-weigner. The CB is something like:

Piecewise[{{norm*Exp[-(1/2) ((x - μ)/σ)^2], (
x - μ)/σ > -α},
{norm*(n/Abs[α])^n*
 Exp[-(1/2) α^2]*((n/Abs[α] - Abs[α]) - (
   x - μ)/σ)^-n, (x - μ)/σ <= -α}}]

I've done this in C++ but I thought I try it in Mathematica and use it to fit some data. So please tell me if I have to make a numerical integration routine in Mathematica or there's more to the analytic integration.

share|improve this question

1 Answer 1

up vote 6 down vote accepted

Here is a way to solve this problem using the convolution theorem:

l = Assuming[{\[Gamma] > 0 && \[Sigma] > 0 && \[Mu] > 0 && 
    k \[Element] Reals},
  FourierTransform[PDF[CauchyDistribution[\[Mu], \[Gamma]], x], x, k]
  ]

$\frac{\left(\theta (-k) e^{2 \gamma k}+\theta (k)\right) e^{-k (\gamma -i \mu )}}{\sqrt{2 \pi }}$

g = Assuming[{\[Gamma] > 0 && \[Sigma] > 0 && \[Mu] > 0 && 
    k \[Element] Reals},
  Expectation[Exp[I k x], 
   x \[Distributed] NormalDistribution[\[Mu], \[Sigma]]]
  ]

$e^{-\frac{k^2 \sigma ^2}{2}+i k \mu }$

Assuming[{\[Gamma] > 0 && \[Sigma] > 0 && \[Mu] > 0 && 
   k \[Element] Reals}, InverseFourierTransform[g l, k, x]]

$\frac{e^{-\frac{(-i \gamma -2 \mu +x)^2}{2 \sigma ^2}} \left(1-i \text{erfi}\left(\frac{-i \gamma -2 \mu +x}{\sqrt{2} \sigma }\right)\right)+e^{-\frac{(i \gamma -2 \mu +x)^2}{2 \sigma ^2}} \left(1+i \text{erfi}\left(\frac{i \gamma -2 \mu +x}{\sqrt{2} \sigma }\right)\right)}{2 \sqrt{2 \pi } \sigma }$

This is the result for the convolution as a function of x.

To take the Fourier transform of the Gaussian, I used Expectation with NormalDistribution (which is the Gaussian), applied to Exp[I k x]. That's faster than using FourierTransform directly - but the latter would also work.

Finally, I just take the inverse Fourier transform of the product of Fourier transforms, and fortunately we get an analytic result. It's still got some imaginary things in it, but they cancel. One could do some more work to get rid of them explicitly, but at least this should get you started.

By replacing the functions in your question with the built in distribution functions, I get slightly different prefactors. That can be amended by looking at the functional forms of the PDFs and comparing with your desired functions.

share|improve this answer
1  
If you're considering the convolution as an integral over [-Infinity, Infinity], isn't it easier to use residues ? –  b.gatessucks May 28 '13 at 6:47
    
Wow. This is great. Thank you Jens. This is a very good starting point. –  user2285967 May 28 '13 at 7:01
2  
One can also get this result in version 8 using PDF[ParameterMixtureDistribution[CauchyDistribution[µ, δ], µ \[Distributed] NormalDistribution[0, σ]], z], and VoigtDistribution is new in 9. However, using the convolution theorem seems to be more reliable (e.g. it still works if the means are both nonzero, and isn't sensitive to the order in which the distributions appear). I wouldn't have thought to try it since I wouldn't expect it to be so much better behaved, but will remember this in future. +1! –  Oleksandr R. May 28 '13 at 9:38
    
@OleksandrR. Thanks for the idea. But I'm afraid I can only remember commands with fewer than 10 syllables... –  Jens May 28 '13 at 15:33
    
@b.gatessucks Yes - actually I thought that Mathematica would use that to do the integral in Convolve automatically... –  Jens May 28 '13 at 15:35

Your Answer

 
discard

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.