MATLAB Answers

Translated by

このページのコンテンツは英語から自動翻訳されています。自動翻訳をオフにする場合は「<a class="turn_off_mt" href="#">ここ</a>」をクリックしてください。

L
0

Where am I going wrong in my for loop? I am getting NaNs

L
さんによって質問されました 2018 年 10 月 10 日
最新アクティビティ dpb
さんによって 回答されました 2018 年 10 月 10 日
I just added the for loop labeled "%surface boundary flux" (the first for loop within the main for loop) and am now getting NaN. The code worked with the previous boundary condition. New to MATLAB, am I doing something clearly wrong? Thanks!
% Parameters
T_av = 250; %[K] initial temperature
Tamp = 150; %[K] impulse temperature at surface layer
conductivity=.0033; %W m-1 K-1
P = 2.55024e6; %[s] period of sinusoidal temperature change
synodic_frequency=(2*pi)/P;
K = 3.7786e-09; %[m2.s-2] thermal diffusivity
dz = 0.001; %[m] depth step
dt = (0.50)*0.5*dz^2/K; %[s] time step where first fraction must be less than 1, comp phy book
J = 700; %number of depth steps
N = 1000000; %number of time steps
flux=0;
% z,t,T
z = 0:dz:dz*J; %depth profile
t = 0:dt:dt*N; %time profile
T_num = nan(length(z),length(t)); %temperature mesh
T_num(:,1) = T_av; %initial temperature is Tavg
for n = 1:N-1 %iterating in time
for sur_flux = 0:1350
T_num(1,n+1) = -sur_flux*dz/conductivity + 250; %surface boundary flux
end
for j = 2:J-1
T_num(j,n+1)=T_num(j,n)+(K*dt/dz^2)*(T_num(j+1,n)-2*T_num(j,n)+T_num(j-1,n)); %Iterature temperature over j and n
end
T_num(J,n+1) = -flux*dz/conductivity+T_num(J-1,n+1); %bottom boundary flux
end

  0 件のコメント

サインイン to comment.

タグ

3 件の回答

回答者: Stephen Cobeldick
2018 年 10 月 10 日
編集済み: Stephen Cobeldick
2018 年 10 月 10 日

The only problem I can see is that you forgot to put the semi-colon on the end of the line.
That will display T_num on every loop iteration, which will look as if T_num contains mostly NaN's... because at that point in the calculation it does mostly contain NaN's! You should add the semi-colon and check the matrix after the code has finished running.
When I try your code T_num's last row consists of NaN's, but this is due to a bug(?) in your code, unrelated to the sur_flux loop.
Note that the loops are probably not required anyway:

  2 件のコメント

dpb
2018 年 10 月 10 日
I edit'ed to clean up formatting/indention a little to more easily read while looking and inserted it, Stephen...
L
2018 年 10 月 10 日
Thank you, I see now I forgot the semicolon -- silly mistake!

サインイン to comment.


回答者: Walter Roberson
2018 年 10 月 10 日

You have
for sur_flux = 0:1350
T_num(1,n+1) = -sur_flux*dz/conductivity + 250; %surface boundary flux
end
which writes to the same location every time, since the index at the left does not change. Perhaps you mean
for sur_flux = 0:1350
T_num(1,sur_flux+1) = -sur_flux*dz/conductivity + 250; %surface boundary flux
end
If you did, then
T_num(1, 1:1351) = -(0:1350)*dz/conductivity + 250;
with no loop.
However, the relationship between 1350 and the other sizes is not at all clear, so it is not obvious what you are trying to do in the statement.

  0 件のコメント

サインイン to comment.


回答者: dpb
2018 年 10 月 10 日

...
T_num = nan(length(z),length(t)); % initialized to NaN
T_num(:,1) = T_av; %initial temperature is Tavg
for n=1:N-1
for sur_flux=0:1350
T_num(1,n+1)=-sur_flux*dz/conductivity + 250; %surface boundary flux
end
for j = 2:J-1
T_num(j,n+1)=T_num(j,n)+(K*dt/dz^2)*(T_num(j+1,n)-2*T_num(j,n)+T_num(j-1,n));
end
T_num(J,n+1) = -flux*dz/conductivity+T_num(J-1,n+1); %bottom boundary flux
end
The above first sets first column to T_av
The loop on sur_flux isn't what you intend I think; it sets the same location of (1,n+1) sequentially to the computed result but does that without changing n so all that happens before the j loop is to set
T_num(1,n+1)=-1350*dz/conductivity + 250;
all the other values are immaterial because the same node is overwritten every pass. Not exactly sure what your intent is, but that undoubtedly isn't it... :)
Then with n=1, for j=2,
T_num(j,n+1)=T_num(j,n)+(K*dt/dz^2)*(T_num(j+1,n)-2*T_num(j,n)+T_num(j-1,n));
is
T_num(2,2)=T_num(2,1)+(K*dt/dz^2)*(T_num(3,1)-2*T_num(2,1)+T_num(1,1));
Hmmmm...that looks ok; I thought there was going to be a reference to an uninitialized location that was going to propagate through.
I still think that must be so and I missed it...set a breakpoint and step through with the debugger and see where the NaN comes from...

  0 件のコメント

サインイン to comment.



Translated by