Clear Filters
Clear Filters

How to replace two for loops with matrix expression

1 view (last 30 days)
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 Comment
Todd Leonhardt
Todd Leonhardt on 18 May 2016
Edit your post to format the code by indenting all code by at least 2 spaces. That will make it more readable.

Sign in to comment.

Answers (1)

Nobel Mondal
Nobel Mondal on 18 May 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

Categories

Find more on Elementary Math in Help Center and File Exchange

Tags

No tags entered yet.

Products

Community Treasure Hunt

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

Start Hunting!