Linear Convection - Solving 1D linear convection equation for different values of time steps (dt)

Objective:-

Effect of time step size on the solution

Write a function that accepts time step as an argument and solves the same problem. You will compare the effect of time step on the numerical solution.

1. Set n = 80

2. Time steps to use = 1e-4, 1e-3, 1e-2, and 1e-1

Expected output

1. In a single graph compare the numerical solution at t = 0.4 seconds

2. Plot the simulation time as a function of the time step

Software used:- MATLAB

Solution:-

Matlab code

clear all
close all
clc
% using tic & toc command to calculate the simulation time at different
% value of dt
tic
[u0,x0] = solution_dt(1e-1) % calling the function and finding the solution at dt=1e-1
toc
simulation_time1 = toc

tic
[u1,x1] = solution_dt(1e-2) % calling the function and finding the solution at dt=1e-2
toc
simulation_time2 = toc

tic
[u2,x2] = solution_dt(1e-3) % calling the function and finding the solution at dt=1e-3
toc
simulation_time3 = toc

tic
[u3,x3] = solution_dt(1e-4) % calling the function and finding the solution at dt=1e-4
toc
simulation_time4 = toc

% plotting the graph to compare the numerical solution at t = 0.4 seconds,
% at different values of dt
figure(1)
plot(x0,u0,'b','linewidth',2)
hold on
plot(x1,u1,'r','linewidth',2)
hold on
plot(x2,u2,'bla','linewidth',3)
hold on
plot(x3,u3,'y','linewidth',2)
legend('Profile at dt=1e-1','Profile at dt=1e-2','Profile at dt=1e-3','Profile at dt=1e-4','location','southeast')
xlabel('x')
ylabel('velocity')
title('Effect of time step size on the solution')
axis([0 1 0 2.5])
grid on

% plotting the bar graph for simulation time as a function of time step (dt)
figure(2)
c = categorical({'dt=1e-1','dt=1e-2','dt=1e-3','dt=1e-4'})
simulation_time = [simulation_time1 simulation_time2 simulation_time3 simulation_time4]
bar(c,simulation_time)
legend('Simulation time','location','best')
xlabel('time steps')
ylabel('Simulation time')
title('Simulation time as a function of the time step')

Function for linear convection

function [u,x] = solution_dt(dt)

L = 1; % Domain lenght in meter
n = 80; % no. of grid points
c = 1; %linear convection velocity
x = linspace(0,L,n);% mesh
dx = x(2) - x(1);
time = 0.4; % total time in seconds
nt = time/dt; % total no. of time steps

%start point for square wave
x_start = 0.1; 
n_start = round((x_start/dx) + 1); % finding the integer value of the node for x_start

%stop point for square wave
x_stop = 0.3;
n_stop = round((x_stop/dx) + 1); % finding the integer value of the node for x_stop

u = ones(L,n);
u(:,n_start:n_stop) = 2;
uold = u;

%time loop
  for k = 1:nt
    % space loop
     for i = 2:n
       u(i) = uold(i) - (c*(dt/dx))*(uold(i) - uold(i-1)); % solving linear convection equation
     end
   %update old velocities
   uold = u;
   end

end

Report:

• First we have to create a function for liner convection, which we can use later in the main programe to calculate linear convection solution at different values of "dt".
• inside the function we define the value of Domain length (L), no. of grid points (n), linear convetion velocity (c), mesh, dx, time and no. of time steps (nt)
• After this, we calculate the node value at which the square wave starts and stops for x=0.1 to x=0.3
• After that we define the value of "u" for the square wave
• Then using for loop we start the space loop from i=2:n. This space loop is defined inside the time loop, which is from k=1:nt.
• Both the for loops together find the solution for linear convection.
• Then we start the main program by calling the function created eralier, we call the function at four given values of "dt"  which creates four different velocity profile.
• we use "tic" and "toc" command while calling the function to calculate the simulation time. Which will be later used to plot a bar graph between simulation time vs dt.
• Then we plot figure (1) all four velocity profiles in one single graph.
• after that we plot figure (2) which is the bar graph for simulation time vs dt.

Figure(1)

• For dt=1e-1, the velocity profile is not an accurate solution. There are lot of wiggles and oscillations in the profile. These errors are due to Gibbs phenomenon.
• For dt=1e-2, the peofile is much closer to a square profile. The valocity profile has its peak value at 2, which is the peak value of initial velocity profile.
• For dt=1e-3 and 1e-4, the velocity profile is almost coinciding eachother. The paek value fo the velocity profile is less than the initial velocity profile.

Conclusion:- Large value of dt gives totaly inaccurate solution. With the decreas in time step (dt), The velocity profile moves closer to the initial velocity profile in its shape and paek value but by further decreasing the time step the solution deviates from the original solution.

figure(2)

• From the above bar graph we can clearly see that simulation time is highest for   dt=1e-1.
• Simulation time is lowest for lowest value of time step dt=1e-4.
• Hence we can say that for large value of time step (dt) the simulation time is high.