View source for WENO methods

Prof. Chi-Wang Shu accepted the invitation on 4 July 2009 (self-imposed deadline: 4 January 2010). 

<strong>WENO methods</strong>   refers to a class of nonlinear finite volume or finite difference methods which can numerically approximate solutions of hyperbolic conservation laws and other convection dominated problems with high order accuracy in smooth regions and essentially non-oscillatory transition for solution discontinuities.

High order accurate weighted essentially non-oscillatory (WENO)
schemes have gained popularity in
numerical solutions of hyperbolic partial differential
equations (PDEs) and other convection dominated problems.  The main 
advantage of such schemes is their
capability to achieve arbitrarily high order formal accuracy in
smooth regions while maintaining stable, non-oscillatory
and sharp discontinuity transitions.  The schemes are thus
especially suitable for problems containing both strong
discontinuities and complex smooth solution features.

At the heart of the WENO schemes is actually an approximation procedure,
not directly related to PDEs, hence the WENO procedure can also be used in
many non-PDE applications, including computer vision and image
processing.

The first WENO scheme was introduced in 1994 by Liu, Osher and
Chan in their pioneering paper, in which a third order 
accurate finite volume
WENO scheme was designed.  In 1996, Jiang and Shu
provided a general framework to construct arbitrary order accurate finite
difference WENO schemes, which are more efficient for multi-dimensional
calculations.  Most of the applications use the fifth order accurate WENO
scheme designed by Jiang and C.-W. Shu 1996.  

As mentioned before, at the heart of the WENO schemes is actually an approximation procedure,
not directly related to PDEs.  We use the following simple 
problem of interpolation to describe this approximation procedure.
Assume that we
are given a uniform mesh <math>\cdots < x_1 < x_2 < x_3 < \cdots</math>
and the point values of a function
<math>u_i = u(x_i)\ .</math>
We would like to find an approximation of the function <math>u(x)</math> at a point
other than the nodes <math>x_i\ ,</math> for example at the half nodes <math>x_{i+\frac 12}\ .</math>

This can be handled by the traditional approach of interpolation.  
For example, we could find a unique polynomial of degree at most two,
denoted by <math>p_1(x)\ ,</math> which interpolates the function <math>u(x)</math> at the
mesh points in the stencil <math>S_1 = \{ x_{i-2}, \, x_{i-1}, \, x_{i} \}\ .</math>
We could then use <math>u^{(1)}_{i+\frac 12} \equiv p_1(x_{i+\frac 12})</math>
as an approximation to the value <math>u(x_{i+\frac 12})\ ,</math> which is given
explicitly as 
:<math>\label{eq1}
u^{(1)}_{i+\frac 12} = \frac{3}{8} u_{i-2} - \frac{5}{4} u_{i-1} + \frac{15}{8} u_{i}
</math>

and is third order accurate
:<math>
u^{(1)}_{i+\frac 12} - u(x_{i+\frac 12}) = O(\Delta x^3)
</math>
if the function <math>u(x)</math> is smooth in the stencil <math>S_1\ .</math>
Similarly, if we choose a different stencil <math>S_2 =\{ x_{i-1}, \, x_{i}, \, x_{i+1} \}\ ,</math>
we would obtain a different interpolation polynomial <math>p_2(x)</math>
satisfying <math>p_2(x_j)=u_j</math> for <math>j=i-1, \, i, \, i+1\ .</math>
We then obtain a different approximation <math>u^{(2)}_{i+\frac 12} \equiv p_2(x_{i+\frac 12})</math>
to <math>u(x_{i+\frac 12})\ ,</math> given explicitly as
:<math>\label{eq2}
u^{(2)}_{i+\frac 12} = -\frac{1}{8} u_{i-1} + \frac{3}{4} u_{i} + \frac{3}{8} u_{i+1} ,
</math>

which is also third order accurate
:<math>
u^{(2)}_{i+\frac 12} - u(x_{i+\frac 12}) = O(\Delta x^3)  
</math>
provided that the function <math>u(x)</math> is smooth in the stencil <math>S_2\ .</math>
Finally, a third stencil <math>S_3 =\{ x_{i}, \, x_{i+1}, \, x_{i+2} \}</math> would lead to
yet another different interpolation polynomial <math>p_3(x)\ ,</math>
satisfying <math>p_3(x_j)=u_j</math> for <math>j=i, \, i+1, \, i+2</math> and giving another
approximation <math>u^{(3)}_{i+\frac 12} \equiv p_3(x_{i+\frac 12})\ ,</math> or explicitly as
:<math>\label{eq3}
u^{(3)}_{i+\frac 12} = \frac{3}{8} u_{i} + \frac{3}{4} u_{i+1} - \frac{1}{8} u_{i+2} ,
</math>

