using trapz in a loop or not?

Hi,

I've been using trapz to calculate numerical integral of SL as below.

flip_MT  = [212 434 639 843];
tmp_MT   = repelem(flip_MT',50);
df = repmat(logspace(0,5)',length(flip_MT),1);
theta  = 0:0.01:pi/2;
T2B    = 18.3E-6;
thetam = meshgrid (theta , df);  
dfm      = meshgrid (df , theta)';
f           = 3*cos(thetam).^2-1;     
SL1  = sin(thetam)*sqrt(2/pi)*T2B./abs(f).*exp(-2.*(2*pi*dfm*T2B./f).^2); 
ret1 = trapz(theta,SL1,2)'; 

for i = 1:length(df)
    SL2 (i,:)   = sin(theta).*sqrt(2/pi).*T2B./abs(f(i,:)).*exp(-2*(2*pi*df(i).*T2B./f(i,:)).^2);
    ret2(i,:)   = trapz(theta,SL2(i,:));
end
SLd = SL1-SL2; %SL1 and SL2 are equal
retd = ret1-ret2';%different

I have noticed that ret1 and ret are different, while basically they have to be the same, since SL1 and SL2 are the same.

I figured trapz behaves differently inside and outside a loop.

My question is that which answer is correct now? ret1 or ret2?

[Roya Afshari]

I noticed that the order of magnitude of ret1 and ret2 differences(retd) is neglectable compared to the ret values. Keep in mind that there is still a difference about 1E-10 when your numbers are in the order of 1E-5. Otherwise both methods work fine.

[Looky]

This is most likely a common numerical problem. In both cases what matlab calculates is:

z = diff(x,1,1).' * (y(1:m-1,:) + y(2:m,:))/2;

In your loop(ret2), this becomes a vector * vector operation, most likely using the underlying Blas/Lapack libs. Such a vector vector operation is Blas Level 1 and to say it simple, there isn't much magic happening. Most likely straight forward computation, maybe such sharing between your cpu cores, you get what you see in most cases.

However, if you compute it as a vector matrix computation it gets a bit different. It's now Blas Level 2 and thats where you can do some computational/numerical magic. Basically the compiler can optimize your problem much better and you have a better balance between memory access and computational effort. However, this leads to a somewhat different way of dealing with the problem in terms of the math involved. Basically you sum/multiply stuff in a different order. This wouldn't make any difference in mathematical terms, but since your computer has to round stuff after most operations (to fit the result in a double type), the round off errors sum in a different manner.

What is wrong or right? No idea, after all, it is just a different way of summing the pieces. What is the faster? The matrix method, by a long shot most likely! Eitherway, the differences should not be very significant.

[Roya Afshari]

Thank you Looky. I agree. In the end, the difference doesn't influence my results that much. I was just wandering what happened there.