Advection equation in 2D using finite differences - the scheme works, but the pulse loses "energy"

Question

I am trying to solve the following equation

$ \partial_t g(x,y,t) = - v\left(\partial_x + \partial_y\right) g(x,y,t)$

using finite differences (here $v>0$). The equation is also solvable analytically, and every function of the form $g(x,y,t) = g(x-vt,y-vt)$ is a solution. So, an initial "shape" $g(x,y,0)$ will simply translate along the $(1,1)$ direction in the x-y plane for increasing times.

I am using the following scheme, with backward differences for space and forward differences for time,

$g(i,j,n+1) = g(i,j,n) - v\;dt\frac{ \displaystyle g(i,j,n) - g(i-1,j,n) + g(i,j,n) - g(i,j-1,n)}{\displaystyle dx} $

with $dx = 2\;v\;dt$. I typically use a 2D gaussian pulse as initial conditions.

The scheme works fairly well, but the amplitude of the pulse decreases as the pulse moves. Typically I lose about 10% of amplitude when the pulse moves over a distance of about $L=10^4 v\;dt$. This improves a bit (i.e., I have less decrease of the pulse amplitude) when I increase the time and spatial resolutions, but the increase in accuracy seems quite "slow".

Do you know if I can optimize my scheme in any way (other than increasing resolution)? Would it help using a central difference scheme for the time?

Answer 1

To answer your question, I will for the sake of simplicity let $v=1$ and consider the one-dimensional case \begin{equation} g_t + g_x = 0 \tag{1} \end{equation} where subscripts denote partial derivatives. The extension to 2D and $v \neq 1$ is trivial but algebraically a little more tedious.

As you say, the solution to (1) is a 'shape' propagating through the domain, i.e. it is a wave. We can decompose this shape into a series of plane waves (a Fourier series). To see what (1) does to each component of the shape, we simply assign a plane wave solution to our problem, \begin{equation} g(x,t) = \exp{i(\kappa x - \omega t)}, \tag{2} \end{equation} for some wavenumber $\kappa$ and frequency $\omega$. We see indeed that (2) is a solution to (1) if $$ \omega = \kappa. \tag{3} $$ Relation (3) is called the analytic dispersion relation and tells us the frequency of a wave with given (real) wavenumber.

Now let's discretise in space with backward differences as you have done: $$ (g_i)_t + \frac{g_i - g_{i-1}}{dx} = 0. \tag{4} $$ Plugging in a plane wave solution $$ g_i = \exp{i(\kappa x_i - \bar{\omega} t)} \tag{5} $$ now yields a numeric dispersion relation $$ \bar{\omega} = -\frac{i}{dx}(1 - \exp{(-i \kappa dx)}), \tag{6} $$ providing an expression for the numeric frequency $\bar{\omega}$ of a wave with a given (real) wavenumber $\kappa$. (6) is an approximation of (3), and indeed, as $dx \rightarrow 0$ they become identical. However, note that the numerical frequency has an imaginary and non-positive component for any choice of $\kappa$. Inserting (6) into (5) thus gives, not just a propagating wave, but an exponentially decaying one. This is what you have observed in your simulation. The decay is suppressed by choosing a small $dx$ but it is nonetheless exponential.

There are several ways to remedy this situation. A higher order backward stencil will give a numerical dispersion relation that approximates (3) better, and will have a smaller decay rate. However, if you want to rid yourself of the decay all together, choose a central stencil in space, e.g. $$ g_t + \frac{g_{i+1} - g_{i-1}}{2dx} = 0. \tag{7} $$ Inserting the plane wave (5) into (7) gives the numeric dispersion relation $$ \bar{\omega} = \frac{1}{dx} \sin{(\kappa dx)}, $$ i.e. the numeric frequency is real and no decay will be seen. Observe however that the central stencil is unstable with the forward difference in time, so use it with e.g. the 4th order Runge-Kutta method instead.