6
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ note that the mu in f and g should most likely be different variables, e.g. the mu in f is the mass of a resonance while g (I guess you have a resolution function in mind) would be centered at zero. $\endgroup$ Commented Feb 7, 2018 at 13:43

1 Answer 1

8
$\begingroup$

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

l = Assuming[{γ > 0 && σ > 0 && μ > 0 && 
    k ∈ Reals},
  FourierTransform[PDF[CauchyDistribution[μ, γ], x], x, k]
  ]

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

g = Assuming[{γ > 0 && σ > 0 && μ > 0 && 
    k ∈ Reals},
  Expectation[Exp[I k x], 
   x \[Distributed] NormalDistribution[μ, σ]]
  ]

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

Assuming[{γ > 0 && σ > 0 && μ > 0 && 
   k ∈ 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.

$\endgroup$
5
  • 1
    $\begingroup$ If you're considering the convolution as an integral over [-Infinity, Infinity], isn't it easier to use residues ? $\endgroup$ Commented May 28, 2013 at 6:47
  • $\begingroup$ Wow. This is great. Thank you Jens. This is a very good starting point. $\endgroup$ Commented May 28, 2013 at 7:01
  • 2
    $\begingroup$ 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! $\endgroup$ Commented May 28, 2013 at 9:38
  • $\begingroup$ @OleksandrR. Thanks for the idea. But I'm afraid I can only remember commands with fewer than 10 syllables... $\endgroup$
    – Jens
    Commented May 28, 2013 at 15:33
  • $\begingroup$ @b.gatessucks Yes - actually I thought that Mathematica would use that to do the integral in Convolve automatically... $\endgroup$
    – Jens
    Commented May 28, 2013 at 15:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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