Convolution of two functions (symbolic) -> closed-form expression

23 visualizzazioni (ultimi 30 giorni)
I want to solve the convolution integral of two functions (fx and fy) obtaining a closed-form expression
syms x t l1 l2 C positive
fx(x) = (l1*l2/(l1-l2))*(exp(-l2*x)-exp(-l1*x));
fy(x) = (C/sqrt(2*pi))*(1/(x^(3/2)))*exp(-C^2/(2*x));
fz = int(fx(x)*fy(t-x),'x',-Inf,+Inf);
Result:
simplify(fz) =
int(-(2383560469838099*exp(-11/(50*(t - x)))*(exp(-x)/9 - exp(-x/10)/9))/(9007199254740992*(t - x)^(3/2)), x, -Inf, Inf)
No luck using Fourier transform
simplify(ifourier(fourier(fx,x,t)*fourier(fy,x,t),t,x))
=
-(2383560469838099*fourier(fourier(exp(-11/(50*x))/x^(3/2), x, t)*fourier(exp(-x), x, t), t, -x) - 2383560469838099*fourier(fourier(exp(-11/(50*x))/x^(3/2), x, t)*fourier(exp(-x/10), x, t), t, -x))/(162129586585337856*pi)
Is this a limitation of the symbolic toolbox or there is no closed-form of this convolution integral?
  11 Commenti
Paul
Paul il 1 Ago 2022
Modificato: Paul il 1 Ago 2022
I'm unaware of any "tricks" with ifourier.
IIRC correctly, the CF of a pdf is not defined the same as with the default parameters that Matlab uses for fourier. So either adjust the parameters that Matlab uses via sympref or use the Matlab conventions for the CF definitions.
syms l1 l2 C positive
syms x t real
fx(x) = (l1*l2/(l1-l2))*(exp(-l2*x)-exp(-l1*x))*heaviside(x);
Matlab default fourier:
Fx(t) = simplify(fourier(fx(x),x,t),100)
Fx(t) = 
Using the expressions above
Fx1(t) = (1/(1+t^2*(1/l1)^2))+1i*((t/l1)/(1+t^2*(1/l1)^2));
Fx2(t) = (1/(1+t^2*(1/l2)^2))+1i*((t/l2)/(1+t^2*(1/l2)^2));
[num,den] = numden(simplify(Fx1(t)*Fx2(t),100));
Fxp(t) = simplify(num/den)
Fxp(t) = 
These are the conjugates of each other
simplify(Fx(t) - conj(Fxp(t)))
ans = 
0
I'll be curious to see the responses at stackexchange. However, I don't think that post is correct. If the upper limit of integration is inf, the definitions of p1 and p2 each must be multiplied by heaviside(x) (actually, with the lower limit of integration being 0, I guess it's really only needed for p2, but I'd show it for both). Or make x the upper limit of integration.
You can also try Wolfram.
If a closed form expression is obtained, would you mind posting the solution back here?
j l
j l il 1 Ago 2022
I edited the post to include the current integral expression as you suggested.
Of course I will post here any solution. I tried integral-online and Wolfram without any success (the latter stated "standard computational time exceeded"...)

Accedi per commentare.

Risposta accettata

j l
j l il 29 Ago 2022
Modificato: j l il 30 Ago 2022
The closed-form solution of the convolution integral was obtained using Wolfram Mathematica (see https://math.stackexchange.com/questions/4501888/closed-form-convolution-integral-of-two-pdfs-hypoexponential-and-l%c3%a9vy/4520901#4520901)
assuming x>=0, l1>0, l2>0, C>0 and integrating between 0 and x (as suggested by @Paul)
The resulting function is complex:
  14 Commenti
David Goodmanson
David Goodmanson il 4 Set 2022
Hi jl,
I used erfw because of its better numercal properties, and I guess the expression in terms of erfc was unstated. In your expression above, the factor in front of the first erfc term should be
exp(-L2*x)*exp(i*C*sqrt(2*L2))
Hopefully you can go back and verify that. Similarly with L1 for the second erfc term, and the whole works is
g1 = (L1*L2/(L1-L2))* ...
( exp(-L2*x).*real(exp(i*C*sqrt(2*L2))*erfc1(C./sqrt(2*x)+i*sqrt(L2*x))) ...
-exp(-L1*x).*real(exp(i*C*sqrt(2*L1))*erfc1(C./sqrt(2*x)+i*sqrt(L1*x))) );
Since Matlab does not provide erfc with complex double precision variables, erfc1 is used:
function y = erfc1(z)
% erfc with complex argument, using erfw
%
% y = erfc1(z)
y = exp(-z.^2).*erfw(i*z);
But don't expect this version to be good for as large values of x as the original erfw version. That's because this version contains the factors like exp(-L2*x), which cancel out the large exponential growth of the erfc part. For large x, erfc will overflow before exp(-L2*x) has a chance to cancel it out. The original version cancels things out algebraically and the exponential growth does not occur.
j l
j l il 6 Set 2022
Thanks @David Goodmanson that was the initial expression I was looking for.
Bests

Accedi per commentare.

Più risposte (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by