Combining trapezoidal rule with upwind scheme

Question

I want to discretize and numerically solve the equation: \begin{equation} v(k)\dfrac{\partial f}{\partial z} + F(z)\dfrac{\partial f}{\partial k} = \alpha f(z,k) \,\,. \end{equation}

Discretization with central differences, \begin{equation} v(k_j)\left(\dfrac{f(z_{i+1},k_j) - f(z_{i-1},k_j)}{2\Delta z}\right) + F(z_i)\left(\dfrac{f(z_{i},k_{j+1}) - f(z_i,k_{j-1})}{2\Delta k}\right) = \alpha f(z_i,k_j) \end{equation} results in spurious oscillations. Using upwind scheme, \begin{equation} v(k_j)\left(\dfrac{f(z_{i},k_j) - f(z_{i-1},k_j)}{\Delta z}\right) + F(z_i)\left(\dfrac{f(z_{i},k_{j}) - f(z_i,k_{j-1})}{\Delta k}\right) = \alpha f(z_i,k_j) \qquad {\color{blue}{\bigg[ v(k_j)>0 \,,\, F(z_i) > 0 \bigg]}} \end{equation}

\begin{equation} v(k_j)\left(\dfrac{f(z_{i},k_j) - f(z_{i-1},k_j)}{\Delta z}\right) + F(z_i)\left(\dfrac{f(z_{i},k_{j+1}) - f(z_i,k_{j})}{\Delta k}\right) = \alpha f(z_i,k_j) \qquad {\color{blue}{\bigg[ v(k_j)>0 \,,\, F(z_i) < 0 \bigg]}} \end{equation}

\begin{equation} v(k_j)\left(\dfrac{f(z_{i+1},k_j) - f(z_{i},k_j)}{\Delta z}\right) + F(z_i)\left(\dfrac{f(z_{i},k_{j}) - f(z_i,k_{j-1})}{\Delta k}\right) = \alpha f(z_i,k_j) \qquad {\color{blue}{\bigg[ v(k_j)<0 \,,\, F(z_i) > 0 \bigg]}} \end{equation}

\begin{equation} v(k_j)\left(\dfrac{f(z_{i+1},k_j) - f(z_{i},k_j)}{\Delta z}\right) + F(z_i)\left(\dfrac{f(z_{i},k_{j+1}) - f(z_i,k_{j})}{\Delta k}\right) = \alpha f(z_i,k_j) \qquad {\color{blue}{\bigg[ v(k_j) < 0 \,,\, F(z_i) < 0 \bigg]}} \end{equation} removes these oscillations but produces inaccurate results.

Is there any other simple scheme to improve the accuracy while keeping on to stability of discretizaion? Can I combine the trapezoidal rule with upwind scheme? What does the resulting discretization look like?

Answer 1

I would recommend to try a second order accurate upwind scheme. There are many of them. Let me show two of them, I use them for solving nonstationary advection equation when the solution is approaching steady state. I call them like in the book of LeVeque on hyperbolic problems.

I write them as the finite difference scheme for the term $$ v \frac{\partial f(z)}{\partial z} \,, $$ and the extension to your 2D equation shall be straightforward.

The first scheme is derived from (time dependent) Fromm scheme and it takes the form $$ v \, \frac{f(z_{i+1}) + 3 f(z_i)-5 f(z_{i-1}) +f(z_{i-2})}{4 \Delta z} \hbox{ , if } v>0 , $$ and $$ v \, \frac{-f(z_{i+2}) + 5 f(z_{i+1})-3 f(z_{i}) - f(z_{i-1})}{4 \Delta z} \hbox{ , if } v<0 . $$ For time dependent problems the related Fromm scheme gives for many problems the most accurate numerical solution.

The second scheme can be derived from Beam-Warming second order upwind scheme and it takes the form $$ v \, \frac{3 f(z_{i}) - 4 f(z_{i-1}) + f(z_{i-2})}{2 \Delta z} \hbox{ , if } v>0 , $$ and $$ v \, \frac{-f(z_{i+2}) + 4 f(z_{i+1})-3 f(z_{i}) }{2 \Delta z} \hbox{ , if } v<0 . $$ The related scheme for nonstationary advection equation is less accurate than the Fromm scheme for many examples. The advantage is that it has a stencil that is one-sided, and it might produce less oscillations than the Fromm scheme.

You must be aware that any second order scheme can produce unphysical oscillations in numerical solution. Depending on your problem (e.g. the smoothness of initial condition and of given velocity field) for enough fine grids the oscillations can be negligible. If it is not the case you must use so-called limiting procedures to reduce the order of accuracy in some grid points. The resulting scheme is then nonlinear, because the coefficients are depending on the unknown numerical solution.

I repeat that I recommend to try it, it is not a definite solution, but if you need to care about the accuracy, it is a direction you should follow.

Answer 2

If at all possible don't use a forward difference derivative formula since it doesn't always provide convergence in all cases since it leads to a form of extrapolation. Instead always use a backward difference derivative. Remember: consistency + stability = convergence and forward difference derivatives do not guarantee consistency. Its a general rule of thumb I always use and a mistake thats not easy to always see in formulae. BTW, central difference second derivative formulas are always consistent.