Caught in loop
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Jan
2012-5-7
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
更多回答(1 个)
Richard Brown
2012-5-7
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.
2 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!