1D Linear Convection in Matlab using first-order discretization schemes and the effect of time-step on the stability of the solution

INTRODUCTION: To solve the 1D linear equation to observe the propagation of wave with time.

Governing Equation

`{partial u}/{partial t} + C{partial u}/{partial x} = 0`

C = Velocity of Propagation

The above equation is discretized using first order forward differencing scheme for time derivative and backward differencing scheme for space derivative.

The explicit numerical solution for the above equation

`u_i^(n+1) = u_i^n - {C ∆t}/{∆x}(u_i^n- u_(i-1)^n)`

where 'n' is time step

'i' is node no.

To define a domain by using some inputs

L = 1 m

n = 80

C =1

set time step, dt = [1e-4, 1e-3, 1e-2, 1e-1]

u(t=0)= 2 m/s if `0.1m lexge0.3`

u(t=0)= 2 m/s if otherwise

The Matlab code for linear convection:

The above numerical solution into Matlab, The programme outputs the grid, the initial velocity profile and final velocity profile at t = 0.4 seconds. The code is shown below.

clear all
close all
clc

% Inputs for the linear convection problem
       % Length of the domain in 'm'    
           L = 1;
       % Number of grid points
           n = 80;
       % Linear convection constant
           C = 1
       % Time step in seconds
          dt = [1e-4 1e-3 1e-2 1e-1]
       %  End time
          t_end = 0.4;
% calculations
       % x array along the length of the domain 
          x = linspace(0,L,n);
       % Grid spacing, dx = L/(n-1) or dx =x(2)-x(1)
          dx = x(2) - x(1);
       % time 
             nt1 = t_end/dt(1); nt2 = t_end/dt(2)
             nt3 = t_end/dt(3); nt4 = t_end/dt(4)  
       % Define the start and end positions for the wave
              x_start = 0.1;    x_end = 0.3;
       % n_start = f(x_start)
              n_start = 0.1/dx +1
       % n_end = f(x_end)
              n_end  = 0.3/dx +1
              
        % Initialize arrays and apply boundary conditions
                u = ones(1,n); u(n_start:n_end) = 2  
        % u(1) = 1;
                u_old = u;    u_initial = u;

        % Time of Execution
                 t_exe = zeros(1,4)
tic
% Time loop(for time step 1e-4, and t_end = 0.4 )
   for k = 1:nt1
         % space loop
             for i = 2:n
                 u(i) = u_old(i) - (C*dt(1)/dx)*(u_old(i)- u_old(i-1));
             end
         % update old velocities
             u_old = u;   t_exe(1)=toc
   end

         % Initialize arrays and apply boundary conditions
              u1 = ones(1,n);    u1(n_start:n_end) = 2  
         % u1(1) = 1;
              u1_old = u1;     u1_initial = u1;
tic
% Time loop(for time step 1e-3, and t_end = 0.4 )
   for k = 1:nt2
        % space loop
        for i = 2:n
            u1(i) = u1_old(i) - (C*dt(2)/dx)*(u1_old(i)- u1_old(i-1));
        end
        % update old velocities
             u1_old = u1;   t_exe(2)=toc
   end

        % Initialize arrays and apply boundary conditions
              u2 = ones(1,n);     u2(n_start:n_end) = 2  
        % u2(1) = 1;
             u2_old = u2;      u2_initial = u2;
tic
% Time loop(for time step 1e-2, and t_end = 0.4 )
   for k = 1:nt3
        % space loop
          for i = 2:n
             u2(i) = u2_old(i) - (C*dt(3)/dx)*(u2_old(i)- u2_old(i-1));
          end
          % update old velocities
             u2_old = u2;        t_exe(3)=toc
   end


          % Initialize arrays and apply boundary conditions
              u3 = ones(1,n);    u3(n_start:n_end) = 2  
          % u3(1) = 1;
              u3_old = u3;    u3_initial = u3;
tic
% Time loop(for time step 1e-1, and t_end = 0.4 )
    for k = 1:nt4
        % space loop
          for i = 2:n
            u3(i) = u3_old(i) - (C*dt(4)/dx)*(u3_old(i)- u3_old(i-1));
          end
         % update old velocities
            u3_old = u3;          t_exe(4)=toc
          end
% plotting the velocities
figure(1)
plot(x,u,x,u1,x,u2,x,u3,'Linewidth',2)
grid on
axis([0 1 0 3])
xlabel('Domain Length')
ylabel('velocity')
lgd = legend(['dt=', num2str(dt(1))], ['dt=', num2str(dt(2))], ['dt=', num2str(dt(3))],['dt=', num2str(dt(4))]);
title('comparion of outputs for different time steps for T =0.4')

figure(2)
z = categorical({'1e-4', '1e-3', '1e-2', '1e-1'})
bar(z,t_exe)
grid on
xlabel('Timestep(dt)')
ylabel('Time of Execution')
title('comparion of time taken to solve 1D Linear convective Equation')
          
       

Plot 1:

Plot 2:

Comparison and Conclusion:

1. From the above plots, it is clear that the solution for dt =0.1 is unstable, where as for other values of dt the solution is stable. This is because from the numerical formulation of the solution if we take K= `{C ∆t}/{∆x}`, the solution modifies to:

`u_i^(n+1) = u_i^n - K(u_i^n -u_(i-1)^n)` = `(1-K)u_i^n + u_(i-1)^n`

The (1-K) part on the RHS affects the stability of numerical discretization. If K > 1 the solution diverges to instability. The stability and accuracy of the solution increases as the value of K approaches to unity. For dt = 0.01 , the solution is accurate and smoother. If the solution is unstable, then it can be analyzed by the von Neumann of stable analysis.

2. The Explicit solution is preferred over the implicit, because of the cost of computation time and cost in order to attain stability. The accuracy of the implicit solution is to be compromised over the explicit solution.

3. The time taken for each value of dt also compared in bar graph, with the decrease of dt the number of time steps required to reach same end time increases, which inturn increase the number of iterations and time taken for solving the equation.