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.


$\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);

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;

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

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)

% 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);

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

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';


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

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

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

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

1 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.

  • $\begingroup$ 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)$? $\endgroup$
    – bosco98
    Jan 26, 2022 at 2:30
  • $\begingroup$ 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$. $\endgroup$ Jan 26, 2022 at 3:01
  • $\begingroup$ 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$? $\endgroup$
    – bosco98
    Jan 26, 2022 at 3:13
  • $\begingroup$ 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. $\endgroup$ Jan 26, 2022 at 3:16
  • 3
    $\begingroup$ 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. $\endgroup$ Jan 26, 2022 at 6:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.