How do I vectorize a nested for loop with different sized steps?

Down below is an iterative formula. I have previously defined dx and dy as the resolution.

for iy=2:ny-1

for ix=2:nx-1

T_new(iy,ix)=(((dy^2)*(T_old(iy,ix+1)+T_old(iy,ix-1)))+((dx^2)*(T_old(iy+1,ix)+T_old(iy-1,ix))))/(2*((dx^2)+(dy^2)));

end

end

I wrote this code to replace my loops but I believe it doesnt't work as dx is not equal to dy i.e. my iteration sizes are different. Any suggestions?

ix=2:nx-1;

iy=2:ny-1;

T_new(iy,ix)=(((dy^2)*(T_old(iy,ix+1)+T_old(iy,ix-1)))+((dx^2)*(T_old(iy+1,ix)+T_old(iy-1,ix))))/(2*((dx^2)+(dy^2)));

[Bruno Luong]

"it doesnt't work as dx is not equal to dy i.e. my iteration sizes are different."

Why do you think dx and dy must be equal for one code to work and not another?

>> ny=4;
>> nx=6;
>> T_old=rand(ny,nx)

T_old =

    0.7431    0.7060    0.0971    0.9502    0.7655    0.4456
    0.3922    0.0318    0.8235    0.0344    0.7952    0.6463
    0.6555    0.2769    0.6948    0.4387    0.1869    0.7094
    0.1712    0.0462    0.3171    0.3816    0.4898    0.7547

>> T_new=zeros(ny,nx);
>> dx=rand; dy=rand;
>> for iy=2:ny-1
    for ix=2:nx-1
        T_new(iy,ix)=(((dy^2)*(T_old(iy,ix+1)+T_old(iy,ix-1)))+((dx^2)*(T_old(iy+1,ix)+T_old(iy-1,ix))))/(2*((dx^2)+(dy^2)));
    end
end
>> disp(T_new)
         0         0         0         0         0         0
         0    0.5914    0.0845    0.7931    0.3596         0
         0    0.5851    0.3879    0.4079    0.5837         0
         0         0         0         0         0         0

>> T_new=zeros(ny,nx);
>> ix=2:nx-1;
>> iy=2:ny-1;
>> T_new(iy,ix)=(((dy^2)*(T_old(iy,ix+1)+T_old(iy,ix-1)))+((dx^2)*(T_old(iy+1,ix)+T_old(iy-1,ix))))/(2*((dx^2)+(dy^2)));
>> disp(T_new)
         0         0         0         0         0         0
         0    0.5914    0.0845    0.7931    0.3596         0
         0    0.5851    0.3879    0.4079    0.5837         0
         0         0         0         0         0         0

Answer 1

[solofthedibs]

What makes you think that it doesn't work?

Edit for clarity: I just copy-pasted your two code snippets, and generated random dimensions and step sizes to try them out.

>> dx = rand() % random step size
dy = rand() % random step size
nx = randi([100 1000], [1 1]) % random 1st dimension size
ny = randi([100 1000], [1 1]) % random 2nd dimension size
T_old = rand([ny nx]); % random starting matrix
% compute T_new the loopy way as T_new1
for iy=2:ny-1
for ix=2:nx-1
T_new1(iy,ix)=(((dy^2)*(T_old(iy,ix+1)+T_old(iy,ix-1)))+((dx^2)*(T_old(iy+1,ix)+T_old(iy-1,ix))))/(2*((dx^2)+(dy^2)));
end
end
T_new2(iy,ix)=(((dy^2)*(T_old(iy,ix+1)+T_old(iy,ix-1)))+((dx^2)*(T_old(iy+1,ix)+T_old(iy-1,ix))))/(2*((dx^2)+(dy^2)));
% compute T_new the MATLAB way
ix=2:nx-1;
iy=2:ny-1;
T_new2(iy,ix)=(((dy^2)*(T_old(iy,ix+1)+T_old(iy,ix-1)))+((dx^2)*(T_old(iy+1,ix)+T_old(iy-1,ix))))/(2*((dx^2)+(dy^2)));
% check agreement - a value of zero indicates they are identical
max(reshape(abs(T_new2 - T_new1), 1, []))

dx =

    0.9722


dy =

    0.4396


nx =

   938


ny =

   748


ans =

     0

>> % ans = 0 -> identical

[Jon]

I guess this was also Bruno Luong's point in his comments

[solofthedibs]

Oh, yes. Exactly.

When I posted my answer there was only the first two lines of his comment, asking why dx~=dy should matter.

At least, that I noticed.

/shrug

Answer 2

[Jon]

I think you may be able to utilize the filter2 function for this purpose. For example

A = [1 2 3 4;5 6 7 8;9 10 11 12;13 14 15 16]
A =

     1     2     3     4
     5     6     7     8
     9    10    11    12
    13    14    15    16
    
H = [0 1 0;1 0 1;0 1 0]
Af = filter2(H,A,'valid')

Af =

    24    28
    40    44

You will still need to apply your scaling with dx and dy but I think this should get you close.

Note the filter matrix

H = [0 1 0;1 0 1;0 1 0]

Tells it to compute the filtered element as the sum of the neighbors to the left and right and above and below, which I think is what you want. The 'valid' computes only interior points, as for example A(m-1,2) is not defined for m = 1. Actually they zero pad the edges and then compute the whole thing if you want that.

I realize this isn't exactly what you asked for, you wanted to vectorize it yourself. I'm not sure how you would do this without the double loop, maybe someone else has an idea, but in anycase I would assume that this built in MATLAB function is efficient, and keeps your code clean of complicated double loops (maybe they are in the built in function but you don't have to see them anyhow)