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?

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

Answer 1

[Torsten]

You can define several stopping criteria.

Try

max(abs(1-CN(:)./COLD(:))) < tol

with tol = 1e-4 ... 1e-6

[Torsten]

I don't understand you code, so I cannot help you in this respect.

Surely an error is the line

dt=0.2+dt; 

It should read

t = t + dt;

because you have to increment t, not dt.

Usually, what I described can be done like

t = 0;
tol = 1e-6;
max_difference = 1.0;
C(:,1) = initial values of solution;
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

[Torsten]

@Gayathri comment moved here:

I have doubt in your suggested code. In my equation, the variable 'dt' for time level. but you taking for increment time level is the variable 't' and I'm not taking the variable 't' in my equation as well. then How it will vary for each time level?. that's why my code also not varying in each grid when I use t=t+dt;. moreover, after i saw you suggested code, I understood steady state reached criteria better., neverthless, i want show steady state time level in command window as well. kindly advice me this content too.

[Torsten]

"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]

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.