How to terminate my code?

2 May 2023
1 Answer
Updated 6 May 2023
9 Views (30 days)

0 votes

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
Show 6 older comments Hide 6 older comments

Walter Roberson
Walter Roberson on 2 May 2023
You have the expression
(Dt*/2*DY^2
which has * followed immediately by /
Two lines above that the expression is instead
(Dt/2*DY^2
with no *
Nathi
Nathi on 2 May 2023
Edited: Nathi on 2 May 2023
Ran in:
Ran in:
@Walter Roberson My apologise!., I debbugged my code and use updated one.
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
Unrecognized function or variable 'crank'.
%k=k+1;
% end
Walter Roberson
Walter Roberson on 2 May 2023
If you need further assistance, we would need your crank function.
Nathi
Nathi on 3 May 2023
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
per isakson
per isakson on 3 May 2023
% HOW TO USE THIS iterate loop to converge 1 to 0
What should converge?
Walter Roberson
Walter Roberson on 3 May 2023
Ran in:
Ran in:
It does converge. But there is a lot of needless output that confuses matters.
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))));
error(m)
if error(m)<tol
break
end
end
1 = 0.000000 2 = 0.000000 3 = 0.000000 4 = 0.000000 5 = 0.000000 6 = 0.000000 7 = 0.000000 8 = 0.000000 9 = 0.000000 10 = 0.000000 11 = 0.000000 12 = 0.000000 13 = 0.000000 14 = 0.000000 15 = 0.000000 16 = 0.000000 17 = 0.000000 18 = 0.000000 19 = 0.000000 20 = 0.000000 21 = 0.000000 22 = 0.000000 23 = 0.000000 24 = 0.000000 25 = 0.000000 26 = 0.000000 27 = 0.000000 28 = 0.000000 29 = 0.000000 30 = 0.000000 31 = 0.000000 32 = 0.000000 33 = 0.000000 34 = 0.000000 35 = 0.000000 36 = 0.000000 37 = 0.000000 38 = 0.000000 39 = 0.000001 40 = 0.000001 41 = 0.000003 42 = 0.000006 43 = 0.000013 44 = 0.000028 45 = 0.000060 46 = 0.000129 47 = 0.000279 48 = 0.000604 49 = 0.001306 50 = 0.002825 51 = 0.006110 52 = 0.013213 53 = 0.028575 54 = 0.061795 55 = 0.133638 56 = 0.289005
ans = 0.2890
1 = 0.000000 2 = 0.000000 3 = 0.000000 4 = 0.000000 5 = 0.000000 6 = 0.000000 7 = 0.000000 8 = 0.000000 9 = 0.000000 10 = 0.000000 11 = 0.000000 12 = 0.000000 13 = 0.000000 14 = 0.000000 15 = 0.000000 16 = 0.000000 17 = 0.000000 18 = 0.000000 19 = 0.000000 20 = 0.000000 21 = 0.000000 22 = 0.000000 23 = 0.000000 24 = 0.000000 25 = 0.000000 26 = 0.000000 27 = 0.000000 28 = 0.000000 29 = 0.000000 30 = 0.000000 31 = 0.000000 32 = 0.000000 33 = 0.000000 34 = 0.000000 35 = 0.000000 36 = 0.000000 37 = 0.000000 38 = 0.000000 39 = 0.000000 40 = 0.000001 41 = 0.000002 42 = 0.000005 43 = 0.000011 44 = 0.000024 45 = 0.000053 46 = 0.000116 47 = 0.000254 48 = 0.000554 49 = 0.001210 50 = 0.002644 51 = 0.005773 52 = 0.012606 53 = 0.027524 54 = 0.060092 55 = 0.131182 56 = 0.286349
ans = 0.0027
1 = 0.000000 2 = 0.000000 3 = 0.000000 4 = 0.000000 5 = 0.000000 6 = 0.000000 7 = 0.000000 8 = 0.000000 9 = 0.000000 10 = 0.000000 11 = 0.000000 12 = 0.000000 13 = 0.000000 14 = 0.000000 15 = 0.000000 16 = 0.000000 17 = 0.000000 18 = 0.000000 19 = 0.000000 20 = 0.000000 21 = 0.000000 22 = 0.000000 23 = 0.000000 24 = 0.000000 25 = 0.000000 26 = 0.000000 27 = 0.000000 28 = 0.000000 29 = 0.000000 30 = 0.000000 31 = 0.000000 32 = 0.000000 33 = 0.000000 34 = 0.000000 35 = 0.000000 36 = 0.000000 37 = 0.000000 38 = 0.000000 39 = 0.000000 40 = 0.000001 41 = 0.000002 42 = 0.000005 43 = 0.000011 44 = 0.000025 45 = 0.000054 46 = 0.000117 47 = 0.000255 48 = 0.000557 49 = 0.001215 50 = 0.002651 51 = 0.005785 52 = 0.012624 53 = 0.027551 54 = 0.060127 55 = 0.131222 56 = 0.286380
ans = 4.0001e-05
1 = 0.000000 2 = 0.000000 3 = 0.000000 4 = 0.000000 5 = 0.000000 6 = 0.000000 7 = 0.000000 8 = 0.000000 9 = 0.000000 10 = 0.000000 11 = 0.000000 12 = 0.000000 13 = 0.000000 14 = 0.000000 15 = 0.000000 16 = 0.000000 17 = 0.000000 18 = 0.000000 19 = 0.000000 20 = 0.000000 21 = 0.000000 22 = 0.000000 23 = 0.000000 24 = 0.000000 25 = 0.000000 26 = 0.000000 27 = 0.000000 28 = 0.000000 29 = 0.000000 30 = 0.000000 31 = 0.000000 32 = 0.000000 33 = 0.000000 34 = 0.000000 35 = 0.000000 36 = 0.000000 37 = 0.000000 38 = 0.000000 39 = 0.000000 40 = 0.000001 41 = 0.000002 42 = 0.000005 43 = 0.000011 44 = 0.000025 45 = 0.000054 46 = 0.000117 47 = 0.000255 48 = 0.000556 49 = 0.001214 50 = 0.002650 51 = 0.005784 52 = 0.012624 53 = 0.027550 54 = 0.060127 55 = 0.131221 56 = 0.286380
ans = 6.9662e-07
%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
Nathi
Nathi on 3 May 2023
@Walter RobersonThank you for assit me. Actually, there are two types of convergence here. one is convergence at each time steps and another is convergence at each space steps. In my code, converged at time discretize. That's why i'm using error/tol condition. @per isakson Now I want convergence at each space step i.e., begins with 1 and increase gradually, then diminish gradually to 0 in each 'm'th row. here, my boundary condition are
if m==1
for n=1:nodey
if n==1 %(top boundary condtion i.e (1,1) 1st row and 1st column)
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 %(bottom boundary condition i.e(nodey,1) last row and 1st col)
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 %(at top boundary condition i.e (1,2:NT)
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 %(at bottom boundary condtion i.e (nodey,2:NT)
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
If this S get gradually values, then it will store in T. and then if I plot at each column wise in T, it must comes like this "begins with 1 and increase gradually, then diminish gradually to 0 ". For the gradual values and convergence, I'm using "ITERATION".
Nathi
Nathi on 6 May 2023
Edited: Nathi on 6 May 2023
@Walter Roberson and @KSSV Please try and include the below condition . It might work.
TT=crank(P,Q,R,S);
T(:,m)=[initial value,TT,endpoint];
t_prev is initial point and t_surf is end point
If get any idea about iteration with boundary condition using while loop., please suggest me.

Sign in to comment.

Sign in to answer this question.

Answers (1)

KSSV
KSSV on 3 May 2023

0 votes

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

1 Comment
Show -1 older comments Hide -1 older comments

Nathi
Nathi on 3 May 2023
@KSSV I tried your idea, it's better but not reached with my boundary condition. I'm using iteration for convergence at each space steps and gradual values with my boundary condition.
if m==1
for n=1:nodey
if n==1 %(top boundary condtion i.e (1,1) 1st row and 1st column)
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 %(bottom boundary condition i.e(nodey,1) last row and 1st col)
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 %(at top boundary condition i.e (1,2:NT)
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 %(at bottom boundary condtion i.e (nodey,2:NT)
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
this boundary condtion must satisfy in T. I need the output like If this S get gradually values, then it will store in T. and then if I plot at each column wise in T, it must comes like this "begins with 1 and increase gradually, then diminish gradually to 0 ".

Sign in to comment.

Categories

Find more on Logical in Help Center and File Exchange

Tags

Asked:

Nathi
Nathi
on 2 May 2023

Edited:

Nathi
Nathi
on 6 May 2023

Community Treasure Hunt

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

Start Hunting!