Finite difference method works for $\frac{\partial u}{\partial t} = \frac{du}{dz}$ but not for $\frac{\partial u}{\partial t} = - \frac{du}{dz}$?

Question

I am using the method of lines with forward differences to solve the transport equation $$\frac{\partial u}{\partial t} = \frac{du}{dz}$$

with initial condition $u(z, 0) = z$

and boundary condition $u(z = 1, t) = 1$

My code below works fine for this situation and produces the expected result.

If however I change the equation to

$$\frac{\partial u}{\partial t} = - \frac{du}{dz}$$

(by changing the relevant code to dt(i) = - (u(i+1) - u(i))/dz;

I get errors as can be seen here

So what is going wrong and how can I fix it? Would I be correct in thinking that in the second equation the boundary condition $(z = 1, t) = 1$ is inconsistent with the way the solution "travels" from the initial condition in this case? If not what is the cause of the error?

function mol

global n

clear
clc

% Time and spatial steps
n = 21;

% Time start and finish
t0 = 0.0;
tf = 1.0;
t_grid_points = linspace(t0,tf,n);

% Spatial start and finish
z0 = 0;
zf = 1;
z_grid_points = linspace(z0,zf,n);

dz = (zf-z0)/(n-1);

% Initial condition
for i=1:n
    u0(i) = (i-1)/(n-1);
end

reltol=1.0e-04; abstol=1.0e-04;
options=odeset('RelTol',reltol,'AbsTol',abstol);
[t,solved]=ode15s(@solve, t_grid_points, u0, options);

% --------------------%
% Plot
% --------------------%
subplot(1,1,1)
mesh(t_grid_points, z_grid_points, solved');
title('u(zf,t) vs t'); xlabel('t'); ylabel('u(zf,t)')
rotate3d on;

function dt = solve(t, u)

    for i=1:n
        if i == n
         dt(i) = 0;
        else
         dt(i) = (u(i+1) - u(i))/dz;
        end
    end

    dt = dt';
end

end

Answer 1

The equation $u_t + v u_x = 0$ describes the evolution of a density of particles each of which is moving with a velocity of $v$. As a result, for positive $v$, information is transported from left to right, and for negative $v$, information is transported from right to left. So if you want to impose a boundary condition on the right, then information must be coming in from there, and thus $v$ must be negative. Putting a boundary condition on the right with positive $v$ is an ill-posed problem at the PDE level (numerical issues aside).

A related issue, which I don't think is causing your current problems but will still cause some problems, is that your spatial differencing must be modified to deal with the sign of the velocity. Here is why.

Put yourself at position $x$ and time $t$. Move forward to time $t+h$. Back at time $t$, where were the particles which are now at $x$? Well, they moved by $vh$, so they were at $x-vh$. That means that your numerical method must "look" in the opposite direction of $v$. If it looks the other way, then you will only have hope of convergence if $u_{i-1}$ and $u_{i+1}$ stay close together at all times. Concretely, consider the forward Euler-type method:

$$u^{n+1}_i=u^n_i - hv\frac{u^n_i-u^n_{i-1}}{k} = \left ( 1 - \frac{hv}{k} \right ) u^n_i + \frac{hv}{k} u^n_{i-1}$$

This is looking for material which was at $(i-1)k$ to now be at $ik$. This is all well and good if the particles are moving to the right, but if they are moving to the left, then you're looking in the wrong direction. You should instead use the opposite spatial difference to estimate the derivative, i.e. you should use

$$u^{n+1}_i=\left ( 1+\frac{hv}{k} \right ) u^n_i - \frac{hv}{k} u^n_{i+1}$$

(where in this case $v<0$).

Note that if you take $\frac{h}{k}=\frac{1}{|v|}$ then this is actually exact (basically because your grid is properly aligned with the characteristics).