Iteration coded weird behaviour

Question

Dejan Cvijanovic 2012-5-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/37802-iteration-coded-weird-behaviour

Sorry to bring up my old code again. The problem I'm having is the radius [a(i)] decreases to 0 (as it should), but the concentration [c(i,j)] decreases, in general, but unfortunately increase on some individual time steps (e.g; 0.5 0.4 0.45 0.3 0.32 0.12 0. 18 ...), when it should be strictly decreasing (e.g; 0.5 0.4 0.3 ....)

...I'm thinking the problem is either in c_ah or the definition for c(i,j+1). Any thoughts?

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

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

r = (0:N_i) * dr; c = zeros(N_i, N_t); c(:, 1) = c_inf; c(r < a_0, 1) = c_s;

a = zeros(1,N_t); a(1) = a_0;

dr3 = dr ^ 3;

 for j = 1:N_t -1
     if a(j) <= 0
         a(j) = 0;
         ajp1 = 0;
         a(j+1) = 0;
     else
  b_1 = floor (a(j)/dr) + 3;
  b_2 = b_1 + 1;
  b_1s  = (b_1 - 1)*dr;
  b_2s = (b_2 - 1)*dr;
%   c_ah = (a(j)+h-b_1s)*(a(j)+h-b_2s)*(c_a)/((a(j)-b_1s)*(a(j)-b_2s)) + ...
%       h*(a(j)+h-b_2s)*c(b_1,j)/((b_1s-a(j))*(-1)) + ...
%       h*(a(j)+h-b_1s)*c(b_2,j)/((b_2s-a(j))*(1));
  c_ah = (a(j)+h-b_1s)*(a(j)+h-b_2s)*(c_a)/((a(j)-b_1s)*(a(j)-b_2s)) + ...
      h*(a(j)+h-b_2s)*c(b_1,j)/((b_1s-a(j))*(b_1s-b_2s)) + ...
      h*(a(j)+h-b_1s)*c(b_2,j)/((b_2s-a(j))*(b_2s-b_1s));
  dc_dr = ((c_ah)-c_a)/h;
  ajp1   = a(j) +dt*beta*dc_dr;
  a(j+1) = ajp1;
     end
  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(' num2str(i) ',' num2str(j) ') = ' num2str(c(i,j))]);
    disp(['a(j) = ' num2str(a(j))]);
  end
  c(1, j+1) = c(2, j+1);
  c(N_i, j+1) = c_inf;
 end