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

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

Answer 1

[Stephen23]

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:

https://www.mathworks.com/help/matlab/matlab_prog/vectorization.html

[dpb]

I edit'ed to clean up formatting/indention a little to more easily read while looking and inserted it, Stephen...

[L]

Thank you, I see now I forgot the semicolon -- silly mistake!

Answer 2

[Walter Roberson]

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.

Answer 3

[dpb]

...
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...