which is of course also third order accurate
:<math>
u^{(3)}_{i+\frac 12} - u(x_{i+\frac 12}) = O(\Delta x^3)
</math>
provided that the function <math>u(x)</math> is smooth in the stencil <math>S_3\ .</math>

If the function <math>u(x)</math> is globally smooth, all three approximations <math>u^{(1)}_{i+\frac 12}\ ,</math>
<math>u^{(2)}_{i+\frac 12}</math> and <math>u^{(3)}_{i+\frac 12}</math> obtained above are third order accurate.
One could choose one of them based on other considerations, e.g. smallest
local errors or stability.

If we use the large stencil <math>S= \{x_{i-2}, \, x_{i-1}, \,  x_{i}, \, x_{i+1}, \, x_{i+2} \}\ ,</math>
which is the union of all three third order stencils <math>S_1\ ,</math> <math>S_2</math> and <math>S_3\ ,</math>
then we would be able to obtain an interpolation polynomial <math>p(x)</math> of degree at most
four, satisfying <math>p(x_j)=u_j</math> for <math>j=i-2, \, i-1, \, i, \, i+1, \, i+2\ ,</math> and
giving an approximation <math>u_{i+\frac 12} \equiv p(x_{i+\frac 12})\ ,</math> or explicitly as
:<math>\label{eq4}
u_{i+\frac 12} = \frac{3}{128} u_{i-2}  - \frac{5}{32} u_{i-1}
 + \frac{45}{64} u_{i} + \frac{15}{32} u_{i+1} - \frac{5}{128} u_{i+2} ,
</math>

which is fifth order accurate
:<math>
u_{i+\frac 12} - u(x_{i+\frac 12}) = O(\Delta x^5)
</math>
provided that the function <math>u(x)</math> is smooth in the large stencil <math>S\ .</math>

An important observation, which is at the heart of the WENO procedure,
is that the fifth order approximation <math>u_{i+\frac 12}\ ,</math> defined in
\eqref{eq4} and based on the large stencil <math>S\ ,</math> can be written as a linear
convex combination of the three third order approximations
<math>u^{(1)}_{i+\frac 12}\ ,</math>
<math>u^{(2)}_{i+\frac 12}</math> and <math>u^{(3)}_{i+\frac 12}\ ,</math> defined by
\eqref{eq1}, \eqref{eq2}, \eqref{eq3} and based on the three
small stencils <math>S_1\ ,</math> <math>S_2</math> and <math>S_3</math> respectively:
:<math>\label{eq5}
u_{i+\frac 12} = \gamma_1 u^{(1)}_{i+\frac 12} + \gamma_2 u^{(2)}_{i+\frac 12}
+ \gamma_3 u^{(3)}_{i+\frac 12}
</math>

where the constants <math>\gamma_1\ ,</math> <math>\gamma_2</math> and <math>\gamma_3\ ,</math> satisfying
<math>\gamma_1 + \gamma_2 + \gamma_3=1\ ,</math> are usually referred to as the
''linear weights'' in the WENO literature. In this case they are given as
:<math>
\gamma_1= \frac{1}{16}, \qquad \gamma_2= \frac{5}{8}, 
\qquad \gamma_3= \frac{5}{16} .
</math>

If the function <math>u(x)</math> is smooth in the big stencil <math>S\ ,</math> then all
three third order approximations
<math>u^{(1)}_{i+\frac 12}\ ,</math>
<math>u^{(2)}_{i+\frac 12}</math> and <math>u^{(3)}_{i+\frac 12}</math> can be used, as well
as the fifth order approximation <math>u_{i+\frac 12}</math> given by 
\eqref{eq4} or by \eqref{eq5}.
However, if the function <math>u(x)</math> has a discontinuity point in <math>[x_{i-2}, x_{i+2}]\ ,</math>
then not all three approximations are equally good.
The classical ENO idea by Harten, Engquist, Osher and Chakravarthy (1987) is to
choose ''one of'' the three approximations
<math>u^{(1)}_{i+\frac 12}\ ,</math>
<math>u^{(2)}_{i+\frac 12}</math> and <math>u^{(3)}_{i+\frac 12}\ ,</math> defined by
\eqref{eq1}, \eqref{eq2}, \eqref{eq3} and based on the three
stencils <math>S_1\ ,</math> <math>S_2</math> and <math>S_3</math> respectively, using the information of the local smoothness
of the given data <math>u_j</math> for <math>i-2 \leq j \leq i+2\ ,</math> measured by divided differences.
This would guarantee third order accuracy and essentially non-oscillatory
performance if <math>u(x)</math> is smooth in at least one of the three
stencils <math>S_1\ ,</math> <math>S_2</math> and <math>S_3\ .</math>

