# Trapezoidal Rule of Convolution with Non-Uniform Intervals

23 views (last 30 days)

Show older comments

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

##### 0 Comments

### Accepted Answer

Catalytic
on 12 Oct 2023

##### 0 Comments

### More Answers (1)

Matt J
on 12 Oct 2023

Edited: Matt J
on 12 Oct 2023

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

##### 8 Comments

Paul
on 13 Oct 2023

Edited: Paul
on 13 Oct 2023

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)

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

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!