# Advection equation with finite difference: importance of forward, backward or centered difference formula for the first derivative

I am a newbie in finite difference methods, so I apologize in advance if the question is trivial.

I am trying to solve the advenction equation, i.e. $\frac{\partial \phi(x,t)}{\partial t} + v \frac{\partial \phi(x,t)}{\partial x}=0$, where $v$ is a velocity, with the finite difference method. The system is of course also solvable analytically, and it turns out that every function of the kind $f(x,t) = f(x-vt)$ is a solution. Therefore, given an arbitrary initial condition $\phi(x,0)$, the "shape" of $\phi(x,0)$ will move towards positive x (for $v>0$) or negative x (for $v<0$).

I first formulated the problem with the backward formula for the spatial derivative,

$\phi(x_i,t_{j+1}) = \phi(x_i,t_j) - v\cdot\Delta T\cdot(\phi(x_i,t_j)-\phi(x_{i-1},t_j))/\Delta x$.

The initial condition is given by a Gaussian packet centered at the origin,

$\phi(x_i,t_1) = exp(-(x_i/2\sigma)^2)$.

With $v>0$ the gaussian packet does indeed move towards right with increasing times, and maintaining always the same shape. However, if I set $v<0$, I obtain large instabilities, with an exponentially increasing amplitude. If I reformulate the problem in terms of the forward formula for the spatial derivative

$\phi(x_i,t_{j+1}) = \phi(x_i,t_j) - v\cdot\Delta T\cdot(\phi(x_{i+1},t_j)-\phi(x_{i},t_j))/\Delta x$

I obtain the opposite behavior: with a negative $v$, the packet propagates nicely towards left, but a positive $v$ leads to a large divergence of the system.

Finally, I tried to use the centered formula for the spatial derivative, $\phi(x_i,t_{j+1}) = \phi(x_i,t_j) - v\cdot\Delta T\cdot(\phi(x_{i+1},t_j)-\phi(x_{i-1},t_j))/(2\Delta x)$

and in this case I obtain divergences for both signs of $v$.

Is there a particular reason for this kind of behaviour? What should be the "most correct" way to implement this problem?

I attach here a short Matlab script that I made, and that reproduces what I described.

            L = 1; %Length
v = 1; %Speed
i=sqrt(-1);
DeltaT = ((L/abs(v))/1000)
DeltaX = 1*abs(v)*DeltaT
x=-L:DeltaX:L;
t = 0:DeltaT:(L/abs(v));
sigma = L/50;
Tbreak = 1e-5; %Parameter for animation

%% v = 1, backward formula
v=1;
phi = zeros(length(x),length(t));phi(:,1) = exp(-((x)./(2*sigma)).^2);
for n=1:(length(t)-1)
for k=2:(length(x)-1)
phi(k,n+1)=phi(k,n) - v*DeltaT*(phi(k,n)-phi(k-1,n))/(DeltaX);
end
end
figure(1);
for n=1:(length(t))
plot(x,phi(:,n))
pause(Tbreak)
end

%% v = -1, backward formula
v=-1;
phi = zeros(length(x),length(t));phi(:,1) = exp(-((x)./(2*sigma)).^2);
for n=1:(length(t)-1)
for k=2:(length(x)-1)
phi(k,n+1)=phi(k,n) - v*DeltaT*(phi(k,n)-phi(k-1,n))/(DeltaX);
end
end
figure(1);
for n=1:(length(t))
plot(x,phi(:,n))
pause(Tbreak)
end

%% v = 1, forward formula
v=1;
phi = zeros(length(x),length(t));phi(:,1) = exp(-((x)./(2*sigma)).^2);
for n=1:(length(t)-1)
for k=2:(length(x)-1)
phi(k,n+1)=phi(k,n) - v*DeltaT*(phi(k+1,n)-phi(k,n))/(DeltaX);
end
end
figure(1);
for n=1:(length(t))
plot(x,phi(:,n))
pause(Tbreak)
end

%% v = -1, forward formula
v=-1'
phi = zeros(length(x),length(t));phi(:,1) = exp(-((x)./(2*sigma)).^2);
for n=1:(length(t)-1)
for k=2:(length(x)-1)
phi(k,n+1)=phi(k,n) - v*DeltaT*(phi(k+1,n)-phi(k,n))/(DeltaX);
end
end
figure(1);
for n=1:(length(t))
plot(x,phi(:,n))
pause(Tbreak)
end

