Error for the finite differences scheme -- Advection equation

Question

Consider the advection equation (1D in space) $$ \frac{\partial u}{\partial t} + V\, \frac{\partial u}{\partial x}=0 $$ and we solve it numerically on $[0,1]\times [0,1]\ni (t,x)$ using a forward scheme in time and a backward decentered in space. This scheme is of order $1$ in space and of order $1$ in time and is stable iff $V\frac{dt}{dx}\le 1$.

Let $N$ be the number of points in the subdivision of $[0,1]$ and $dt$ be the time step. Numerically, to observe the order, if I fix $N$ and $dt$ is varying the error between the analyitc solution and the numerical solution is represented in the log-plot figure:

The red curve is the linear regression between the points, but useless here because we don't observe a straight line so the order can be checked. More importantly when $dt$ is decreasing, the error is increasing ... What is the reason for that ?

• problem in the implementation ? maybe
• the fact is when we say of order $1$ in space and time it means we approximate the equation by an error which can be written: $$ C_1\, dt + C_2\, dx$$ where $C_1$ and $C_2$ are constant. But $C_1$ may depend on $dx$ and $C_2$ may depend on $dt$. So when $dt$ is decreasing, $C_2$ is not constant and the error can increase, which could explain the plot.

Note that when the ratio $V\frac{dt}{dx}$ is constant equal to $1$ and $dt,dx \to 0$, the error is constant and very small (no diffusion)

Answer 1

Let's define $$ \lambda = V \frac{dt}{dx}. $$ Can the error increase if you reduce $dt$ but not $dx$, i.e. if $\lambda < 1$? It seems so.

If the initial data is smooth, the solution can be expanded into a Fourier series. We can investigate what happens to each component of the series by assigning the solution $$ u = \exp{i(\tilde{\kappa} x - \tilde{\omega}t)}, $$ where $\tilde{\omega}$ and $\tilde{\kappa}$ denote a generic frequency and wavenumber respectively. This is a useful ansatz since the true solution will be a (probably infinite) sum of such Fourier modes.

Inserting $u$ into the PDE, the result is $$ (-\tilde{\omega} + V \tilde{\kappa}) u = 0 \quad \Rightarrow \quad \tilde{\omega} = V \tilde{\kappa}. $$ This expression is the dispersion relation for the advection equation and tells us the relation between the frequency and wavenumber of any Fourier mode of the solution. This relation tells us that a Fourier mode will retain its amplitude for all time ($\tilde{\omega}$ is real if $\tilde{\kappa}$ is real). It also tells us that the speed of each Fourier mode is $V$ (the wave speed is given by $\tilde{\omega} / \tilde{\kappa}$). This is enough to characterise the behaviour of the exact solution.

Now, after discretisation there are two qualities of the numerical solution that will keep it accurate: correct amplitude (no growth or decay of the numerical solution) and correct wavespeed (it won't drift out of phase), both of which are encoded in the numerical dispersion relation. The numerical scheme is $$ \frac{u_j^{n+1} - u_j^n}{dt} + V \frac{u_j^n - u_{j-1}^n}{dx} = 0. $$ We may derive the numerical dispersion relation much like we did for the PDE: Assign a solution $$ u_j^n = \exp{i(\kappa x_j - \omega t_n)}, $$ where $\omega$ and $\kappa$ are the numeric frequency and wavenumber respectively. Inserting $u_j^n$ into the discretisation gives $$ \left( \frac{e^{-i\omega dt} - 1}{dt} + V \frac{1 - e^{-i\kappa dx}}{dx} \right) u_j^n = 0, $$ which holds for any Fourier mode. Using $dt = \lambda dx / V$ we obtain the discrete dispersion relation $$ \frac{e^{-i \omega \lambda dx / V} - 1}{\lambda} + 1 - e^{-i \kappa dx} = 0. $$ Now, if $\lambda = 1$, this reduces to $$ e^{-i(\omega dx / V)} = e^{-i(\kappa dx)} \quad \Rightarrow \quad \omega = V \kappa, $$ which is exactly the same dispersion relation as we had for the PDE. In other words, the numerical solution will preserve its initial amplitude and it will travel with the wavespeed $V$. This is as close to a perfect numerical solution as we can hope for since there is not much else that can cause error for the linear advection equation.

If $\lambda \neq 1$ we will get a different numerical dispersion error. I'll leave the calculation to you, but (ignoring the imaginary part) a Maclaurin expansion seems to give something like $$ \omega = \frac{V \kappa}{\sqrt{\lambda}} + \mathcal{O}(dx). $$ The point is that choosing $\lambda < 1$, i.e. $dt < dx / V$ is detrimental to the approximation quality of the numerical dispersion relation. The smaller $\lambda$ is, the worse the approximation becomes, leading to larger errors in amplitude and phase. This may explain your observations.