How to stop time loop when steady state is reached?
2 次查看(过去 30 天)
显示 更早的评论
I'm using unsteady case. so, it will reach a steady state at a certain time level. I fixed time at 'j'th column wise and it ran at maximum time level, but I want to stop the 'time loop' when steady state is reached. Even though I'm using a steady criteria i.e error >tolerance along dt>0 using 'while' loop, it didn't work and using dt>0, my code didn't stop.
clc; close all; clear all;
ymax=20; n=80; dy=ymax/n;
xmax=1; m=20; dx=xmax/m;
tmax=100; nt=500; dt=tmax/nt; t=0:dt:tmax;
%tol=1e-5; max_difference=1.0;
UOLD=zeros(n,nt); VOLD=zeros(n,nt);
TNEW=0; TOLD=ones(n,nt);
A=zeros([1,m]);
B=A;
C=A;
D=A;
T=TOLD;
tic
%while dt>0 && max_difference>tol
for j=1:m
for i=1:n
if j>i
C(i)=(dt*VOLD(i,j)/4*dy)-((dt)/(2*dy^2));
elseif i>j
A(i)=(-dt*VOLD(i,j)/4*dy)-((dt)/(2*dy^2));
elseif i==j
B(i)=1+(dt*UOLD(i,j)/2*dx)+((dt)/(dy^2));
end
end
end
for j=2:m
for i=1:n
D(i)=(-dt*UOLD(i,j)*((-TNEW+TOLD(i,j)-TNEW)/(2*dx)))+(((dt)/(2*dy^2))*(TOLD(i+1,j)-(2*TOLD(i,j))+TNEW))-(((dt*VOLD(i,j))/(4*dy))*(TNEW-TOLD(i+1,j)-TOLD(i,j)));
end
end
T(:,j)=TriDiag(A,B,C,D); %call tridiagonal
dt=0.2+dt;
Thanks for any advice.
3 个评论
Mathieu NOE
2023-4-20
ok
I am not sure to understand all the details of your code , but something disturbs me
if we consider that the metric (error) that defines how your computation converge to the supposed steady state values
difference = (abs(T(:,j)-T(:,j-1))./max(1,abs(T(:,j-1))));
then if I plot it for the 80 iterations , we see the difference increasing and not decreasing .... so I wonder if you have 100% confidence in your code
ploted with y log scale to make it clear :
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1362188/image.png)
Image Analyst
2023-4-20
Which loop, the one over i or the one over j, is the "time loop"?
I don't see any line like "error >tolerance". Where exactly is that line?
采纳的回答
Mathieu NOE
2023-4-21
hello
think I have understood
your difference computation must be indexed with j , so it must be inside the for loop
I used a break to exist the for loop when the difference is below tol
now, tol = 1e-5 will not be achieved with only 500 iteration
for the sake of the demo I put tol = 1e-2 and it will do only 103 iterations
plot of difference vs j with tol = 1e-5
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1362693/image.png)
plot of difference vs j with tol = 1e-2
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1362698/image.png)
code :
clc; close all; clear all;
ymax=20; n=80; dy=ymax/n;
xmax=1; m=20; dx=xmax/m;
tmax=100; nt=500; dt=tmax/nt; t=0:dt:tmax;
tol=1e-2;
UOLD=zeros(n,nt); VOLD=zeros(n,nt);
TNEW=0; TOLD=TNEW*ones(n,nt); TWALL=ones(1,length(t));
A=zeros([1,m]);
B=A;
C=A;
D=A;
T=TOLD;
difference(1) = 1;
tic
for j=1:nt
for i=1:n
if j>i
C(i)=(dt*VOLD(i,j)/4*dy)-((dt)/(2*dy^2));
elseif i>j
A(i)=(-dt*VOLD(i,j)/4*dy)-((dt)/(2*dy^2));
elseif i==j
B(i)=1+(dt*UOLD(i,j)/2*dx)+((dt)/(dy^2));
end
end
end
for j=2:nt
if j==1
for i=1:n
if i==1
D(i)=(-dt*UOLD(i,j)*((-TNEW+TOLD(i,j)-TNEW)/(2*dx)))+(((dt)/(2*dy^2))*(TOLD(i+1,j)-(2*TOLD(i,j))+TNEW))-(((dt*VOLD(i,j))/(4*dy))*(TNEW-TOLD(i+1,j)-TOLD(i,j)))-((dt/(4*dy))*VOLD(i,j))-((dt/2*dy^2)*TNEW);
elseif i==n
D(i)=(-dt*UOLD(i,j)*((-TNEW+TOLD(i,j)-TNEW)/(2*dx)))+(((dt)/(2*dy^2))*(TWALL(j)-(2*TOLD(i,j))+TOLD(i-1,j)))-(((dt*VOLD(i,j))/(4*dy))*(TOLD(i-1,j)-TWALL(j)+TOLD(i,j)))-(((-dt)/(4*dy))*VOLD(i,j))-((dt/2*dy^2)*TWALL);
else
D(i)=(-dt*UOLD(i,j)*((-TNEW+TOLD(i,j)-TNEW)/(2*dx)))+(((dt)/(2*dy^2))*(TOLD(i+1,j)-(2*TOLD(i,j))+TOLD(i-1,j)))-(((dt*VOLD(i,j))/(4*dy))*(TOLD(i-1,j)-TOLD(i+1,j)+TOLD(i,j)));
end
end
else
for i=1:n
if i==1
D(i)=(-dt*UOLD(i,j)*((-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1))/(2*dx)))+(((dt)/(2*dy^2))*(TOLD(i+1,j)-(2*TOLD(i,j))+TNEW))-(((dt*VOLD(i,j))/(4*dy))*(TNEW-TOLD(i+1,j)+TOLD(i,j)))-((dt/(4*dy))*VOLD(i,j))-((dt/2*dy^2)*TNEW);
elseif i==n
D(i)=(-dt*UOLD(i,j)*((-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1))/(2*dx)))+(((dt)/(2*dy^2))*(TWALL(j)-(2*TOLD(i,j))+TOLD(i-1,j)))-(((dt*VOLD(i,j))/(4*dy))*(TOLD(i-1,j)-TWALL(j)+TOLD(i,j)))-(((-dt)/(4*dy))*VOLD(i,j))-((dt/2*dy^2)*TWALL(j));
else
D(i)=(-dt*UOLD(i,j)*((-T(i,j-1)+TOLD(i,j)-TOLD(i,j-1))/(2*dx)))+(((dt)/(2*dy^2))*(TOLD(i+1,j)-(2*TOLD(i,j))+TOLD(i-1,j)))-(((dt*VOLD(i,j))/(4*dy))*(TOLD(i-1,j)-TOLD(i+1,j)+TOLD(i,j)));
end
end
end
T(:,j)=TriDiag(A,B,C,D); %call tridiagonal
dt=0.2+dt;
TOLD=T;
difference(j) = max(abs(T(:,j)-T(:,j-1))./max(1,abs(T(:,j-1))));
if difference(j)<tol
break
end
end
4 个评论
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!