How to replace two for loops with matrix expression

1 Ansicht (letzte 30 Tage)
Fei Lai
Fei Lai am 18 Mai 2016
Beantwortet: Nobel Mondal am 18 Mai 2016
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 Kommentar
Todd Leonhardt
Todd Leonhardt am 18 Mai 2016
Edit your post to format the code by indenting all code by at least 2 spaces. That will make it more readable.

Melden Sie sich an, um zu kommentieren.

Antworten (1)

Nobel Mondal
Nobel Mondal am 18 Mai 2016
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

Kategorien

Mehr zu Fortran with MATLAB finden Sie in Help Center und File Exchange

Tags

Noch keine Tags eingegeben.

Produkte

Community Treasure Hunt

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

Start Hunting!

Translated by