Trapezoidal Rule of Convolution with Non-Uniform Intervals

Hi all,

I'm supposed to do the numerical integration of a convolution for t, which is given by specific non-uniform timepoints. For t, there is also a function Cp which is a value associated with the t value.

Here's what I have so far. It doesn't work but I hope I'm on the right track?

Thanks in advance!

t=4.80;
k1=0.102;
k2=0.130;
k3=0.062;
k4=0.0068;
alpha1=(k2+k3+k4+(sqrt((k2+k3+k4)^2-(4*k2*k4))))/2;
alpha2=(k2+k3+k4-(sqrt((k2+k3+k4)^2-(4*k2*k4))))/2;
A=(k1*(k4+k3-alpha1))/(alpha2-alpha1);
B=(k1*(alpha2-k4-k3))/(alpha2-alpha1);

function trape=traperule(t)
        time=[0,1.08,1.78,2.30,2.75,3.30,3.82,4.32,4.80,5.28,5.95,6.32,6.98,9.83,16.30,20.25,29.67,39.93,58,74,94,100,200,300,400,500,591];
        [~,pos] = ismembertol(t, time, 1E-8)
        tp=time(1:pos)
        conc=[0,84.9,230,233,220,236.4,245.1,230.0,227.8,261.9,311.7,321,316.6,220.7,231.7,199.4,211.1,190.8,155.2,140.1,144.2,139.737,111.3006,82.8375,54.3744,25.9113,0.0098];
        C=conc(1:pos);
        k1=0.102;
k2=0.130;
k3=0.062;
k4=0.0068;
alpha1=(k2+k3+k4+(sqrt((k2+k3+k4)^2-(4*k2*k4))))/2;
alpha2=(k2+k3+k4-(sqrt((k2+k3+k4)^2-(4*k2*k4))))/2;
A=(k1*(k4+k3-alpha1))/(alpha2-alpha1);
B=(k1*(alpha2-k4-k3))/(alpha2-alpha1);
trape=0;
        for n=1:pos
            tprime=tp(n)
            if n>1
            told=tp(n-1)
            else
                told=0;
            end
            Cp=C(n);
            fun = @(tprime) 0.0675*Cp.*exp(-alpha1*(t-tprime)+Cp.*0.0345*exp(-alpha1*(t-tprime)));
        Ci=integral(fun,told,tprime);
        h=tprime-told;
        area = Ci * h / 2 % area of the trapezoid
        trape = trape + area;
        end
    end

Answer 1

[Catalytic]

I don't believe there is such a thing as the "Trapezoidal Rule of Convolution". Undoubtedly, you are thinking of the trapezoidal rule of integration, however, that will not help you deal with the non-uniform time sampling. Trapezoidal integration assumes you can approximate the integrand with linear interpolation. However, you do not have samples of the convolution integrand with which to interpolate. You could linearly interpolate the two signals you are trying to convolve and then multiply them together, but that would be a different kind of integral apprxoimation.

Answer 2

[Matt J]

It is not clear from your post which two functions you are convolving, so I will assume here that it is the same as your other recent post here. It does not give significantly different results as the other post, however, because in both posts the signals being convolved are approximated at some stage by trapezoids.

t=[0,1.08,1.78,2.30,2.75,3.30,3.82,4.32,4.80,5.28,5.95,6.32,6.98,9.83,16.30,20.25,29.67,39.93,58,74,94,100,200,300,400,500,591];

%v denotes the exponential function
v=myfunction(t);
%conce is the function for Cp(t)
cp=Cp(t);

vlin=@(tq) interp1(t,v,tq,'linear',0);
cplin=@(tq) interp1(t,cp,tq,'linear',0);

dt=min(diff(t));
tnew=min(t):dt:2*max(t);

fun=@(tau) integral( @(t) vlin(t).*cplin(tau-t), min(t),max(t) );

Ci=arrayfun(fun,tnew);

%CC is convolution of Cp(t)*the exponential function

plot(t,cp,tnew,Ci);
title("Concentration over Time");
xlabel('Time (mins)');
ylabel('Concentration')
legend("Cp(t)","Ci(t)");

%Based on plot, it appears that the highest signal strength in the brain is
%around 39.93 minutes, which should be the best time to run the PET scan.
function conce=Cp(a)
t=[0,1.08,1.78,2.30,2.75,3.30,3.82,4.32,4.80,5.28,5.95,6.32,6.98,9.83,16.30,20.25,29.67,39.93,58,74,94,100,200,300,400,500,591];
conc=[0,84.9,230,233,220,236.4,245.1,230.0,227.8,261.9,311.7,321,316.6,220.7,231.7,199.4,211.1,190.8,155.2,140.1,144.2,139.737,111.3006,82.8375,54.3744,25.9113,0.0098];
conce=conc(t==a);
end
function v=myfunction(t);
k1=0.102;
k2=0.130;
k3=0.062;
k4=0.0068;
alpha1=(k2+k3+k4+(sqrt((k2+k3+k4)^2-(4*k2*k4))))/2;
alpha2=(k2+k3+k4-(sqrt((k2+k3+k4)^2-(4*k2*k4))))/2;
A=(k1*(k4+k3-alpha1))/(alpha2-alpha1);
B=(k1*(alpha2-k4-k3))/(alpha2-alpha1);
v=A*exp(-alpha1*t)+B*exp(-alpha2*t);
%Ci=Cp(t).*(A.*exp(-1*alpha1*t)+B.*exp(-1*alpha2*t));
end

