Iteration coded weird behaviour

1 次查看(过去 30 天)
Dejan Cvijanovic
Dejan Cvijanovic 2012-5-8
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

回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Surface and Mesh Plots 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by