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

Question

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

Answer 1

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?

Answer 2

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.

Answer 3

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.

Answer 4

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.