Caught in loop

My code is as follows:

alpha = 0.5;
beta = 0.1;
a_0 = 1;
dr = 0.1;
dt = input('Value of dt:');
h = 0.001;
c_inf = 0;
c_a = 1;
c_s = 2;
mesh_points = 3;
rho = dt/(dr)^2;
R= 30;
N_i = R/dr;
t_diss = 5.8;
N_t = t_diss/dt;

if rho >= 0.5;
    disp('Invalid dt value');
    return;
end

for i=1:N_i;
    r(i)=(i-1)*dr;
    if r(i)< a_0
        c(i,1)= c_s;
    else 
        c(i,1)= c_inf;
    end
end
a(1) = a_0;

for j = 1:N_t -1;
    b_1 = floor (a(j)/dr) + 2;
    b_2 = b_1 + 1;
    c_ah = (a(j)+h-b_1)*(a(j)+h-b_2)*(c_a)/((a(j)-b_1)*(a(j)-b_2))+h*(a(j)+h-b_2)*c(b_1,j)/((b_1-a(j))*(b_1-b_2)) + h*(a(j)+h-b_1)*c(b_2,j)/((b_2-a(j))*(b_2-b_1));
    dc_dr = ((c_ah)-c_a)/h;
    a(j+1) = a(j) +dt*beta*dc_dr;

     for i=2:N_i - 1;
         if r(i) > a(j+1);
             d = rho - (rho/(i-1)) - (a(j)^2)*(alpha-1)*(a(j+1)-a(j))/(2*(i-1)^2*(dr)^3);
             e = 1-2*rho;
             f = rho + (rho/(i-1)) + (a(j)^2)*(alpha-1)*(a(j+1)-a(j))/(2*(i-1)^2*(dr)^3);
             c(i, j+1) = d*c(i-1,j) + e*c(i,j) + f*c(i+1, j);

         elseif r(i) == a(j+1);
                 c(i,j+1) = c_a;
             else c(i,j+1) = c_s;
         end
             disp(c(i,j));
         end
         c(1, j+1) = c(2, j+1);
  end

It basically calculates concentration near a dissolving sphere but when it gets to disp(c(i,j)) I get caught in loop of displaying numbers that change periodically. Any ideas?

Answer 1

[Jan]

And some cleanup:

...
% Vectorization of the "for i=1:N_i;" loop:
r = (0:N_i) * dr;
c = zeros(N_i, N_t - 1);  % Pre-allocate!
c(:, 1)       = c_inf;
c(r < a_0, 1) = c_s;
...
dr3 = dr ^ 3;
...
  % Reduce repeated calculations inside the "for j = 1:N_t -1" loop:
  ajp1   = a(j) +dt*beta*dc_dr;
  a(j+1) = ajp1;
  tmp1   = (a(j)^2) * (alpha-1) * (ajp1-a(j)) / (2 * dr3);
  for i = 2:N_i - 1
    if r(i) > ajp1
      tmp2 = rho / (i-1) + tmp1 / ((i-1)^2);
      d = rho - tmp2;
      e = 1 - 2 * rho;
      f = rho + tmp2;
      c(i, j+1) = d*c(i-1,j) + e*c(i,j) + f*c(i+1, j);
    elseif r(i) == ajp1
      c(i, j+1) = c_a;
    else
      c(i, j+1) = c_s;
    end
    % disp(c(i,j));
  end

[Richard Brown]

What service!

[Jan]

It is relaxing to improve some Matlab code instead of rererepeating "please format your code" and "please post any details".

[Dejan Cvijanovic]

Thanks. Cleaned up nicely indeed.

[Dejan Cvijanovic]

Got an error: Attempted to access a(1); index out of bounds because numel(a)=0.

Error in CleanDiss (line 36)

ajp1 = a(j) +dt*beta*dc_dr;

Ideas?

[Dejan Cvijanovic]

I fixed the error and get these numbers (which make sense):

0.4011

0.3037

0.3882

0.2886

0.3634

0.2649

0.3293

0.2350

0.2886

0.2014

0.2447

0.1667

0.2006

0.1333

0.1589

0.1028

0.1217

0.0766

0.0900

0.0550

0.0643

0.0381

0.0443

0.0254

0.0294

0.0164

0.0189

0.0101

They generally decrease, which is correct but some of them increase, which shouldn't be true?

Answer 2

[Richard Brown]

Your code works fine. If you remove the disp statement (which is causing it to slow down heaps), you'll find that it runs completely in around 10 seconds or so.

[Richard Brown]

Oh, and a hint: Select all the code, then Text->Smart Indent

[Image Analyst]

Cool. Nice trick that I didn't know about.