• You could also try Runge Kutta method, that can handle centre discretization. From Von Newmann stability plot we can say it will not diverge but I don't know how it handles directional independence in derivative stencil biasing. – AGN Aug 31 '17 at 0:24
• I found your question very carefully and clearly posed! If you've found come to an understanding of what puzzled you, it would be great if you wrote it up in an answer (with that care and clarity). It would help the site and consolidate your own understanding! BTW I think the responses here have the right answers, textbook-quality, but I'm not sure they address the particular aspects that puzzled you. (Of course, if one does, you should select it as the answer!). – hyperpallium Oct 9 '17 at 1:47

Okay. Let's begin with the first situation.

Your equation is: $$\frac{\partial f}{\partial t}+v\frac{\partial f}{\partial x}=0\tag{1}$$

$\textbf{Previous comments}$

There are plenty of webs that can explain you what happens with the error when you does not choose properly the values of your mesh to approximate the derivatives in an equation (Von Neumann stability analysis), but I prefer writing the following for your better understanding based on physical grounds.

You mentioned that the solutions of that equation are of the form $f(x-vt)$, therefore as you can prove that the function $f$ along the characteristic lines ($x=x_0+vt$) is constant a of the motion: $$\frac{df}{dt} = \frac{\partial f}{\partial t}+\frac{\partial f}{\partial x}\frac{dx}{dt}=0$$ because $dx/dt=v$.

$\textbf{Case 1}:v>0$

In this case, physical information propagates through time and space according these lines: $$x=x_0+vt \tag{2}$$

Your in the $n$th timestep you have the following available information for the node $i$ from $(2)$ $$x_{i} = x_{i-1}+v\,n\Delta t$$ i.e. the information at each node $i$ depends on the previous value and moves from left o right. To see this, the function $f$ evaluated at this point is $$f(x_i-v\,n\Delta t) = f(x_{i-1})$$

For this reason, the backward derivative is stable, because it uses the values $x_i$ and $x_{i-1}$, values in which the function $f$ has information.

$\textbf{Case 2}:v<0$

On the other hand if $v$ is negative, information is transmitted in other way. Now the point $x_i$ depends on which was the point $x_{i+1}$ from $(2)$: $$x_i=x_{i+1}-v\,\Delta t$$ Analogously the the previous case we have: $$f(x_i+v\,n\Delta t) = f(x_{i+1})$$

Therefore a stable scheme for this sign of $v$ is the one which uses the values at $x_{i}$ and $x_{i+1}$, because information is transmitted from right to left.

$\textbf{Last comments}:$

For the centred scheme, you use prohibited points $x_{i-1}$ and $x_{i+1}$ for the $v<0$ and $v>0$ cases respectively, where you have no available information. That is the reason of the instabilities for this centred scheme which always uses a prohibited value whatever the sign of $v$.

$\textbf{Solution}$:Two in one

For simplicity, you can try this scheme: $$v\frac{\partial f}{\partial x}\approx \frac{1}{2}(v+|v|)\frac{f_i-f_{i-1}}{\Delta x}+ \frac{1}{2}(v-|v|)\frac{f_{i+1}-f_{i}}{\Delta x}$$ This will give you the right flux whatever the sign of $v$. Can you see it?

