week-5 :: Stability analysis in unsteady problem..

Objective :: Detail analysis of Stability analysis methods for numerical methods.

1. Types of Stability Analysis :

• The discrete Pertubation stability analysis.
• The Matrix method stability analysis.
• The Von-Neumann stability analysis.

2. Importance to understand stability analysis :

• As known , for convergence of the solution consistency and stability conditions should satisfy.
• As stability depends on the Finite precesion of the computation , so we made a certain condition , we can may conlusions stability in solution satisfy.
• Find the CFL number for different schemes for different physical problems , as domain of dependence varies with type of schemes.

3. Von-Neumann staility analysis :

• Limitations ::

1. Linear PDE.
2. Constant Co-efficinets.
3. Assume Periodic Boundary conditions (Effect of boundary conditions is neglected).

`phi^(n+1) = A*phi^(n)`

`phi` = Notation of error at time step (n+1).

`phi` = Notation of error at time step (n).

A = Amplification factor (or) growth term (or) decay term.

4. Amplification factor tell about the numerical method.

• `absA` < 1 , Stable and Damping.
• `absA` = 1 , Stable.
• `abs A`> 1 , Unstable.

5. Types of Errors :

• Discretization error (De):

Difference between Exact solution of pde(E) and Exact solution of FDE(D).

• Round off Error (`epsilon`):

Difference between Exact solution of FDE and Numerical solution of FDE(N).

