Why is the convolution calculated wrong? I need to do interp1 and multiply by dx,

Hi, I have a PDE that to run in time, I need to at each time step calculate the convolution of two functions.

First: I notice that we need to multiply by the dx, the discretization of space. Can anyone explain this to me?

Second: I find that convolution needs to be re-interpolated in some cases, and I wonder why this need to happen to be more accurate?

Third: The following Code performs the convolution of two gaussians in different ways but for N=odd, that is (N+1 even) number of points, we have less accurate. Thanks in advance

FIG=240;
N=1024,   L=30;
LxU = L; LxL =(-L); 
dx=(LxU-LxL)/N;  x=linspace(LxL,LxU,N+1);

% Convolution of two gaussian is another gaussian

mu_1=1; mu_2=2; sigma_1=2; sigma_2=1;
f=exp(-(x-mu_1).^2/(2*sigma_1^2))/(sigma_1*sqrt(2*pi));
g=exp(-(x-mu_2).^2/(2*sigma_2^2))/(sigma_2*sqrt(2*pi));

g2=exp(-(-(x)-mu_2).^2/(2*sigma_2^2))/(sigma_2*sqrt(2*pi));
% Convolution answer
fg=1/( sqrt( 2*pi*(sigma_1^2+sigma_2^2) ) )*...
      exp( -[x-(mu_1+mu_2)].^2 / [2*(sigma_1^2+sigma_2^2)] );
% By hand
FG1=[];
for h=1:length(x);
  gh=exp(-((-1*(x-x(h)))-mu_2).^2/(2*sigma_2^2))/(sigma_2*sqrt(2*pi));
  FG1(h)=trapz(x,f.*gh);    
end

 % Conv: All need the dx mult. Why?
 FG2  =dx*conv2(f,g,'same');
 FG2p =dx*conv(f,g,'same');
 FG2pp=dx*cconv(f,g,N+1);

 % fft Conv Theorem: 
 % needs to multiply by dx, AND fftshift
 ft  = fftshift(fft( f ));    gt  = fftshift(fft( g ));
 FG3 = ifft(ifftshift( ft.*gt ));
 FG3 = dx*fftshift(FG3);

% Circular conv need to be FFTShifted FG2pp=fftshift(FG2pp);

 if mod(N,2)==1
     FG2b  =interp1(x+dx/2,FG2,x)';
     FG2pb =interp1(x+dx/2,FG2p,x)';
     FG2ppb=interp1(x+dx/2,FG2pp,x)';
     FG3b  =interp1(x+dx/2,FG3,x)';
 else % circ conv and the conv theorem need full shift
     FG2b  =FG2;
     FG2pb =FG2p;
     FG2ppb=interp1(x+dx,FG2pp,x,)';
     FG3b  =interp1(x+dx,FG3,x,)';
 end

%% Plotting

mov=0.0001;
figure(FIG), clf, hold on
plot(x,f,'b'), plot(x,g,'r'), plot(x,fg,'k')
%By hand:
plot(x,FG1+mov*5,'c')
%By conv commands
% plot(x,FG2 -mov,'g')
plot(x,FG2b+mov,'g-.')
% plot(x,FG2p -2*mov,'y')
plot(x,FG2pb+2*mov,'y-.')
% plot(x,FG2pp -3*mov,'y','Linewidth',1.1)
plot(x,FG2ppb+3*mov,'y-.','Linewidth',1.1)
% Convolution theorem
% plot(x,FG3-1.5*mov,'m')
plot(x,FG3b+1.5*mov,'m-.')
% Particular case at f*g[t=0]
plot(0,trapz(x,f.*g2),'r*')

%% Error

fprintf('Errors:\n');
%By hand
error_hand=max(abs(FG1-fg))
% conv commands
%   max(abs(FG2-fg))   
error_conv2=max(abs(FG2b-fg))    
%   max(abs(FG2p-fg))  
error_conv=max(abs(FG2pb-fg))   
%   max(abs(FG2pp-fg)) 
error_cconv=max(abs(FG2ppb-fg))  
%   max(abs(FG3-fg))  
error_theo=max(abs(FG3b-fg)) 
% In case we do not have exact answer
fprintf('Number of points N+1, with N=%d\n',N);
fprintf('Approx of all convolution calculations \n')
%   max(abs(FG2-FG1))   %a little off N odd
approx_conv2=max(abs(FG2b-FG1))    
%   max(abs(FG2p-FG1))  %a little off N odd
approx_conv=max(abs(FG2pb-FG1))   
%   max(abs(FG2pp-FG1)) %a little off N odd
approx_cconv=max(abs(FG2ppb-FG1))  
%   max(abs(FG3-fg))  %a little off N odd/even
approx_theo=max(abs(FG3b-FG1)) % large doamin N odd
% Between Each other Conv
compare1_2bvs2pb=max(abs(FG2pb-FG2b))
compare2_2bvs2ppb=max(abs(FG2ppb-FG2b))
compare3_2pbvs2ppb=max(abs(FG2pb-FG2ppb))
compare4_3bvs2b=max(abs(FG3b-FG2b))
compare5_3bvs2pb=max(abs(FG3b-FG2pb))
compare6_3bvs2ppb=max(abs(FG3b-FG2ppb))

Answer 1

[Matt J]

First: I notice that we need to multiply by the dx, the discretization of space. Can anyone explain this to me?

Continuous convolutions are integrals while discrete convolutions are sums. When approximating integrals by sums, you multiply the sum by a sampling interval, thus obtaining a Riemann sum.

Second: I find that convolution needs to be re-interpolated in some cases, and I wonder why this need to happen to be more accurate?

The interpolation you are doing is equivalent to shifting by dx/2. If you think accuracy improved, it probably means your axis origin is dx/2 from where you think it is. Note, for example, that when mod(N,2)=1 and mu=0, you are not sampling the peak of the Gaussians. If you are pretending one of your samples is the peak regardless, that would explain why you are dx/2 off.

Third: The following Code performs the convolution of two gaussians in different ways but for N=odd, that is (N+1 even) number of points, we have less accurate. Thanks in advance

Again, possibly because when N+1 is odd, there are cases when the peak is not sampled.

[Carlos Prieto]

Hi Matt,

thanks for your input. I really made me think a little about the fact not sampling the peak correctly.

If set mu_1 and mu_2 equal to zero to have the gaussian with peaks at zero.

For N odd (where my issue lies), the number of points are even, and it is clear that I am not sampling the peak. But with mu's equals to dx/2, both functions are sampled at the peak.

I still obtain a larger error 10^-5, compare to N even that has errors of 10^-16. Isn't that odd (as in weird)? Something else must be going on.

Carlos.

[Matt J]

Dunno. Some questions to consider:

1. Do all the convolution methods show the same large error?
2. Have you checked that all convolution methods put the peak in the same place?
3. Have you checked that all convolution results are symmetric about the peak?
4. What happens to the error as you increase N (not from N to N+1, but to something much larger)?