MATLAB: Conv not giving the expected results. Why

convconvolutionintegral

My mathematical model consists of a convolution. I tried the the MATLAB "conv" function and it does not seem to return the expected results. I also tried to implement the convolution using "integral" (using formula below) and it works.
Please advise on how to make "conv" work because it is much faster and I need to calculate the convolution quite a lot.
t = 0:5/60:100/6;
% Method 1: conv
y1 = 0.3*conv(helper(t)/0.6,exp(-0.3*t));
%% Method 2: integral
fun = @(tt,Tau) helper(tt)/0.6.*exp(-0.3*(Tau-tt));
for i=1:length(t)
    y2(i) = 0.3*integral(@(ttt)fun(ttt,t(i)),0,t(i),'RelTol', 1e-4, 'AbsTol', 1e-6);
end
%% Plotting results
figure; plot(t,y1(1:length(t)),t,y2);
%% Helper function
function y = helper(t) 
    syms x(t_fun)
    x(t_fun) =    5.7326*exp(-(t_fun-0.17046)^2/0.0063) ...
                         + 0.9974*exp(-(t_fun-0.365)^2/0.0348)  ...
                         + 1.050 *exp(-0.1685*t_fun)/(1+exp(-38.078*(t_fun-0.483)));
    y = double(x(t));
end

Best Answer

Hello Tara,
You want to approximate an integral by using a sum. In general, for equally spaced t values with spacing deltat, then
Int f(t) dt --->  Sum{k} f(t_k)*deltat
i.e. deltat is the width of the rectangles that are summed up to approximate the Riemann integral. Conv is doing the sum, so If you multiply that by (in this case) 5/60 the two methods agree fairly well. The integral does a bit better because it can dynamically reduce the spacing. In the example below I reduced deltat to 1/60 and the two methods agree better.
In your helper function, there is absolutely no need to go to symbolic variables. All you do is turn the result into a double anyway. It's better to just calculate the function directly, in which case all you need to do is replace
^ by .^,  * by .*,  /by ./
in the appropriate places so operations on the vector t are done term-by-term as below.
deltat = 1/60;
t = 0:deltat:100/6;
% Method 1: conv
y1 = 0.3*conv(helper(t)/0.6,exp(-0.3*t))*deltat;
% Method 2: integral
fun = @(tt,Tau) helper(tt)/0.6.*exp(-0.3*(Tau-tt));
y2 = zeros(size(t));
for i=1:length(t)
    y2(i) = 0.3*integral(@(ttt)fun(ttt,t(i)),0,t(i),'RelTol', 1e-4, 'AbsTol', 1e-6);
end
% Plotting results
figure(1)
plot(t,y1(1:length(t)),t,y2);
%Helper function
function y = helper(t)
    y =    5.7326*exp(-(t-0.17046).^2/0.0063) ...
                         + 0.9974*exp(-(t-0.365).^2/0.0348)  ...
                         + 1.050 *exp(-0.1685*t)./(1+exp(-38.078*(t-0.483)));
end
Related Question