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?
采纳的回答
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
5 个评论
Richard Brown
2012-5-7
What service!
Jan
2012-5-7
It is relaxing to improve some Matlab code instead of rererepeating "please format your code" and "please post any details".
Dejan Cvijanovic
2012-5-7
Dejan Cvijanovic
2012-5-7
Dejan Cvijanovic
2012-5-7
更多回答（1 个）
Richard Brown
2012-5-7
1 个投票
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 个评论
Richard Brown
2012-5-7
Oh, and a hint: Select all the code, then Text->Smart Indent
Image Analyst
2012-5-7
Cool. Nice trick that I didn't know about.
类别
在 帮助中心 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
2012-5-7
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
