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

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):

• How does the solution's maximum decay over time? I can't see that from that last graph. Jan 26, 2022 at 0:18
• @WolfgangBangerth I have edited in a bigger/better picture Jan 26, 2022 at 0:22

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.
• I think I'm misunderstanding something about the semidiscretized problem. Shouldn't we have the first component of $y'(t)$ equal to $y_0'(t)$? Which should then equal (based on what I have in my code, but this may be "bad" to do) $y_1(t) - y_0(t)$ divided by $dx$, ie the first row of $A$ is [-1 1 0 ... 0], and the first component of $y(t)$ is $y_0(t)$ instead of $y_1(t)$? Jan 26, 2022 at 2:30
• While possible to have $y_0(t)$ as the first component, it is unnecessary. It's value is already specified by the boundary condition, so there is no need to pose it as an ODE and solve for it. Note that the first component of the ODE I wrote is $y'_1(t) = (y_0(t)-y_1(t))/dx$. Jan 26, 2022 at 3:01
• Ah, I see, you are correct we should know $y_0(t)$ from the boundary condition... But I am still confused about the equation for $y_1'(t)$. That doesn't look like a backward or forward difference, did you perhaps mean $y_1'(t) = (y_1(t) - y_0(t))/dx$? Jan 26, 2022 at 3:13
• No, I think the signs are correct, and it's a normal backward difference. Plug in $j=1$ to the first equation in your original question. Jan 26, 2022 at 3:16
• And if you want a different view on this phenomenon, insert the Taylor expansion of $y_{k-1}$ to see that the discretized equation is close to $u_t=-u_x+\frac{\Delta x}{2}u_{xx}+O(\Delta x^2)$, so that there indeed is a dissipation/diffusion term. Jan 26, 2022 at 6:24