How to terminate my code?
2 views (last 30 days)
Show older comments
clc;
close all;
clear all;
%-----------SPACEGRID---------%
Y=14; nodey=56; DY=Y/nodey; DX=0.05;
%--------TIMEGRID-------------%
TMAX=100; NT=500; Dt=TMAX/NT; t=0:Dt:TMAX;
error(1)=1; tol=1e-5;
MI=100;
u_prev=zeros(nodey,NT); v_prev=zeros(nodey,NT);
t_cur=0; t_prev=t_cur*ones(nodey,NT); T_surf=ones(1,length(t));
P=zeros([1,nodey]); Q=P; R=P; S=P;
T=t_prev;
k=0;
%while k<MI %HOW TO USE THIS iterate loop to converge 1 to 0
for m=1:NT
for n=1:nodey
if m>n
R(n)=(Dt*v_prev(n,m)/4*DY)-(Dt/(2*DY^2));
elseif n>m
P(n)=(-Dt*v_prev(n,m)/4*DY)-(Dt/(2*DY^2));
elseif n==m
Q(n)=1+(Dt*u_prev(n,m)/2*DX)+(Dt/(DY^2));
end
end
end
for m=2:NT
if m==1
for n=1:nodey
if n==1
S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)-t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
elseif n==nodey
S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf);
else
S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
end
end
else
for n=1:nodey
if n==1
S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)+t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
elseif n==nodey
S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf(m));
else
S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
end
end
end
T(:,m)=crank(P,Q,R,S);
Dt=0.2+Dt;
t_prev=T;
error(m) = max(abs(T(:,m)-T(:,m-1))./max(1,abs(T(:,m-1))));
if error(m)<tol
break
end
end
%k=k+1;
% end
function x = crank(e,f,g,r)
%input
% e = subdiagonal vector
% f = diagonal vector
% g = superdiagonal vector
% r = right hand side vector
% output:
% x = solution vector
n=length(f);
for k = 2:n
factor = e(k)/f(k-1);
f(k) = f(k) - factor*g(k-1);
r(k) = r(k) - factor*r(k-1);
end
x(n) = r(n)/f(n);
for k = n-1:-1:1
x(k) = (r(k)-g(k)*x(k+1))/f(k);
end
for i = 1:length(x)
fprintf('\n%d = %f\n', i, x(i));
end
8 Comments
Answers (1)
KSSV
on 3 May 2023
You can use for loop instead of while....You know the number of counts right.
clc;
close all;
clear all;
%-----------SPACEGRID---------%
Y=14; nodey=56; DY=Y/nodey; DX=0.05;
%--------TIMEGRID-------------%
TMAX=100; NT=500; Dt=TMAX/NT; t=0:Dt:TMAX;
error(1)=1; tol=1e-5;
MI=100;
u_prev=zeros(nodey,NT); v_prev=zeros(nodey,NT);
t_cur=0; t_prev=t_cur*ones(nodey,NT); T_surf=ones(1,length(t));
P=zeros([1,nodey]); Q=P; R=P; S=P;
T=t_prev;
E = cell(MI,1) ;
for k = 1:MI %HOW TO USE THIS iterate loop to converge 1 to 0
for m=1:NT
for n=1:nodey
if m>n
R(n)=(Dt*v_prev(n,m)/4*DY)-(Dt/(2*DY^2));
elseif n>m
P(n)=(-Dt*v_prev(n,m)/4*DY)-(Dt/(2*DY^2));
elseif n==m
Q(n)=1+(Dt*u_prev(n,m)/2*DX)+(Dt/(DY^2));
end
end
end
for m=2:NT
if m==1
for n=1:nodey
if n==1
S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)-t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
elseif n==nodey
S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf);
else
S(n)=(-Dt*u_prev(n,m)*-t_cur+t_prev(n,m)-t_cur/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
end
end
else
for n=1:nodey
if n==1
S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_cur)-(Dt*v_prev(n,m)/4*DY*t_cur-t_prev(n+1,m)+t_prev(n,m))-(Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*t_cur);
elseif n==nodey
S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*T_surf(m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-T_surf(m)+t_prev(n,m))-(-Dt/4*DY*v_prev(n,m))-(Dt/2*DY^2*T_surf(m));
else
S(n)=(-Dt*u_prev(n,m)*-T(n,m-1)+t_prev(n,m)-t_prev(n,m-1)/2*DX)+(Dt/2*DY^2*t_prev(n+1,m)-2*t_prev(n,m)+t_prev(n-1,m))-(Dt*v_prev(n,m)/4*DY*t_prev(n-1,m)-t_prev(n+1,m)+t_prev(n,m));
end
end
end
T(:,m)=crank(P,Q,R,S);
Dt=0.2+Dt;
t_prev=T;
error(m) = max(abs(T(:,m)-T(:,m-1))./max(1,abs(T(:,m-1))));
if error(m)<tol
break
end
E{k} = error ;
end
%k=k+1;
end
function x = crank(e,f,g,r)
%input
% e = subdiagonal vector
% f = diagonal vector
% g = superdiagonal vector
% r = right hand side vector
% output:
% x = solution vector
n=length(f);
for k = 2:n
factor = e(k)/f(k-1);
f(k) = f(k) - factor*g(k-1);
r(k) = r(k) - factor*r(k-1);
end
x(n) = r(n)/f(n);
for k = n-1:-1:1
x(k) = (r(k)-g(k)*x(k+1))/f(k);
end
for i = 1:length(x)
fprintf('\n%d = %f\n', i, x(i));
end
end
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!