# Square pulse test of Upwind Finite Differences

Im 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))*

• Could you add more context to your question? What does it mean for a pulse to "move in domain and time"? Also: please edit your question so the code is readable; see help here: dsp.stackexchange.com/help/formatting
– MBaz
Apr 1, 2021 at 12:53
• I start with a concentration - c(x,t) - with a square shape, this is the initial condition, it is supposed that this square will have some ripples at edges , I think that is Gibbs phenomenon. SO now Im looking for the IC function with simulates the movement of this square pulse Apr 2, 2021 at 17:11