I`m analyzing the numerical methods for the 1D convection equation for stability, consistency, and accuracy. I want to see How does this square pulse move in domain and time? Here is my code
% Convection equation - Ct + u * Cx = 0;
% IC: C(x,0)= x, 0<=x<=0.5; C(x,0) =(1-x), 0.5<=x<=1;
% u = 0.1 m/s
% Explicit Upwind Method
% nu = udelta_t/delta_x <= 1
% nu is CN - Courant Number = Beta
clear vars
clf
clear plot
clear
clc
u = 0.1; % convection velocity, m/s eg 0.1
L = 100; % Distance, m
m1 = 10; % number of spatial segments
delta_x = L/m1 % Domain discretization
m = 10; % Graphing distance
nt = 500; % Total time, s
n = 10 % number of time steps
delta_t = nt/n
nu = udelta_t/delta_x; % Beta (for stability purposes)
sf = 11; % graphical scaling factor; note sfL, eg 11
x = zeros(m+3,1); % location along the x direction
x(1) = 0;
for i = 1:m
x(i+1) = x(i) + delta_x;
end
x
t = zeros(n+1,1); % Time at different steps
t(1) = 0;
for i = 1:n
t(i+1) = t(i) + delta_t;
end
t
C = zeros(m+3,n+1);
% Initial Condition
for i = 1:m
if (x(i) >= 0) && (x(i) <=L)
C(i,1) = 1;
elseif (x(i) > 0.5) && (x(i) <= L)
C(i,1) = x(i);
end
end
% Boundary Conditions
% Interior nodes
for j = 1:n
for i = 3:m
C(i,j+1) =nuC(i-1,j) + ((1-nu)) (C(i,j));
end
end
C(1:m+3,1);
C(1:m+3,n+1);
% C(2,:) = (1.2)(C(1,:) + C(3,:));
disp([' courant number, nu = ',num2str(abs(nu))]);
%{
plot(x,u(1:m+1,n))
grid on
title(['Upwind method, Initial Displacement profile, time = ', num2str(nt), ' s, at Courant Number, nu = ', num2str(abs(nu))])
xlabel('Distance (m)')
ylabel('Velocity (m/s)')
%}
% Plot
Tvector = C(:);
Tmax = max(Tvector);
% Concentration Profile Initial
subplot(2,2,1)
plot(x,C(1:m+3,1))
grid on
axis([0 sfL 0 Tmax])
title(['Upwind, Initial Concentration profile, time = ', num2str(t(1)), ' s, at CN, nu = ', num2str(abs(nu))])
xlabel('Distance (m)')
ylabel('Concentration (mg/L)')
% Concentration Profile Final
subplot(2,2,3)
plot(x,C(1:m+3,n+1))
grid on
axis([0 sfL 0 Tmax])
title(['Upwind, Final Concentration profile, time = ', num2str(nt), ' s, at CN, nu = ', num2str(abs(nu))])
xlabel('Distance (m)')
ylabel('Concentration (mg/L)')
%plot(t,heaviside(t+1).heaviside(2-t))*
