Solving Advection Equation PDE with 2nd order Accuracy

I am solving the advection equation ∂t/∂u = -c * ∂t/∂x and imposing a small pulse of material in the beginning of the simulation, imposed in the boundary. I'm using ODE15s as a solver and am discretizing the equation spatially with a backward fininte difference scheme. I have a code that can implement a first order or second order backward discretization scheme. For the first order implementation, I can simulate and do not get osciallations. For the 2nd order scheme, I get oscillations, regardless of how fine my mesh is. Is there a way to implement this with 2nd order accuracy and not get oscillations?

clear
% Parameters
Nx = 1000;        % Number of spatial points
Lx = 10;        % Length of domain
dx = Lx / Nx;   % Spatial step
c = 1.0;        % Advection speed
tspan = [0, 20]; % Time span for simulation

% Switch between 'first' and 'second' order backward difference
order = 'second';  % Options: 'first' or 'second'
% Spatial grid
x = linspace(0, Lx, Nx);
% Initial Condition: All zeros
u0 = zeros(Nx, 1);
% Solve ODE system using ODE15s
[t, U] = ode15s(@(t, u) advection_ode(t, u, dx, c, Nx, order), tspan, u0);
% Choose a fixed observation point in space (e.g., midpoint of domain)
obs_index = round(Nx / 2);
concentration_over_time = U(:, obs_index); % Track concentration at obs_index
% Plot Concentration vs. Time
figure;
plot(t, concentration_over_time, 'b', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('Concentration at x = Lx/2');
title(['Advection Equation: ', order, ' Order Backward Difference with Pulse Injection']);
grid on;

% -----------------------
% Function for ODE System
% -----------------------
function dudt = advection_ode(t, u, dx, c, Nx, order)
dudt = zeros(Nx, 1);
% Inject a small pulse at x = 0 for a short time (0 <= t <= 0.5)
if t >= 0 && t <= 0.5
    u(1) = 1;  % Set a pulse at the inlet
else
    u(1) = 0;  % No pulse after injection time
end
if strcmp(order, 'first')
    % First-order backward difference
    for i = 2:Nx
        dudt(i) = -c * (u(i) - u(i-1)) / dx;
    end
elseif strcmp(order, 'second')
    % Second-order backward difference
    for i = 3:Nx
        dudt(i) = -c * (3*u(i) - 4*u(i-1) + u(i-2)) / (2*dx);
    end
    % Approximate first-order difference for the second point (i=2)
    dudt(2) = -c * (u(2) - u(1)) / dx;
end
% Apply the boundary condition at x = 0 (Pulse Injection)
dudt(1) = 0;  % Fix inlet value
end

Answer 1

[Torsten]

Usual second-order schemes for the advection equation will always produce oscillations. Section 3.1.3 of the enclosed article suggests a scheme that adds "artificial viscosity" near the discontinuity to avoid oscillations and is second-order accurate elsewhere. But it costs quite a bit of effort to implement it.

I suggest you use a fine grid together with the first-order scheme - this will be almost as exact as a second-order scheme. Or - if the equation doesn't become more complicated - you can solve it analytically.

At least you should integrate in two steps (with a restart from the solution of the first step) avoiding the if-statement on time

if t >= 0 && t <= 0.5
    u(1) = 1;  % Set a pulse at the inlet
else
    u(1) = 0;  % No pulse after injection time
end

[Nick]

That makes sense, thanks!