Beam and Warming scheme

In numerical mathematics, Beam and Warming scheme or Beam–Warming implicit scheme introduced in 1978 by Richard M. Beam and R. F. Warming, is a second order accurate implicit scheme, mainly used for solving non-linear hyperbolic equations. It is not used much nowadays.

Introduction

This scheme is a spatially factored, non iterative, ADI scheme and uses implicit Euler to perform the time Integration. The algorithm is in delta-form, linearized through implementation of a Taylor-series. Hence observed as increments of the conserved variables. In this an efficient factored algorithm is obtained by evaluating the spatial cross derivatives explicitly. This allows for direct derivation of scheme and efficient solution using this computational algorithm. The efficiency is because although it is three-time-level scheme, but requires only two time levels of data storage. This results in unconditional stability. It is centered and needs the artificial dissipation operator to guarantee numerical stability.

The delta form of the equation produced has the advantageous property of stability (if existing) independent of the size of the time step.

The method

Consider the inviscid Burgers' equation in one dimension

${\frac {\partial u}{\partial t}}=-u{\frac {\partial u}{\partial x}}\quad {\text{with }}x\in R$

Burgers' equation in conservation form,

${\frac {\partial u}{\partial t}}=-{\frac {\partial E}{\partial x}}$

where :$E={\frac {u^{2}}{2}}$

Taylor series expansion

The expansion of :$u_{i}^{n+1}$

$u_{i}^{n+1}=u_{i}^{n}+{\frac {1}{2}}\left[\left.{\frac {\partial u}{\partial t}}\right|_{i}^{n}+\left.{\frac {\partial u}{\partial t}}\right|_{i}^{n+1}\right]\,\Delta t+O(\Delta t^{3})$

This is also known as the trapezoidal formula.

$\therefore {\frac {u_{i}^{n+1}-u_{i}^{n}}{\Delta t}}=-{\frac {1}{2}}\left(\left.{\frac {\partial E}{\partial x}}\right|_{i}^{n+1}+\left.{\frac {\partial E}{\partial x}}\right|_{i}^{n}+{\frac {\partial }{\partial x}}\left[A(u_{i}^{n+1}-u_{i}^{n})\right]\right)$

$\because {\frac {\partial u}{\partial t}}=-{\frac {\partial E}{\partial x}}$

Note that for this equation,

$A={\frac {\partial E}{\partial u}}=u$

Tri-diagonal system

Resulting tri-diagonal system:

${\begin{aligned}&-{\frac {\Delta t}{4\,\Delta x}}\left(A_{i-1}^{n}u_{i-1}^{n+1}\right)+u_{i}^{n+1}+{\frac {\Delta t}{4\,\Delta x}}\left(A_{i+1}^{n}u_{i+1}^{n+1}\right)\\[6pt]={}&u_{i}^{n}-{\frac {1}{2}}{\frac {\Delta t}{\Delta x}}\left(E_{i+1}^{n}-E_{i-1}^{n}\right)+{\frac {\Delta t}{4\,\Delta x}}\left(A_{i+1}^{n}u_{i+1}^{n}-A_{i-1}^{n}u_{i-1}^{n}\right)\end{aligned}}$

This resulted system of linear equations can be solved using the modified tridiagonal matrix algorithm, also known as the Thomas algorithm.

Dissipation term

Under the condition of shock wave, dissipation term is required for nonlinear hyperbolic equations such as this. This is done to keep the solution under control and maintain convergence of the solution.

$D=-\varepsilon _{e}(u_{i+2}^{n}-4u_{i+1}^{n}+6u_{i}^{n}-4u_{i-1}^{n}+u_{i-2}^{n})$

This term is added explicitly at level $n$ to the right hand side. This is always used for successful computation where high-frequency oscillations are observed and must be suppressed.

Smoothing term

If only the stable solution is required, then in the equation to the right hand side a second-order smoothing term is added on the implicit layer. The other term in the same equation can be second-order because it has no influence on the stable solution if

$\nabla ^{n}(U)=0$

The addition of smoothing term increases the number of steps required by three.

Properties

This scheme is produced by combining the trapezoidal formula, linearization, factoring, Padt spatial differencing, the homogeneous property of the flux vectors (where applicable), and hybrid spatial differencing and is most suitable for nonlinear systems in conservation-law form. ADI algorithm retains the order of accuracy and the steady-state property while reducing the bandwidth of the system of equations. Stability of the equation is

$L^{2}$-stable under CFL : $|a|\,\Delta t\leq 2\,\Delta x$

The order of Truncation error is

$O((\Delta t)^{2}+(\Delta x)^{2})$

The result is smooth with considerable overshoot (that does not grow much with time).