How to replace two for loops with matrix expression
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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
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.
Antworten (1)
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
0 Kommentare
Siehe auch
Kategorien
Mehr zu Fortran with MATLAB finden Sie in Help Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!