MATLAB Answers

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

조회 수: 1(최근 30일)
Rebecca Thomas
Rebecca Thomas 2020년 11월 20일
댓글: solofthedibs 2020년 11월 20일
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)));
  댓글 수: 1
Bruno Luong
Bruno Luong 2020년 11월 20일
"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

댓글을 달려면 로그인하십시오.

채택된 답변

solofthedibs
solofthedibs 2020년 11월 20일
편집: solofthedibs 2020년 11월 20일
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
  댓글 수: 2
solofthedibs
solofthedibs 2020년 11월 20일
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

댓글을 달려면 로그인하십시오.

추가 답변(1개)

Jon
Jon 2020년 11월 20일
편집: Jon 2020년 11월 20일
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)

태그

Community Treasure Hunt

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

Start Hunting!

Translated by