How to replace two for loops with matrix expression

조회 수: 1 (최근 30일)
Fei Lai
Fei Lai 2016년 5월 18일
답변: 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
Todd Leonhardt
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.

댓글을 달려면 로그인하십시오.

답변 (1개)

Nobel Mondal
Nobel Mondal 2016년 5월 18일
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

카테고리

Help CenterFile Exchange에서 Fortran with MATLAB에 대해 자세히 알아보기

태그

아직 태그를 입력하지 않았습니다.

제품

Community Treasure Hunt

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

Start Hunting!

Translated by