On the other hand, the WENO idea is to choose the final approximation as a convex 
combination of the three third order approximations
<math>u^{(1)}_{i+\frac 12}\ ,</math>
<math>u^{(2)}_{i+\frac 12}</math> and <math>u^{(3)}_{i+\frac 12}\ :</math>
:<math>\label{eq6}
u_{i+\frac 12} = w_1 u^{(1)}_{i+\frac 12} + w_2 u^{(2)}_{i+\frac 12}
+ w_3 u^{(3)}_{i+\frac 12}
</math>

where <math>w_j \geq 0\ ,</math> <math>w_1 + w_2 + w_3 =1\ .</math>  
We would hope that the ''nonlinear weights'' <math>w_j</math> satisfy the
following requirements

* <math>w_j  \approx \gamma_j</math> if <math>u(x)</math> is smooth in the big stencil <math>S\ .</math>
* <math>w_j \approx 0</math> if <math>u(x)</math> has a discontinuity in the stencil <math>S_j</math>
but it is smooth in at least one of the other two stencils.

The choice of the nonlinear weights <math>w_j</math> relies on the ''smoothness indicator'' 
<math>\beta_j\ ,</math> which measures the relative smoothness of the function
<math>u(x)</math> in the stencil <math>S_j\ .</math>  The larger this smoothness indicator <math>\beta_j\ ,</math>
the less smooth the function <math>u(x)</math> is in the stencil <math>S_j\ .</math>  In most of the
WENO papers, this smoothness indicator is chosen as in the paper by  Jiang and Shu (1996),
:<math>\label{eq7}
\beta_j = \sum_{l=1}^{k} \Delta x^{2l-1} \int_{x_{i-\frac 12}}^{x_{i+\frac 12}}
\left( \frac{d^l}{dx^l} p_j(x) \right)^2 \, dx
</math>

where <math>k</math> is the polynomial degree of <math>p_j(x)</math> (in our
example <math>k=2</math>).  This is
clearly just a scaled sum of the square <math>L^2</math> norms of all the derivatives
of the relevant interpolation polynomial <math>p_j(x)</math> in the relevant 
interval <math>[x_{i-\frac 12}, x_{i+\frac 12}]</math> where the interpolating point
is located.  
In our example, we can easily work out these explicit formulas as
:<math>\label{eq8}
\begin{array}{l}
\beta_1 = \frac{1}{3} \left( 4 u_{i-2}^2-19u_{i-2}u_{i-1}+25u_{i-1}^2+11u_{i-2}u_i
               -31u_{i-1}u_i+10u_i^2 \right) \\
\beta_2 = \frac{1}{3} \left( 4u_{i-1}^2-13u_{i-1}u_i+13u_i^2+5u_{i-1}u_{i+1}
               -13u_iu_{i+1}+4u_{i+1}^2 \right) \\
\beta_3 = \frac{1}{3} \left( 10u_i^2-31u_iu_{i+1}+25u_{i+1}^2+11u_iu_{i+2}
               -19u_{i+1}u_{i+2}+4u_{i+2}^2 \right) 
\end{array}
</math>

Equipped with these smoothness indicators, we can now define the nonlinear
weights as
:<math>\label{eq9}
w_j = \frac{\tilde{w}_j}{\tilde{w}_1 + \tilde{w}_2 + \tilde{w}_3},
\qquad \mbox{with} \quad
\tilde{w}_j =\frac{\gamma_j}{(\varepsilon + \beta_j)^2} .
</math>

Here <math>\varepsilon</math> is a small positive number to avoid the denominator to
become zero and is typically chosen as <math>\varepsilon = 10^{-6}</math> in actual
calculations.  

The WENO procedure has been generalized and applied to various types of
schemes (finite difference, finite volume, compact schemes,
residual distribution schemes, limiters for the discontinuous Galerkin
schemes, etc.), and have been applied to various fields including
computational fluid dynamics,
computational astronomy and
astrophysics, semiconductor device simulation,
traffic flow models, and computational biology.
For a recent review of WENO schemes, we refer to the paper by Shu (2009).
For more details of the WENO procedure, we refer to the lecture notes by Shu (1998).

==References==

* A. Harten, B. Engquist, S. Osher and S. Chakravarthy, Uniformly high order essentially non-oscillatory schemes, III, Journal of Computational Physics, 71:231-303, 1987.
* G. Jiang and C.-W. Shu, Efficient implementation of weighted ENO schemes, Journal of Computational Physics, 126:202-228, 1996.
* X.-D. Liu, S. Osher and T. Chan, Weighted essentially non-oscillatory schemes, Journal of Computational Physics, 115:200-212, 1994.
* C.-W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, B. Cockburn, C. Johnson, C.-W. Shu and E. Tadmor (Editor: A. Quarteroni), Lecture Notes in Mathematics, volume 1697, Springer, 325-432, 1998.
* C.-W. Shu, High order weighted essentially non-oscillatory schemes for convection dominated problems, SIAM Review, 51:82-126, 2009.
[[Category:Numerical Analysis]]