• `epsilon = D - N``

6. Fourier series :

Since we considered the periodic boundary conditions , a periodic domain. Can represent the

and error in terms of fourier series.

`phi = sum_(k = -prop)^(prop) A_k * e^(i*k*x)`

k = wave number.

Over increasing the wave number the approximate solution converges.

Here we consider only one wave number , as the pde is linear , superposition of solutions of all

wavenumber also be a solution.

So , `phi_(j)^n = A^n * e^(i*k*j*Deltax)`

Uniform grid : x = `j*Deltax`.

7. Linear 1-D convection equation for FTBS :

Dicretization equation of advection equation -> ```phi_j^(n+1) = phi_j^n - c*(phi_j^(n) - phi_(j-1)^n)`

c = courant number (`frac{c*Deltax}{Deltat}`).

Substitute in `phi_i^n = A^n * e^(i*kj*Deltax)`

By simple mathematical manipulation ,

A = 1 - c*(1 - `e^(-i*k*Deltax)`).

`e^(-i*k*Deltax)` = `cos(k*Deltax) -i*sin(k*Deltax)`.

`implies A = 1 - c*(1-cos(k*Deltax)) - i*c*sin(k*Deltax).`

From complex analysis , `absA^2 = A*A^(**)`

`absA^2 = 1-2*c*(1-cos(k*Deltax)) + c^2*(1-2*cos(k*Deltax) + cos^2(k*Deltax)) + c^2*sin^2(k*Deltax).`

`implies absA^2 = 1-2*c*(1-c)*(1-cos(k*Deltax)).`

Condition , `absA<=1.``iff absA^2 - 1<=0`

`iff -2*c*(1-c)*(1-cos(k*Deltax)) <=0.`

`iff c*(1-c)*(1-cos(k*Deltax)>=0.`

`implies c*(1-c)>=0 iff 0<=c<=1 (because 1-cos(k*Deltax > 0))`

`so , 0<=c<=1` , as domain of dependence of numerical solution present inside the domain of dependence of physical(domain)`

(CTCS) , `(-1<=c<=1)`

Code :::

%% Clearing commands
close all
clear all
clc
% Input variables for the linear convection 
L =1;
n = [20 , 40 , 80 ,160];
%n = 20;
c = 1;
dt = 1e-2;
t = 0.4;
nt = t/dt;
% Calculations
[q,u,dx] = linear_convection_func1(n,L);
 u_initial = u;
 for k = 1:length(n)
 figure(k)
 v = u{1,k};
 cfl = (c*dt)/dx(k);
 % Now for different time steps for velocity  
   % For loop for number of time steps
   for i = 1:nt  
   % For loop for all grid velocity values in the mesh
   for j = 2:length(q{1,k}) 
       
       z(j) = v(1,j) - (c*dt)/dx(k)*(v(1,j) - v(1,j-1)); 
     % Tried relaxation not pretty sure...
       %{
       if cfl<1
       cfl = 1.10*cfl
       end
       if cfl>1
          cfl = cfl - 0.50*cfl 
      end
       %}
   end
       z(1) = 1;
       v = z;
       plot(q{1,k},u{1,k},'color','b','linewidth',3)
       hold on
       plot(q{1,k},z,'r')
       pause(0.3)
       axis([0 1 0 2])
   end
   hold off
 end

Description about code snippet :

1. Used clear command for to clear the data in the present interface.

2. Declared the required variables.

3. Used function defination , for to get data to plot the initial velocity profile for different number of grids considered.

4. Used for loops , to calculate velocity data and plot the results at the simulation time itslef using 'pause' command.

5. axis command , to set the axis bounds.

6. 'hold on' command used , not to clear axes properties for the corresponding future plots.

Function :::

function [q,v,dx] = linear_convection_func1(n,L)
  
for i = 1:length(n)
   q{1,i} = linspace(0,L,n(i));  % Making mesh
   x = q{1,i};
   dx(i) = x(2) - x(1);     % Grid spacing
   u = 1*ones(1,n(i));    % Initial conditions
   A = logical(x>=0.1 & x<=0.3);
   index = find(A==1);
   n_start = index(1);
   n_end = index(end);
   u(n_start:n_end) = 2;  
   v{1,i} = u;
end
  
end
%{
1. find [start,stop] values.
2. U values;
%}

Description of function code snippet :

1. Used function , to find the mesh discretization , grid spacing , velocity profiles data , start and end indices for grid locations (0.1 and 0.3).

2. Using simple steps and properties in matlab.

Results :

1. The first polt , for n = 20. (number of grids)

2. Fixed Time step (dt) = 1e-2.

3. Number of iterations , nt = t/dt = 40 time steps for final velocity profile.

-> dx = 0.006289308176101

-> dt = 1e-2.

-> Co = 1.5900 > 1 , Local truncation error increases instability occurs.

This condition is for courant number is 'one' , so No truncation erroe -> No diffusion...

For implicit method , let's consider the 2-D heat equation :

For heat equation unsteady -> CFL <= 0.5.

-> CFL < 0.5.

%%                   ****** Unsteady using Jacobi method ******

close all
clear all
clc

nx = 20;
ny = 20;
nt = 1000;

x = linspace(0,1,nx);
y = linspace(0,1,ny);

dx = x(2) - x(1);
dy = y(2) - y(1);

y1 = 1;

error = 9e9;
tolerance = 1e-4;
dt = 6e-4;
alpha = 1.4;

T_L = 400;
T_R = 800;
T_T = 600;
T_B = 900;

T = 300*ones(nx,ny);
T(2:ny-1,1) = T_L;
T(2:ny-1,nx) = T_R;
T(1,2:nx-1) = T_T;
T(ny,2:ny-1) = T_B;

T(1,1) = (T_L + T_T)/2;
T(nx,ny) = (T_R + T_B)/2;
T(1,ny) = (T_T + T_R)/2;
T(nx,1) = (T_L + T_B)/2;

CFL = ((alpha*dt)/(dx^2)) + ((alpha*dt)/(dy^2))
k1   =  (1 + ((2*alpha*dt)/(dx^2)) + ((2*alpha*dt)/(dy^2)));
k2   =  ((alpha*dt)/(dx^2));
k3   =  ((alpha*dt)/(dy^2));
Tp = T;
Tn =T;

[X,Y] = meshgrid(x,y);
choose = input('enter 1 for implcit otherwise any number :');

for i = 1:1000% Time loop
    error = 9e9;  
   while error > 1e-4
      for j = 2:nx-1
         for k = 2:ny-1  
             term1 = k1^(-1);
             term2 = k2*term1;
             term3 = k3*term1;
             if choose ==1 
                Tn(j,k) = (term1*(Tp(j,k))) + (term2*(T(j+1,k) + T(j-1,k))) + (term3*(T(j,k+1) + T(j,k-1)));
                title1 = sprintf('implicit');

             else
                Tn(j,k) =  ((1 - 2*k2 - 2*k3)*T(j,k)) + (k2*(T(j+1,k) + T(j-1,k))) +  (k3*(T(j,k+1) + T(j,k-1))); 
                title1 = sprintf('explicit');
             end
         end
      end
      error = max(max(abs(Tn - T)));
      T = Tn;
   end
   Tp = Tn;
   contourf(X,-Y,T)
   title2 = sprintf('Transient state using %s jacobi method CFL = %d' ,title1,CFL);
   title(title2);
   pause(0.000000001)
end

Results :

As can observe the results for CFL number > 0.5 , since CFL < 0.5 for 2-D heat conduction equation.

For implicit , no effect of CFL number , CFL gives the stability information only for explicit.

Implicit methods are always stable for any CFL number.