Stability in discretization of a PDE

Question

Suppose I want to numerically solve for $f(x,k)$ the one-dimensional Boltzmann equation for electrons in steady-state condition, given as: \begin{equation} \left( \dfrac{\hbar k}{m} \right) \dfrac{\partial f}{\partial x} + E \dfrac{\partial f}{\partial k} = \dfrac{f-f_0}{\tau} \end{equation} where all other symbols represent constant values.

In order to discretize the equation, I have defined the $(x,k)$ grid and I identify the coordinates of real space and momentum space with $x[i]$ and $k[j]$, respectively.

\begin{equation} \left( \dfrac{\hbar k_j}{2m} \right) \Big[f(i+1,j) - f(i-1,j)\Big] \Delta k + \dfrac{E}{2} \Big[ f(i,j+1) - f(i,j-1) \Big]\Delta x = \dfrac{f(i,j) - f_0(i,j)}{\tau} \Delta x \Delta k \end{equation} Now, the above matrix equation gives a left-hand side matrix which is not stable for arbitrary values of constants and results in very ugly oscillations.

I have seen that many people use a so-called "staggered grid" (consisting of direct/adjoint grid nodes), or use various splitting methods (even/odd or positive/negative splitting) in order to solve this problem. The more I read on these stabilization schemes, the more confused I get about the main idea. Can somebody help me with an easy-to-understand explanation or reference on the basic idea and the details of its implementation in such cases?

Answer 1

Analogous equations are considered in different applications e.g. stationary advection equation with right hand side. Looking to your equation through this view, you should try one-sided finite differences instead of central differences that are known to produce oscillations.

For instance if your input data are such that $$ \dfrac{\hbar k}{m} \ge 0 \,, \quad E \ge 0 $$ then you should use backward differences obtaining \begin{equation} \left( \dfrac{\hbar k_j}{m} \right) \Big[f(i,j) - f(i-1,j)\Big] \Delta k + E \Big[ f(i,j) - f(i,j-1) \Big]\Delta x = \dfrac{f(i,j) - f_0(i,j)}{\tau} \Delta x \Delta k \end{equation}

If those coefficients have different signs, choose forward differencing. I expect that the oscillations will disappear, otherwise the other data may produce the difficulties.

The only drawback is that the suggested "upwind" method in this answer is only first order accurate, so you may need to use fine grids to obtain satisfactory results. If this is an issue for you, you can use one-sided finite difference formulas of second order accuracy, see e.g. wikipedia.

P.S. The staggered grids are used, to my experience, when you have more than one unknown function, e.g. at least two PDEs.