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

INTRODUCTION:

In this challenge, we will solve the following 1D linear convection equation to observe the propagation of a square wave with time.

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

Here, `C` is the velocity of propagation. The above equation is discretized using first-order forward differencing in time (time-marching) and first-order backward differencing in space. After applying the numerical schemes to the corresponding partial derivatives, we obtain the following explicit numerical solution to the above equation.

`u_i^{n+1} = u_i^n - {C Delta t}/{\Delta x}(u_i^n - u_{i-1}^n)`, where `u_{i \larr \text{node no.}}^{n \larr \text{time step}}`

The following set of values have been used to define the domain and solve the equation for the propagation of a square wave in this domain:

`L = 1\ m`

`C = 1\ m/s`

`t_{end} = 0.4\ s`

`dt = [0.1, 0.01, 0.001, 0.0001]\ s`

`nx = 80`

`u(t=0) = {(2\ m/s, quad text{if} quad 0.1\ m le x \le 0.3\ m),(1\ m/s, quad text{otherwise}):}`

MATLAB FUNCTION: (linear_convection1d.m)

The above numerical solution is programmed into MATLAB as a function. The arguments of the function are the variables defined above. The function outputs the grid, initial velocity profile and final velocity profile at `t = 0.4` seconds. The code is shown below.

function [x, u_init, u_out] = linear_convection1d(C,L,t_end,nx,dt)
    % Initializing the time and space grid
    x = linspace(0,L,nx); dx = x(2) - x(1);
    t = 0:dt:t_end;
    
    % Initial condition for U
    u_init = ones(1,nx);
    u_init(x>=0.1 & x<=0.3) = 2;
    
    % Initializing the output
    u_out = u_init;
    K = C*dt/dx;
    % Time loop
    for k=1:length(t)
        % Space loop
        u_out(2:end) = u_out(2:end) - K*(u_out(2:end) - u_out(1:end-1));
    end
end

MAIN CODE AND RESULTS:

The calculation for different values of `dt` is done in the main code by calling out the above function for different values of `dt`. The time taken for solving the function is recorded for each value of the time-step for comparison. The code and the results are shown below.

clear all; close all;
clc;

% Initializing the constants and variables
L = 1; %m
C = 1; %m/s
dt = [1e-4 1e-3 1e-2 1e-1]; %sec
t_end = 0.4; %sec
nx = 80;

% Solving the 1D Linear Convection equation and recording time taken
teval = zeros(1,4);
tic
[x1,u_init1,u_out1] = linear_convection1d(C,L,t_end,nx,dt(1)); teval(1) = toc;
tic
[x2,u_init2,u_out2] = linear_convection1d(C,L,t_end,nx,dt(2)); teval(2) = toc;
tic
[x3,u_init3,u_out3] = linear_convection1d(C,L,t_end,nx,dt(3)); teval(3) = toc;
tic
[x4,u_init4,u_out4] = linear_convection1d(C,L,t_end,nx,dt(4)); teval(4) = toc;

% Plotting the results
figure(1)
plot(x1,u_out1,x2,u_out2,x3,u_out3,'LineWidth',2)
hold on; grid on; ax = gca; ax.FontSize = 12;
plot(x4,u_out4)
xlim([0 1]); ylim([0.9 2.1]);
xlabel('X (m)','FontSize',14); ylabel('U (m/s)','FontSize',14)
lgd = legend(['dt = ',num2str(dt(1))],['dt = ',num2str(dt(2))],...
            ['dt = ',num2str(dt(3))],['dt = ',num2str(dt(4))]);
lgd.FontSize = 14;
title('Comparison of outputs for different time steps at T = 0.4 seconds');

figure(2)
c = categorical({'1e-4','1e-3','1e-2','1e-1'});
bar(c,teval); grid on; ax = gca; ax.FontSize = 12;
xlabel('dt (sec)','FontSize',14); ylabel('Time for execution (sec)','FontSize',14);
title('Comparison of time taken to solve 1D Linear Convective equation');

From the above plots, it is clear that the solution for `dt = 0.1` is an unstable one, whereas for other values of `dt` we get a stable solution. This is because from the numerical formulation of the solution if we take `K = C{Delta t}/{Delta 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 + Ku_{i-1}^n`

The `(1-K)` part on the RHS affects the stability of the numerical discretization. For `0 lt K lt 1` the solution provides a stable solution, whereas for `K gt 1` the solution diverges. The stability and accuracy of the solution increase as the value of `K` approaches unity which can be observed for `dt = 0.01 implies K = 0.79`. The time taken for each value of `dt` is also compared in a bar graph. As we decrease the value of `dt`, the number of time-steps required to reach the same end-time increases, which in turn increases the number of iterations and the time taken for solving the function.