I am solving the advection problem with high order numerical methods, using the method of lines. The boundary conditions and initial condition are selected in a way where I know that the exact solution should obey $0\le u \le 1$. $$ \frac{\partial u}{\partial t} + \nabla \cdot (u \textbf{v}) = 0.$$

One popular way to deal with oscillations in the numerical solutions is to use slope limiting. But these slope limiters appear to assume that the approximation (cell average?) is non-negative somewhere in each cell. When I use a higher order space discretization and a implicit time integrator method like backward Euler, the approximation in some places can be strictly negative on a cell no matter how much the time step is reduced. I've observed this behavior for advection-diffusion and other equations that obey some maximum principle.

What are common strategies or techniques to deal with unphysical oscillations in implicit numerical solutions? For explicit time integration there appears to be a lot more success. What about for fully implicit methods?

  • 1
    $\begingroup$ Didn't you have the same problem with explicit Euler? $\endgroup$ Mar 17 at 15:50

1 Answer 1


I will suppose that your equation represents a linear advection equation in the conservative form with the velocity ${\bf v}={\bf v}(x,t)$. It is well-known that numerical solutions of such equation using higher-order numerical methods can produce unphysical oscillations especially if the exact solution exhibits a large variation of the gradient $\nabla u$ that is not well resolved by the used computational grid.

The worst case scenario are the exact solutions having discontinuities. Here, only the first order accurate space discretization preserve maximal and minimal value of the solution, and any linear higher order scheme will produce somewhere unphysical oscillations. The bad news is that the magnitude of such oscillations does not decrease with the grid refinement, see below an illustrative picture for the advection in 1D with the speed 1.

enter image description here

Note that the time discretization has not yet entered the game. To avoid oscillations in higher-order space discretization, you must use nonlinear numerical scheme, i.e., having some numerical parameters that depend on the solution itself. To simplify it here as much as possible and to relate it to "slope limiters", I mention that many 2nd order numerical schemes are based somewhere on a (weighted) approximations by the backward and forward finite difference. The weights must depend on the numerical solution, so, formally, the scheme is nonlinear even for the linear advection equation. Now, the explicit time discretization can resolve any nonlinearity by the "linearization" using the available values of the numerical solution from the previous time or the initial condition.

The implicit time discretization can offer significantly improved stability conditions with respect to explicit methods including unconditional stability for the choice of discretization steps. The unavoidable price to pay is that to obtain the values of numerical solution one must solve systems of algebraic equations.

Consequently, the weights in the "slope limiters" also should depend on the unknown values of numerical solution giving you nonlinear algebraic equations to solve even for linear advection equation. That is quite bad news.

If you linearize the weights, e.g., by evaluating them in the old time values, you can simplify the algebraic equations, but you loose very likely the mentioned enhanced stability properties. Nevertheless, the stability properties of implicit schemes are so attractive for many applications that substantial research is done in this field.

In my opinion, good results are obtained if you, e.g., compute (linearize) the weights of the slope using predicted values obtained by a lower order scheme including the first order upwind scheme. In such a case with the "corrected" numerical scheme you have to solve only a system of linear algebraic equations if your origin PDE is linear.

There is another good news in this story. The upwind type of discretization methods for the linear advection equation can give systems of algebraic equations having a matrix with good structural properties, so some iterative solvers can be used in a very efficient way. I mention only fast sweeping method that is based on Gauss-Seidel iterations in different orderings of unknowns that can give you in many situations the solution in a finite number of iterations.

In general, to have non-oscillatory numerical solutions, one can not use strictly upwind type of space discretization and one must use also some "downwind" values of the numerical solution that can destroy the good structural property of the matrix with the fully implicit time discretization.

In our research group we have published recently some results on "implicit-explicit" 2nd order scheme that is unconditionally stable and uses "slope limiting" (more precisely, the TVD flux limiter functions) and it is fully upwinded in its implicit part. Clearly, to stop this self-promoting one has to search for this topic if interested in details. But I want to emphasize that one should definitely not restrict to "fully implicit schemes", but instead to consider "implicit schemes with explicit values in the stencil". If you are familiar with Crank-Nicolson time discretization methods (a trapezoidal rule), this is an example of such approach.

You mentioned only backward Euler implicit time discretization. If you are interested in higher-order implicit time discretizations, the things are getting more complicated, but I would say up to the second and third order accuracy theya are well manageable. You have to choose the so called strong stability preserving time integrators.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.