Effect of time step size on the solution

LINEAR ONE-DIMENSIONAL CONVECTION EQUATION:

`(partialu)/(partialt)+a*(partialu)/(partialx) = 0`

where u(x,t) is the unknown function of (x,t) and a the convection speed.we can rewrite it as `u_t + a*u_x = 0`

this equation has to be substantiated by the following initial and boundary conditions, for a>0:

At t = 0 `u(x,0) = u^(0)*(x)` x [0,L]

At x = 0 `u(0,t) = g(t)` t > 0

In order to apply a finite difference method to this equation, we discretize the space and time domains, with constant steps. That is, the x-axis is discretized with N constant mesh intervals `trianglex,`and the time axis is subdivided in constant time intervals `triangle t`

`x^i = i*trianglex` `t^n = n*trianglet`

`u_i^n = u(itrianglex, ntrianglet)`

We indicate the time level n by a superscript and the space position is indicated by the subscript i.

In order to obtain a numerical scheme, we have to discretize separately the space and the time derivatives.

First order in time, first order in space

With a first order approximation for the space derivative, would be obtained with a backward difference in space, leading to the semi-discrete form:

`(u_t)_i= -a/(trianglex)*(u_t - u_(t-1))`

The left-hand side represents the time derivatives evaluated at point i, and the next step is to define a discretization in time. This implies the replacement of the time derivative by a discrete form but also a decision as to the time level at which the right-hand side will be evaluated.

In a convection problem, we know the solution at t=0 and we look for the solution at later times. This is translated numerically, when we progress in time from time level n to time level (n+1). It is therefore logical to select a forward difference formula for`(u_t)_i` leading to

`(u_i^(n+1) - u_i^n)/(trianglet) + (a/(trianglex))*(u_i - u_(i-1) ) = 0`

This simplest scheme would be obtained with an evaluation of the right-hand side of equation at time step n, for which all quantities are considered as known. This method is known as the Euler method for the time integration of ordinary differential equations, leading to the Explicit numerical scheme:

`(u_i^(n+1) - u_i^n)/(trianglet) + (a/(trianglex))*(u_i^n - u_(i-1)^n ) = 0`

These schemes are the first order in space and in time and are known as the First Order Upwind Schemes for the convection equation.

MATLAB CODE:

%% Effect of time step size on the 1D Linear Convection Equation

% Initializing the variable
dt = [1e-4 1e-3 1e-2 1e-1];

% Solving the 1D LC and recording time taken 
tic 
[x1,u1] = time_step_LC(dt(1));
t1 = toc;

tic 
[x2,u2] = time_step_LC(dt(2));
t2 = toc;

tic 
[x3,u3] = time_step_LC(dt(3));
t3 = toc;

tic 
[x4,u4] = time_step_LC(dt(4));
t4 = toc;

simulation_time = [t1 t2 t3 t4];

% Plotting
figure(1),clf
plot(x1,u1,'r','linewidth',2)
hold on 
plot(x2,u2,'g','linewidth',2)
plot(x3,u3,'b','linewidth',2)
plot(x4,u4,'m')
axis([0 1 0.5 2.5])
grid on 
xlabel('X (m)')
ylabel('Velocity (m/s)')
title('Comparison of effect of time step on the solution')
legend('1e-4','1e-3','1e-2','1e-1')

% Plotting Time Step v/s Simulation Time
figure(2),clf
h = categorical({'1e-4','1e-3','1e-2','1e-1'});
bar(h,simulation_time)
ax = gca;
ax.FontSize = 12;
grid on
xlabel('Time Step')
ylabel('Simulation Time')
title('Time-Step v/s Simulation-Time')
legend('Simulation Time for 80 grid points')

%%

FUNCTION CODE:

%% 1D - Linear Convection Equation for 80 grid points

function [x,u] = time_step_LC(dt)

% Inputs
L = 1;     % Domain Length 
n = 80;    % No. of grid points
C = 1;     % Constant Convection Velocity
T = 0.4;   % Total Time 

% Calculations
x = linspace(0,L,n);
dx = L/(n-1);

% Initialize Arrays 
u = ones(1,n);

% Boundary Conditions 
u(1) = 1;

% Start and End points for SQUARE WAVE
for i = 1:n
    if ( (0.1 <= x(i)) && (x(i) <= 0.3) )
        u(i) = 2;       
    else 
        u(i) = 1;
    end 
end 

u_old = u;

nt = T/dt;  % No. of Time Steps

%  Time Loop ( First Order Forward Difference )
for k = 1:nt
    
    %  Space Loop ( First Order Backward Difference )
    for i = 2:n
        u(i) = u_old(i) - ( C*(dt/dx) )*( u_old(i) - u_old(i-1) );
    end
    
    % Updating Old Velocities
    u_old = u;
    
end

GRAPHS:

As the time step increases, the simulation time increases. This is due to increase in the no. of iterations.

The graph shows the effect of time step on a square wave for a given no. of grid points.

For 1e-4, 1e-3, 1e-2 the square wave is stable, but for 1e-1 the square wave blows up, this is because the time step exceeds the stability condition.

The maximum time step upto which the equation can be stable can be calculated by CFL stability condition.

It is given by `C = u*(triangle t)/(triangle x) <= 1`

`C = 1, u =1 m/s`

`triangle x = L/(n-1)`, where L = 1m, n = 80; we get `triangle t = 0.01265`

After which the square wave blows up.

The 'magenta line' shows the maximum time step line.