using trapz in a loop or not?

Question

Roya Afshari 2020-9-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/597253-using-trapz-in-a-loop-or-not

评论： Roya Afshari 2020-9-22

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?

3 个评论
显示 1更早的评论隐藏 1更早的评论

Looky 2020-9-22

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 2020-9-22

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

请先登录，再进行评论。

请先登录，再回答此问题。

using trapz in a loop or not?

3 个评论
显示 1更早的评论隐藏 1更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

using trapz in a loop or not?

3 个评论 显示 1更早的评论隐藏 1更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

3 个评论
显示 1更早的评论隐藏 1更早的评论