Caught in loop

Question

Dejan Cvijanovic 2012-5-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/37620-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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Jan 2012-5-7

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/37620-caught-in-loop#answer_46935

在 MATLAB Online 中打开

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

Dejan Cvijanovic 2012-5-7

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 2012-5-7

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 2012-5-7

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/37620-caught-in-loop#answer_46924

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.