What is the steady state criteria in crank nicolson scheme, especially impilicit finite difference method? and How to stop time level when steady state is reached?

6 次查看（过去 30 天）
Yanni
Yanni 2023-3-27
编辑： Yanni 2023-4-11
Ran in:
I'm dealing with unsteady case in my code. so, the CN term updated in each time step and it ran upto final time (tmax=100). but i want to stop the time level when steady state is reached. how to include steady state criteria in my code? and what is the steady state criteria?
clc; clear all; close all
xmax=1; ymax=7; tmax=10; m=20; n=28; nt=500;
dx=xmax/m; dy=ymax/n; dt=tmax/nt; y=dy:dy:ymax;
%====================INTIALIZATION===============================%
UOLD=zeros(m,n); VOLD=zeros(m,n);
CNEW=0;COLD=ones(m,n);
tic
CN=COLD;
for j=1:n
for i=1:m
C(i)=((dt*VOLD(i,j))/(4*dy))-(dt/(2*(dy^2)));
A(i)=((-dt*VOLD(i,j))/(4*dy))-(dt/(2*(dy^2)));
B(i)=1+((dt*UOLD(i,j))/(2*dx))+(dt/(dy^2));
end
end
for j=2:n
for i=1:m
D(i)=(-dt*UOLD(i,j)*((-CNEW+COLD(i,j)-CNEW)/(2*dx)))+((dt*(1/(2*dy^2)))*(COLD(i+1,j)-(2*COLD(i,j))+CNEW))-(((dt*VOLD(i,j))/(4*dy))*(CNEW-COLD(i+1,j)-COLD(i,j)));
end
end
Index in position 1 exceeds array bounds. Index must not exceed 20.
dt=0.2+dt;
CN(i,j)=TriDiag(A,B,C,D);
toc

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2023-3-27
移动：Torsten 2023-3-27
You can define several stopping criteria.
Try
max(abs(1-CN(:)./COLD(:))) < tol
with tol = 1e-4 ... 1e-6
  4 个评论
显示 2更早的评论
Torsten
Torsten 2023-3-28
编辑：Torsten 2023-3-28
"dt" is used as the time increment in your code, not as the absolute time.
Thus setting dt = 0.2+dt and thus making the time increment bigger and bigger makes no sense.
If your equations do not explicitly depend on time, you will only see the time increment "dt" in their discretized form, not the absolute time "t" itself.
Torsten
Torsten 2023-3-29
In your double loops over i and j, the arrays A(i), B(i), C(i) and D(i) are permanently overwritten when j changes. So in the end, the A(i), B(i), C(i) and D(i) are the results for j=n before you call TriDiag. So somehow you will have to work with A(i,j), B(i,j), C(i,j) and D(i,j) to capture the j dependence of the coefficients.
Then I guess the CN and COLD arrays should depend on i, j and time (thus they should be three-dimensional instead of two-dimensional matrices), shouldn't they ?
Further, the j in my piece of code refers to time. You will have to choose another index apart from i and j, maybe itime.
j = 1;
while t < tmax && max_difference > tol
j = j + 1;
t = t + dt;
% Compute solution C for new time
...
% Save solution in matrix C
C(:,j) = solution at time step j;
% Compute maximum difference of the solution between the last two
% output times
max_difference = max ( abs(C(:,j)-C(:,j-1))./max(1,abs(C(:,j-1))) );
end
Maybe you should include the partial differential equation you are trying to solve in a mathematical notation so that I get an idea what problem you are trying to discretize here.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by