• Thanks for your very detailed answer. While I got an intuitive idea of the problem, I still have some doubts about it. But for now I will trust you and go on with your suggestion! :) Yes I can see while your last formula would fix the problem. But it still sounds like "cheating" to me, since we are changing the equation according to the sign of $v$, even if the original differential equation does not depend on it. If you know any good resource where I can deepen my knowledge on this, it would be more than welcome! Thanks again! – Michele Cotrufo Aug 29 '17 at 18:29
• It is not about "cheat". For example when you approximate the derivative numerically you are not treating the true value but an approximation. When you solve this equation numerically, you split it into little Riemann problems (basics of finite volume methods) and solve them analytically (thats the reason of the care with the sign of $v$). Then you build up the approximate solution solving the system of as many equations as problems that arose when you divided the original one. For reference this book is excellent: Riemann Solvers and Numerical Methods for Fluid Dynamics written by E.F. Toro – HBR Aug 29 '17 at 18:43
• @MicheleCotrufo HBR From this answer, I found a video lecture applying Von Neumann Stability Analysis to forward (and backward) finite difference advection, exactly the question asked here (though not cemtral differences): 14 minutes. (Though I didn't have the background on "Fourier Modes" and the amplification factor A being complex to follow it fully) – hyperpallium Oct 9 '17 at 6:03

You stumbled across something very fundamental that is important to understand.

For any PDE we can define for each point its "domain of dependence" This is the region of space/time that is able to have an effect on the solution at that point. For your advection problem, it is just the characteristic line through that point. If the point of interest is $x_0,t_0$ the domain of dependence is the line $(x-x_0)=v(t-t_0)$ with $t<t_0$. For more general hyperbolic problems it is usually a region bounded by more than one characteristic.

The same idea applies to the discrete model. The discrete domain of dependence for a grid point P is the set of grid points where a change of data would affect the solution at P. When you use backward differences, the discrete domain of dependence is a pyramid with its vertex at P. One side extends vertically downward and consists of points Q such that $x_Q=x_P$. The other side slopes down to the left, and consists of points Q such that $(x_Q-x_P)/\Delta x=(t_Q-t_P)/\Delta t$ (in both cases $t_Q<t_P$)

A necessary, but not sufficient, condition for stability is that The discrete domain of dependence contains the true domain of dependence. This is not a mathematical result (although often presented as such) but a logical result.

Consider refining the grid in a way that preserves the ratioof $\Delta t/\Delta x$, so the shape of the discrete DoD does not change. We do this with a discrete DoD that does not contain the true DoD

                                      / O
/  O O
/.  O O O
/   O O O O
/    O O O O O


(I hope this conveys the idea. for some reason I could not import an image, x goes to the right, t goes up, the slashes represent the characteristic line) Hoping to make the answer better by refining the grid is hopeless because it will never contain the single relevant piece of information that we need. This much worse that just the error associated with extrapolation. Remember we want to solve for arbitrary (nonsmooth) initial data, and the data inside and outside the discrete DoD may be quite independent. It follows that if $v>0$, we must use (out of the three choices) backward differencing, and that $\Delta t$ must be less than $v\Delta x$. This the Courant (Courant-Friedrichs-Levy, CFL) condition.

The case with $v<0$ follows by symmetry.

The central difference is allowable by this test, but is NOT stable. I do not think that there is a simple intuitive explanation. However, there IS a stable (but very bad) method using this stencil. It is called the Lax-Friedrichs method $$\phi_j^{n+1}=\frac{1}{2}(\phi_{j-1}^n+\phi_{j+1}^n)-\frac{v\Delta t}{2\Delta x}(\phi_{j+1}^n-\phi_{j-1}^n)$$

It can be very hard to investigate the stability of more general cases, but the necessary Courant condition is universal because it is a logical requirement.

• This explains why the wrong answer is obtained, and how to avoid it. Arguably it doesn't matter, but it would be interesting to know why the discrete DoD not containing the true DoD causes it to behave wrongly in that particular way, of instability, with exponentially increasing values. (I think, being a "necessary condition for stability" means that without it, you always get instability.) – hyperpallium Oct 9 '17 at 1:37
• @Hyperpalium. The Lax-Wendroff theorem says that we have convergence if and only if we have consistency + stability. We do have consistency, but we cannot have convergence because we use the wrong data. Therefore we do not have stability. – Philip Roe Oct 9 '17 at 22:33

This is well known. I'll give you the part I understand, and await a better answer

It's to do with the physicality of the system: where information is coming from. In your 1d scenario, when fluid is advecting to the right, what happens at a point can only be affected by new fluid coming to that point from the left. Therefore, the backward difference discretization, using $x_i$ and $x_{i-1}$, is more accurate.

Doing this is called an upwind scheme (that's wind like a breeze, not wind like a wristwatch); aka upstream, because upstream/upwind is where the fluid is coming from.

Mathematically, if one extrapolates a slope beyond the data points defining it, it is not only inaccurate, but increasingly inaccurate. Furthermore, with each timestep, any oscillating effect will be amplified in each direction. This compunding effect over timesteps is a crucial ingredient of instability. However, I don't fully understand this process.

One might expect the greater accuracy of central differences to also be more stable; but note that because it uses data points on either side of the point of interest, with one upstream and one downstream, half its information is irrelevant (the downstream data point is about fluid moving away). It is unconditionally unstable - there's no way to prevent it "blowing up". Again, I don't fully understand this.

• Thanks for your answer! It helped to give me an intuitive idea of the problem! – Michele Cotrufo Aug 29 '17 at 18:31

Central difference with forward Euler in time will be unstable. But if you put in RK4 as mentioned in another comment, it will be L2 stable. This is enough for smooth solutions as in your case with a Gaussian initial condition. But for discontinuous solutions, it will still generate oscillatory solutions. But it should not blow up since it is L2 stable. Of course the solution may be useless since it is has too much error. I have a python code here https://github.com/cpraveen/numpde/blob/master/linhyp1d/rk4cs2.py that uses RK4. You can set the initial condition in ic.py file or define the function uinit with your own initial condition.