Effect of time step size on the numerical solution.

• I am going to write a function in MATLAB that accepts time step as an argument and solves the same problem as mentioned here and compare the effect of time step on the numerical solution.

Problem Setup - Addition of assumptions to the problem setup discussed above:
1. Grid points set to `80`.
2. Time steps to be used are `1e-1`, `1e-2`, `1e-3` and `1e-4`.

PROGRAM

Function file:

The function file 'convection.m' looks like this:

function [u,x] = convection(dt)  

%%Writing down the assumptions and applying Boundary conditions to initialize the velocity profile
L = 1;  % Length of the domain
C = 1;  % Constant speed  for convection
T = 0.4;  % Total time
nt = T/dt + 1;  % Setting up iterations on basis of dt
n = 81;  %Setting grid points to 80
x = linspace(0,L,n);  % Creating 1D grid in x direction
dx = L/(n-1);  % Difference between two consecutive grid points
u = ones(1,n);  % Creating velocity array with all velocity values = 1 including BC
%Using the formula for nth term of AP, I derived the following two results
n_start = 0.1/dx + 1;  % Starting node of the square wave
n_stop = 0.3/dx + 1;  % Ending node of the square wave
%Using the above results obtained
u(n_start:n_stop) = 2;  % Changing velocity values from 1 to 2 for the nodes where square wave will start and end
u_old = u;  %Previous node velocity


%%Time loop - Time marching solution
for i = 1:nt
     %Space loop - Space marching
     for j = 2:n
          u(j) = u_old(j) - (C*dt/dx)*(u_old(j) - u_old(j-1));
     end
     %Updating velocity
     u_old = u;
end     

CFL = dt/dx;  %CFL = C*dt/dx it can help us knowing the stability.

fprintf('The CFL number for time step size of %g is %g. ', dt, CFL)

if CFL <= 1 
   disp('Thus, stable.')
else 
   disp('Thus, unstable.')
end   

Driver Script:

The driver script 'effect_of_time_step_size_on_sol_time.m'looks like this:

% Effect of time step size on solution time

clear all
close all
clc



tic
[u0,x0] = convection(1e-1);      %Time step size of 1e-1
et1 = toc;

tic
[u1,x1] = convection(1e-2);      %Time step size of 1e-2
et2 = toc;

tic
[u2,x2] = convection(1e-3);      %Time step size of 1e-3
et3 = toc;

tic
[u3,x3] = convection(1e-4);     %Time step size of 1e-4
et4 = toc;

y = [et1;et2;et3;et4];

%%Generating plots and saving them
%Compare the numerical solution at t = 0.4 seconds for all of the time steps used above
figure(1)
plot(x0, u0, '-bp','linewidth',2,'markerfacecolor','w')
hold on
plot(x1, u1, '--ko','linewidth',2,'markerfacecolor','w')
plot(x2, u2, ':ys','linewidth',2,'markerfacecolor','w')
plot(x3, u3, '-.rd','linewidth',2,'markerfacecolor','w')
set(gca,'XTick',(0:0.1:1))
xlabel('1D Grid','FontSize',12)
ylabel('Vlocity Profile','FontSize',12)
legend('1e-1','1e-2','1e-3','1e-4')
title('Comparing the numerical solution at T = 0.4 seconds','FontSize',15)
grid on
axis([0 1 0.8 2.2])
saveas(gca,'test.png')
%Plotting the simulation time as a function of the time step
figure(2)       
bar(y,0.4,'k','EdgeColor','b','linewidth',0.8);
j = legend('Solution time for 80 grid points','location', 'northwest');
set(j,'FontSize',12);
ylabel('Solution Time','FontSize',12)
xlabel('Time Steps')
title('Comparing Solution Time For Different Time Steps','FontSize',15)
labels = {sprintf('%0.3f',et1),sprintf('%0.3f',et2),sprintf('%0.3f',et3),sprintf('%0.3f',et4)};
xt = get(gca,'XTick');
text(xt,y, labels, 'HorizontalAlignment','center', 'VerticalAlignment','bottom','FontSize',10)
set(gca,'xtickLabel',{'1e-1', '1e-2','1e-3','1e-4'},'FontSize',12)
saveas(gca,'bar.png')

Results
Bar graph generated:

Velocity profiles for different time steps used above are plotted together for T = 0.4 seconds:

The stability of the velocity profile based on CFL criteria displayed on the command window:

The CFL number for time step size of 0.1 is 8. Thus, unstable.
The CFL number for time step size of 0.01 is 0.8. Thus, stable.
The CFL number for time step size of 0.001 is 0.08. Thus, stable.
The CFL number for time step size of 0.0001 is 0.008. Thus, stable.

Explaining the results:
1. The least solution time with the best velocity profile obtained was for time step size of `1e-2`as CFL number for respective time step size is closest to 1. The shape of the profile is lost due to the artificial diffusion.
2. As CFL number for time step size of `1e-1 > 1` the solution goes unstable after in a short time.
3. The velocity profiles for time step sizes of `1e-3` and `1e-4` respectively are pretty much similar though it is observed that the solution time for `1e-4 > 1e-3`