How to replace two for loops with matrix expression

Question

Fei Lai 2016-5-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/284700-how-to-replace-two-for-loops-with-matrix-expression

回答： Nobel Mondal 2016-5-18

the original code is like this:

 %
for i=2:nhx-1
    for j=2:nhy-1
      Unew(i,j)= U(i,j)-dt*(P(i+1,j)-P(i-1,j))/(2*hx)...
                 +nu*dt*(1/(hx*hx)*(U(i+1,j)-2.*U(i,j)+U(i-1,j))...
                        +1/(hy*hy)*(U(i,j+1)-2.*U(i,j)+U(i,j-1)))...
           -dt*U(i,j)/(hx)*(U(i,j)-U(i-1,j))...
                 -dt*V(i,j)/(2*hy)*(U(i,j+1)-U(i,j-1));
    end
  end

how can I replace two for loops to increase the calculation speed?

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

Todd Leonhardt 2016-5-18

Edit your post to format the code by indenting all code by at least 2 spaces. That will make it more readable.

请先登录，再进行评论。

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

Answer 1

Nobel Mondal 2016-5-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/284700-how-to-replace-two-for-loops-with-matrix-expression#answer_222557

在 MATLAB Online 中打开

You may directly use "vectorization" and "element-wise operation to solve this :

Here is an example code -

%%Assumptions - as I don't have the actual data
U = magic(100);
P = magic(length(U));
V = magic(length(U));
dt = 0.1;
hx = rand;
hy = rand;
nu = rand;
nhx = size(U,1);
nhy = size(U,2);
Unew = zeros(nhx, nhy);
%%Old algo
for i=2:nhx-1
    for j=2:nhy-1
      Unew(i,j)= U(i,j)-dt*(P(i+1,j)-P(i-1,j))/(2*hx)...
                 +nu*dt*(1/(hx*hx)*(U(i+1,j)-2.*U(i,j)+U(i-1,j))...
                        +1/(hy*hy)*(U(i,j+1)-2.*U(i,j)+U(i,j-1)))...
           -dt*U(i,j)/(hx)*(U(i,j)-U(i-1,j))...
                 -dt*V(i,j)/(2*hy)*(U(i,j+1)-U(i,j-1));
    end
end
%%Vectorized algo
myUnew = zeros(nhx, nhy);
myUnew(2:nhx-1, 2:nhy-1)= U(2:nhx-1, 2:nhy-1)-dt.*(P(3:nhx,2:nhy-1)-P(1:nhx-2,2:nhy-1))./(2*hx)...
                 +nu.*dt.*(1/(hx*hx).*(U(3:nhx,2:nhy-1)-2.*U(2:nhx-1, 2:nhy-1)+U(1:nhx-2,2:nhy-1))...
                        +1/(hy*hy).*(U(2:nhx-1,3:nhy)-2.*U(2:nhx-1, 2:nhy-1)+U(2:nhx-1,1:nhy-2)))...
           -dt.*U(2:nhx-1, 2:nhy-1)/(hx).*(U(2:nhx-1, 2:nhy-1)-U(1:nhx-2,2:nhy-1))...
                 -dt.*V(2:nhx-1, 2:nhy-1)/(2*hy).*(U(2:nhx-1,3:nhy)-U(2:nhx-1,1:nhy-2));
%%Now verify functionality
success = isequal(Unew, myUnew);
if success
    disp('Test passed')
else
    disp('Test failed. Please check the algorithm.')
end