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

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?

[David Goodmanson]

Hi Walter/JL

this expression is complex for x>t (likely not intended), and increases without bound for x--> -inf, so there need to be restrictions on fx,fy (such as both being zero for negative argument) or this integral will not work.

[Walter Roberson]

Is there good reason to expect that a closed form expression of the convolution exists?

[Walter Roberson]

As @David Goodmanson indicates, this goes complex,

syms x t 
rng(12345); l1 = rand(), l2 = rand(), C = rand()

l1 = 0.9296

l2 = 0.3164

C = 0.1839

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));

X = linspace(-10, 10, 250);

FZ =(sum(fx(X).*fy(t-X)));
vpa(FZ, 10)

ans = 

FZN = matlabFunction(FZ);
T = linspace(0,15,300);
y = FZN(T);
yr = real(y); yr = min(yr, 100);
yi = imag(y); yi = min(yi, 100);
plot(T, yr, 'k-'); title('real part')

plot(T, yi, 'r-.'); title('imaginary part')

[j l]

Dear Walter and David, thank you for your replies. That's perfectly true, since both function should be 0 for x<0. They are probability density functions defined for positive values of x and all parameters.

I though it was sufficient to set the assumption positive to variables and parameters.

Howevere, I also tried the following based on David's suggestions (notice that I substituted C with C^(1/2) just to simplify the output)

syms x t real

syms l1 l2 C P fx(x) fy(x) fz(x,t) positive

fx(x) = piecewise(x < 0,0,x >= 0,(l1*l2/(l1-l2))*(exp(-l2*x)-exp(-l1*x)) );

fy(x) = piecewise(x < 0,0,x >= 0,(sqrt(C/(2*pi)))*(1/(x^(3/2)))*exp(-C/(2*x)) );

fz = int(fx(t)*fy(x-t),'t',0,+Inf);

simplify(fz)

Result:

piecewise(0 < x, -(2^(1/2)*C^(1/2)*l1*l2*int((exp(C/(2*(t - x)))*(exp(-l1*t) - exp(-l2*t)))/(x - t)^(3/2), t, 0, x))/(2*pi^(1/2)*(l1 - l2)), x <= 0, 0)

I also tried with heaviside and sympref('HeavisideAtOrigin',1), and assuming (x>=t), with no luck.

[Torsten]

And how to solve the problem that fy is not in L^1(0,oo) because of the singularity at 0 and that fy becomes complex for x>t ?

[Walter Roberson]

(x-t)^(3/2) is a problem when x-t is negative.

[Paul]

Based on the statement by the @j l that both fx and fy zero for x<0, we have

syms l1 l2 C positive
syms x t real
fx(x) = (l1*l2/(l1-l2))*(exp(-l2*x)-exp(-l1*x))*heaviside(x);
fy(x) = (C/sqrt(2*sym(pi)))*(1/(x^(3/2)))*exp(-C^2/(2*x))*heaviside(x); % note sym(pi)

fy seems well behaved at the origin, for example

fplot(subs(fy(x),C,1),[-1 5])

limit(subs(fy(x),C,1),x,0,'right')

ans = 

0

The convolution integral is

fz(x) = int(fx(t)*fy(x-t),t,-inf,inf,'Hold',true)

fz(x) = 

We see that if t<0 the integrand is 0 and that if x-t < 0 the integrand is 0. So we don't really need worry about (x-t)^3/2 beocoming complex, and the integral is better expressed by changing the limits of integration

fz(x) = int(fx(t)*fy(x-t),t,0,x,'Hold',true)

fz(x) = 

Of course, we still don't have a closed from expression (at least I couldn't get there)

release(fz(x))

ans = 

Or

fz(x) = int(fy(t)*fx(x-t),t,0,x)

fz(x) = 

But a numerical solution seems like it should be feasible.

[David Goodmanson]

fy is in L^1(0,oo) since as t--> 0, e^(-C^2/(2*x)) gets small a lot faster than (x^(-3/2) gets large. It seems highly unlikely that this integral has a simple analytic solution.

[j l]

Thanks all.

Yes @Paul the fy(x) is the Lévy distribution probability density function thus it has heavy tail but still with finite integral.

I obtained the numerical integration with defined parameters to have an idea of the shape, but I would like to get a closed form since it is part of a biophysical model I developed.

I also posted a question about this at https://math.stackexchange.com/questions/4501888/closed-form-convolution-integral-of-two-pdfs-hypoexponential-and-l%c3%a9vy

I am also trying to solve it with the characteristic functions of the two distributions (Fourier transform of the pdf), which are available in closed form. fx(x) is the pdf of the hypoexponential distribution with 2 parameters l1, l2, i.e. the distribution of the sum of 2 exponential distributions with parameter l1 and l2, respectively.

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));