[Lilly]

Hi, yes it is the same two functions mentioned in that post.

I was trying to do the numerical integration using the trapezoidal rule to determine the overall concentration (Ci) as a value.

What should I change in my code to do so?

For this specific instance, the convolution in question being the integral of fun = @(tprime) A*Cp.*exp(-alpha1*(t-tprime)+Cp.*B*exp(-alpha2*(t-tprime)));

Let me know if this doesnt make sense

[Matt J]

I don't understand what your code is doing. Why isn't my proposed code a suitable replacement?

[Lilly]

My code is trying to find the total concentration via the trapezoidal rule when t=4.80.

I was trying to accomplish that but I'm sure there are some things that I need to change in my code.

[Matt J]

If you are only interested in t=4.80 specifically, my code can be modified to,

Ci=arrayfun(fun,4.80);

Similar to what @Catalytic said, though, the trapezoidal rule doesn't apply to your problem in the usual sense, so whatever your own code is trying to do, it seems futile. To apply the trapezoidal rule, you would need to break the convolution integrand into trapezoids. That is not possible for you because, although you have and t a common set of discrete time points, you do not have and at a common set of time points.

[Lilly]

For this problem, the specific ask is that a method of numerical integration be used to solve for Ci(t). That's the only reason I am hesitant to implement your code.

[Matt J]

Any integration method used here is going to be a numerical method. There is no closed analytical form for the convolution of these signals.

[Lilly]

In that case, thank you so much for your help! I appreciate it

[Paul]

Assuming that Cp(t) is defined by the linear interpolation between the points, and that Cp(t) = 0 for t < 0 and for t > 591, and that v(t) = 0 for t < 0 (which are the same assumptions used above), then a closed form expression for the convolution Cp(t)*v(t) can be obtained.

tconc = [0,1.08,1.78,2.30,2.75,3.30,3.82,4.32,4.80,5.28,5.95,6.32,6.98,9.83,16.30,20.25,29.67,39.93,58,74,94,100,200,300,400,500,591];
conc  = [0,84.9,230,233,220,236.4,245.1,230.0,227.8,261.9,311.7,321,316.6,220.7,231.7,199.4,211.1,190.8,155.2,140.1,144.2,139.737,111.3006,82.8375,54.3744,25.9113,0.0098];

syms t real
v(t) = myfunction(t);

syms Cp(t)
Cp(t) = 0;
for ii = 1:numel(tconc)-1
    Cp(t) = Cp(t) + trapezoidalPulse(tconc(ii),conc(ii),tconc(ii+1),conc(ii+1),t);
end

figure
fplot(Cp(t),[0 600])
hold on

syms tau real
Ci(t) = int(Cp(tau)*v(t-tau),tau,0,t);

% use vpa to make the display marginally readable
vpa(Ci(t),5)

ans = 

fplot(Ci(t),[0 1200])

In the interval 0 < t < 591, this solution for Ci(t) is very close to that developed by @Matt J with a small difference presumably due to the difference between the exact solution returned by int and the approximate solution returned by integral and that Matt J's solution used linear interpolation for v(t). However, for t > 591 the solutions diverge. I believe that's because Matt J's solution assumes that v(t) = 0 for t > 591 whereas this soution uses the definition of v(t) as given, which has infinite duration.

Matt J's solution could be modified so that the integrand would be defined as

@(t) myfunction(t).*cplin(tau-t)

and then the solutions would (presumably, I didn't try it) match up quite well. Or the solution here could (presumably, I didn't try it) be modified to window v(t) over the inteval 0 < t < 591 if that's the appropriate assumption.

function y = trapezoidalPulse(x1,y1,x2,y2,x)

y = (y2 - y1)/(x2 - x1)*(x - x1) + y1;
y = y*rectangularPulse(x1,x2,x);

end

function v=myfunction(t);
k1=0.102;
k2=0.130;
k3=0.062;
k4=0.0068;
alpha1=(k2+k3+k4+(sqrt((k2+k3+k4)^2-(4*k2*k4))))/2;
alpha2=(k2+k3+k4-(sqrt((k2+k3+k4)^2-(4*k2*k4))))/2;
A=(k1*(k4+k3-alpha1))/(alpha2-alpha1);
B=(k1*(alpha2-k4-k3))/(alpha2-alpha1);
v=A*exp(-alpha1*t)+B*exp(-alpha2*t);
end