How to replace two for loops with matrix expression

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?

[Todd Leonhardt]

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]

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