Fy(t) = exp(-sqrt(-2*1i*C*t));

The convolution is the iFourier of the product

Fz(t) = Fx1*Fx2*Fy;

fz(x) = ifourier(Fz,t,x);

simplify(fz(x))

= (l1^2*l2^2*fourier(exp(-2^(1/2)*C^(1/2)*(-t*1i)^(1/2))/(l1^2*l2^2 + l1^2*t^2 + l2^2*t^2 + t^4), t, -x) - l1*l2*fourier((t^2*exp(-2^(1/2)*C^(1/2)*(-t*1i)^(1/2)))/(l1^2*l2^2 + l1^2*t^2 + l2^2*t^2 + t^4), t, -x) + l1*l2^2*fourier((t*exp(-2^(1/2)*C^(1/2)*(-t*1i)^(1/2)))/(l1^2*l2^2 + l1^2*t^2 + l2^2*t^2 + t^4), t, -x)*1i + l1^2*l2*fourier((t*exp(-2^(1/2)*C^(1/2)*(-t*1i)^(1/2)))/(l1^2*l2^2 + l1^2*t^2 + l2^2*t^2 + t^4), t, -x)*1i)/(2*pi)

No closed-form...

Is there maybe a way to get the ifourier by imposing some restraint (e.g. integral of fz(x) = 1 or fz(x<0) = 0?

[Paul]

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]

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

Answer 1

[j l]

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:

[Paul]

@j l

Does the result actually evaluate as complex? Matlab sometimes forms expressions that look complex but really aren't, presumabably because the expressions with "i" have a simpler form. Maybe Mathematica works the same. If it is complex, does that make physical sense? I thought this problem was about the convolution of pdfs to form a new pdf? If so, does a pdf have to be real valued? Also, can the closed form solution be verified against a numerical solution?

[j l]

@Paul I agree, the output must be a probability, thus real valued between 0 and 1, but I was not able to further simplify the expression without complex expressions or special functions up to now. As for the evaluation, using meaningful parameters leads to real output <1 and I verified coincidence with the numerical convolution over a broad range of x. Still I would prefer a better solution, also because evaluation of erfc with large x (>10^3) in those complex arguments produces NaN with all the matlab erfc implementations I tried and with Mathematica as well.

[Paul]

The output is a probability density function, not a probability, and so doesn't have to be between 0 and 1, though of course it could be in this particular problem.

How are you evaluating those erfc functions with complex inputs in Matlab? Can you show an example where erfc returns NaN?

[David Goodmanson]

The function is going to be real. Taking the first part involving i2, and ignoring the real factor e^(i2*x) in front, it's of the form

