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.
Answer
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.
Comments
Post a Comment