Why is my Runge-Kutta 4 solution to the 1-D advection equation decaying so quickly?

Question

I am trying to numerically solve the advection equation $y_t + y_x = 0$ using a the "classical" Runge-Kutta 4 explicit timestepping method, along with a left-hand finite difference approximation for the x-derivative, in MATLAB.

Specifically:

$\frac{d}{dt}y(t_i,x_j)= -\frac{y(t_i,x_j) - y(t_i,x_{j-1})}{dx} = f_j(\vec{y})$, for $j=2,3,...,N$ and

$\vec{y_{n+1}} = \vec{y_n} + \frac{1}{6}dt(k_1 + 2k_2 + 2k_3 + k4),$

$t_{n+1} = t_n + dt$ for $n=0,1,2,...$ using

$k_1 = f(\vec{y_n}),$

$k_2 = f(\vec{y_n} + \frac{1}{2}dt *k_1),$

$k_3 = f(\vec{y_n} + \frac{1}{2}dt*k_2),$

$k_4 = f(\vec{y_n} + dt * k3)$, where

$f(\vec{y_n}) = (f_1(\vec{y_n}), f_2(\vec{y_n}),...,f_N(\vec{y_n}))$.

For this problem, I am looking at the initial condition:

$y(0,x) = e^{-10(x-0.25)^2}$, which has the analytic solution $y(t,x) = e^{-10(x-t-0.25)^2}$.

And I am just using the analytic solution as the boundary condition at $x=0$. The solution 'looks like' a right-travelling wave.

The issue I am having is that my scheme seems to provide a solution which has some (seemingly) large amplitude decay, ie as time goes on, the amplitude of the wave diminishes. I am wondering if there is a bug in my code (which I cannot seem to find...) or if there is something about the initial condition/PDE/RK4 scheme, or the $dt$ and $dx$ values chosen that I am missing which makes it unsuitable for this type of problem.

Apologies if there is any information missing, I am just starting my first course in Numerical Analysis and a lot of this is new to me. I can add info in with an Edit if so.

My MATLAB code (putting everything into one script):

% Define boundaries and boundary condition functions for the problem
f = 'scenario2FunctionIC'; % y(0,x) = exp(-10*(x-0.25)^2)
g = 'scenario2FunctionBC'; % y(t,0) = exp(-10*(-c*t-0.25)^2)
a = 5; % x in [0,a]
b = 2; % time in [0,b]
N = 80; % number of grid nodes in the x-direction
M = 160; % number of grid nodes in the t-direction
c = 1;

[Y] = RK4(f, g, a, b, c, N, M);

% trueSoln will contain the values of the analytic solution, for comparison
% against the numerical solution
trueSoln = zeros(M,N);
dx = a/(N-1);
dt = b/(M-1);
for i = 1:M
    for j = 1:N
        trueSoln(i,j) = exp(-10*((j-1)*dx - c*(i-1)*dt - 0.25)^2);
    end
end

colors_an = zeros(M, N, 3); % make analytic solution graph black
for i=1:M
    for j=1:N
        colors_an(i,j,1) = 0.25;
        colors_an(i,j,2) = 0.25;
        colors_an(i,j,3) = 0.25;
    end
end

x = [0:dx:a];
t = [0:dt:b];
surf(x, t, trueSoln, colors_an, 'FaceAlpha', 0.8);
hold on
surf(x,t,Y);
colorbar

function Y = RK4(f, g, a, b, c, N, M, option)
% This subroutine solves the 1-way wave equation (1st order) given by
% y_t(t, x) + c * y_x(t,x) = 0, on 0<=x<=a, 0<=t<=b
% with IC y(0,x)=f(x) and BC y(t,0)=g(t),
% using a Runge-Kutta 4 explicit timestepping scheme.
%
% Input - f=y(0,x) as a string ’f’
%       - g=y(t,0) as a string 'g'
%       - a and b are the right endpoints of [0,a] and [0,b]
%       - c is the constant in the wave equation
%       - N and M number of grid points over [0,a] and [0,b]
% Output - Y is the solution matrix

% Initialize parameters
dx = a / (N-1);
dt = b / (M-1);
Y = zeros(M, N);

% Compute first row of Y (ie values of y(t,x) at t=0)
for j = 1:N 
    Y(1,j) = feval(f, dx*(j-1)); % Using initial condition y(0, x) = f(x)
end

% Compute first column of Y (ie values of y(t,x) at x=0)
for i = 2:M
    Y(i,1) = feval(g, dt*(i-1), c);
end

F = 'computeLeftDerivative';

y_i = zeros(1,N); % current row of y(t*,x) for each time t*
for i = 1 : M-1
    
    y_i = Y(i,:);
    k1 = feval(F, y_i, dx/(-c));
    k2 = feval(F, y_i + 0.5 * dt * k1, dx/(-c));
    k3 = feval(F, y_i + 0.5 * dt * k2, dx/(-c));
    k4 = feval(F, y_i + dt * k3, dx/(-c));
    
    nextStep = y_i + dt/6 * (k1 + 2*k2 + 2*k3 + k4);
    Y(i+1,2:N) = nextStep(2:N); % keep LHS boundary condition
    
end
end

function Y = computeLeftDerivative(V, delta)
% Given a 1xN row vector V, returns a 1xN vector of left-hand finite 
% difference approximations to the first derivative, Y (where possible).
%
% At left-hand side, uses a right-hand difference.

N = size(V,2);

centrediag = [-1 ones(1,N-1)];
lowerdiag = -ones(1,N-1);
upperdiag = [1 zeros(1,N-2)];

A = zeros(N,N);
A = diag(centrediag) + diag(lowerdiag,-1) + diag(upperdiag,1);
A = A/delta;

Y = V * A';

end

function y = scenario2FunctionIC(x)
% Custom function for Scenario 2 of the solve1WayWaveEqn script.
y = exp(-10*(x-0.25)^2);
end

function y = scenario2FunctionBC(t, c)
% Custom function for Scenario 2 of the solve1WayWaveEqn script.
y = exp(-10*(-c*t-0.25)^2);
end

And finally what the output looks like (analytic solution is coloured black and is slightly transparent for comparison):

Answer 1

The cause of the decay is the highly-diffusive, first order, finite difference method used for the spatial discretization

Let's ignore the time discretization, as that does not appear to be the problem. Consider the semidiscretized advection problem $$ \begin{bmatrix} y'_1(t) \\ y'_2(t) \\ \vdots \\ y'_N(t) \end{bmatrix} = \frac{1}{dx} \begin{bmatrix} -1 \\ 1 & -1 \\ & \ddots & \ddots \\ & & 1 & -1 \end{bmatrix} \begin{bmatrix} y_1(t) \\ y_2(t) \\ \vdots \\ y_N(t) \end{bmatrix} + \frac{1}{dx} \begin{bmatrix} y_0(t) \\ 0 \\ \vdots \\ 0 \end{bmatrix}, $$ or in matrix form $$ y'(t) = A y(t) + b(t). $$ The exact solution is $$ y(t) = e^{t A} y(0) + \int_{0}^{t} e^{(t - \tau) A} b(\tau) d\tau. $$ With the eigenvalue of $A$ real and negative, we can see where the exponential decay of the solution comes from (at least away from the boundary condition).

There are several ways to address the decaying. You could use more points in space, which brings the eigenvalues closer to $0$. Probably a better option is to consider a different spatial discretization such as a higher-order finite difference method or finite volume method.