e^(-ia) ( e^2ia (erfc(b+ic) + erfc(b-ic) )            (1)
=  e^ia erfc(b+ic) + e^(-ia) erfc(b-ic)

Since erfc obeys the Schwarz reflection principle

f(x*) = f(x)*

this becomes

e^ia erfc(b+ic) + [ e^(ia) erfc(b+ic) ]*
= 2 Re( e^ia erfc(b+ic) )                              (2)

so it's real. This means that when you evaluate it you don't have to compute erfc twice, as in (1). You can just go with (2). Same for the second part involving i1.

[j l]

I edited my answer to include the solution from the original integral with C, which is slighly more compact and the subsitution was only aimed at using the standard Levy pdf form.

Thanks for the replies.

@Paul I used many implementations for erfc including double(erfc(sym(z))) which is very slow but produced results for a wider range, faddeeva package https://it.mathworks.com/matlabcentral/fileexchange/38787-faddeeva-package-complex-error-functions and almost all the other functions found.

In the following code aa is the first erfc argument in the final pdf:

clc
erfcc = @(x) double(erfc(sym(x)));

aa = @(x,l1,l2,c) 2.^(-1/2).*C.*x.^(-1/2)+(sqrt(-1)*(-1)).*(l1.*x).^(1/2);

b1 = aa(713,1,0.1,0.6633)
erfc_b1 = erfcc(b1)
b2 = aa(714,1,0.1,0.6633)
erfc_b2 = erfcc(b2)
b3 = aa(715,1,0.1,0.6633)
erfc_b3 = erfcc(b3)


b1 =
   0.0117 -26.7021i

erfc_b1 =
 -5.5258e+307 +7.7106e+307i

b2 =
   0.0116 -26.7208i

erfc_b2 =
 -1.5010e+308 +        Infi

b3 =
   0.0116 -26.7395i

erfc_b3 =
     -Inf +    Infi

Which results in Inf values, while other implementations give NaN.

However I verified that Mathematica can actually compute the erfc for all values but there are other problems related to evaluation for larger x with e^-x.

[j l]

As for @David Goodmanson suggestion, thanks, but if I understood well, that approach would be useful for faster evaluation but cannot help simplifying the expression, right? Actually faddeeva package is super fast compared to symbolic toolbox... still it would be ideal to get a non-complex version of the expression or at least a more compact one.

[Paul]

I wonder why erfc doesn't accept complex inputs?

Isn't aa the last erfc argument?

The reason that erfc_b3 is inf is because erfc(sym(b3)) has both real and imaginary part greater than realmax.

However, the actual value of fz(715) with L1, L2, and C defined here is

0.000014167625983587929920562835962237 ,

give or take a digit or two, as determined by implementing the closed form expression symbolically, subs-ing in the values of the constants, and taking the real part of the vpa of the result. Got the same answer using vpa on the integral expression of fz(x) expressed in terms of int, which I guess actually uses vpaintegral.

How far out into the tail of fz(x) is needed?

[David Goodmanson]

Hello jl,

I'm not sure whether you object more to the compex variable calculation or the fact that Matlab is slow at it. As to the speed, I think it's highly to Mathwork's discredit that they do not have a well organized suite of special functions that work in double precision and take complex input.

In the example that Paul mentioned, why shouldn't there be a fast double precision calculation of erf(x+iy)? In some cases Mathworks shuttles you off to syms, seemingly as an afterthought, where you can get an answer, but of course it's slow.

Fortunately there are contributed files to fill in some of the gaps where Mathworks is lacking. The code below leaves out the (l1 12)/(l1-l2) factor in front and does a calculation with just one of the exponentials, exp(-Lx).

The code uses the erfw function, where

erfc(z) = e^(-z^2) erfw(iz).

It appears to work reliably for input x as least as large as 10^13. (That's for C=1, L = 1. I didn't try any other values). I'm not sure what you consider to be fast, but on my PC I get 10 million values in about 0.8 sec.

% convolution of
% f1(x) = exp(-L*x)
% f2(x) = (C/sqrt(2*pi))*(1/(x^(3/2)))*exp(-C^2/(2*x))

C = 1;
L = 1;
x = 0:.001:1000;
% x = 0:.0001:1000;                % 1e7 values

A = (1/2)*exp(i*C*sqrt(2*L));
a = C./(sqrt(2*x));
b = sqrt(L*x);
E = exp(-(a.^2 + 2*i*a.*b));
f3 = 2*real(A*E.*erfw(i*(a+i*b)));
f3(x==0) = 0;

figure(2)
plot(x,f3)
grid on
xlim([0 40])

function w = erfw(z,n)

% computes w(z)=exp(-z^2)*erfc(-i*z) using a rational series
% with n terms.  n is optional, default is 32.  n = 16 gives
% about six significant digits and n = 32 about 12 significant
% digits (at least in the upper half plane).
% w = erfw(z,n)

% calculation is direct for Im(z)>=0.
% For Im(z)<0, use w(z) = 2*exp(-z^2) - w(-z).
% From Computation of the Complex Error Function, J.A.C Weideman,
% SIAM J. Numerical Analysis, 31, 5, pp 1497-1518.
% provided by Doug DeWolf 

% Although erfw is an integral over the real line, it does
% NOT have a branch cut, because a correction is made when Im(z) < 0  (dg)

if nargin==1; n = 32; end
m = 2*n;
k = (-m+1:m-1)';
L = sqrt(n/sqrt(2));
theta = k*pi/m;
t = L*tan(theta/2);
f = exp(-t.^2).*(L^2 + t.^2);
f = [0; f];
a = real(fft(fftshift(f)))/(2*m);
a = flipud(a(2:n+1));
iminus = (imag(z)<0);
z(iminus) = -z(iminus);
Z = (L+i*z)./(L-i*z);
p = polyval(a,Z);
w = 2*p./(L-i*z).^2 + (1/sqrt(pi))./(L-i*z);
w(iminus) = 2*exp(-z(iminus).^2) - w(iminus);
end

[j l]

@David Goodmanson as I commented above, with external packages such as faddeeva the computation of erfc is fast so speed is not a problem, as those functions do not require symbolic toolbox. The problem is that with the mentioned arguments erfc returns -Inf, probably related to an overflow as suggested by @Paul: for instance, the output of the first erfc for x = 714, l1 = 1, l2 = 0.1, c = 0.44, computed with Mathematica is -2.076642385572023*10^308 + 1.524858881404226*10^308 i likely above realmax.

I tried your implementation with erfw but the problem is the same (for that argument).

I am trying the approach by @Paul with vpa but it is very slow for many x values and it is still computing...

Just to clarify, evaluating the closed-form expression directly was needed only for veryfing its validity compared to the numerical convolution, while it is ok to use the latter for any other purpose... Still it would be nice to solve this computational issue and get at least x = 10^4 to show a broad plot of the heavy tail.

As for your previous comment @David Goodmanson, are you able to manually simplify the closed-form expression based on Schwarz reflection principle?

Thanks for your replies

[j l]

Ok actually the solution by @Paul was very effective although a bit slow:

syms f0s(x)
C = 0.6633
l1 = 1;
l2 = 0.1;
f0s(x) = (1/2).*l1.*(l1+(-1).*l2).^(-1).*l2.*((-1).*exp(1).^((sqrt( ...
  -1)*(-1)).*2.^(1/2).*C.*l1.^(1/2)+(-1).*l1.*x).*(erfc(2.^( ...
  -1/2).*C.*x.^(-1/2)+(sqrt(-1)*(-1)).*(l1.*x).^(1/2))+exp(1) ...
  .^((sqrt(-1)*2).*2.^(1/2).*C.*l1.^(1/2)).*erfc(2.^(-1/2).* ...
  C.*x.^(-1/2)+sqrt(-1).*(l1.*x).^(1/2)))+exp(1).^((sqrt(-1)*( ...
  -1)).*2.^(1/2).*C.*l2.^(1/2)+(-1).*l2.*x).*(erfc(2.^(-1/2).* ...
  C.*x.^(-1/2)+(sqrt(-1)*(-1)).*(l2.*x).^(1/2))+exp(1).^(( ...
  sqrt(-1)*2).*2.^(1/2).*C.*l2.^(1/2)).*erfc(2.^(-1/2).*C.* ...
  x.^(-1/2)+sqrt(-1).*(l2.*x).^(1/2))));

f0vpa = @(x) real(vpa(f0s(x)));

works accurately for x = logspace(-5,5,10^4) and takes 207.73 secs with my mac M1.

[David Goodmanson]

Hi jl,

I forgot to put L = 1 C=1 at the top of the code so it will run (fixed now).

Aside from that, I don't know what you mean, since the code I suggested does not use the symbolic toolbox, uses the Schwarz reflection principle, and works fine for x up to about 10^13 (with L=1, C=1).

The full calculation for l1 and l2 is included below. It uses a maximum x value of 10^10, uses 10^6 points and on my PC takes about 0.2 seconds.

% convolution of
% f1(x) = (L1L2/(2(L1-L2)) (exp(-L2*x) -exp(-L1*x)) 
% f2(x) = (C/sqrt(2*pi))*(1/(x^(3/2)))*exp(-C^2/(2*x))
C = 0.6633
L1 = 1;
L2 = 0.1
x = logspace(-5,10,1e6);
g1 = convf1f2(x,C,L1);
g2 = convf1f2(x,C,L2);
g = (L1*L2/(2*(L1-L2)))*(g2-g1);
figure(2)
loglog(x,g)
grid on

function g = convf1f2(x,C,L)
A = exp(i*C*sqrt(2*L));
a = C./(sqrt(2*x));
b = sqrt(L*x);
E = exp(-(a.^2 + 2*i*a.*b));
g = 2*real(A*E.*erfw(i*(a+i*b)));     % schwarz
g(x==0) = 0;
end

function w = erfw(z,n)
     (same as before, posted already)
end

[j l]

Thanks @David Goodmanson your code works like a charm and is super fast and accurate.

I then tried to derive the full expression using the Schwarz reflection principle and using erfc (not erfw) but I guess there is some error, could you check it?

[David Goodmanson]

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]

Thanks @David Goodmanson that was the initial expression I was looking for.

Bests