\documentstyle[11pt]{article}
\begin{document}

\begin{center}
{\bf Beam and Warming Algorithm for Compressible Navier-Stokes Equations}
\end{center}
\vspace{10pt}

\noindent{\bf I.\ \ Introduction}
\vspace{10pt}

The Navier-Stokes equations in conservative form are given by

\begin{eqnarray}
\partial_t\,\rho + \partial_j \left( \rho\,u_j \right) & = & 0\ , \nonumber \\
&& \nonumber \\
\partial_t \left( \rho\,u_i \right) + \partial_j \left\lbrack \rho u_j u_i +
p \delta_{ij} - \tau_{ij} \right\rbrack & = & \rho F_i\ , \nonumber \\
&& \nonumber \\
\partial_t (\rho E) + \partial_j \left\lbrack \rho u_j H - u_i \tau_{ij} -
k \partial_j T \right\rbrack & = & \rho u_i F_i + \dot{q}\,,
\end{eqnarray}

\noindent where

\begin{equation}
\tau_{ij} = \mu \left\lbrack \partial_j u_i + \partial_i u_j \right\rbrack
+ \lambda\,\delta_{ij}\,\partial_k\,u_k\ ,
\end{equation}

\begin{description}
\item{$\rho$:} density
\item{$u_i$:} velocity component in $i$
\item{$p$:} pressure
\item{$F_i$:} external force in $i$
\item{$E$:} total internal energy $ = e + \frac{u_i\,u_i}{2}$, where $e$ is
the specific internal energy
\item{$H$:} total enthalpy $ = h + \frac{u_i\,u_i}{2}$, where $h$ is the
specific enthalpy
\item{$T$:} temperature
\item{$k$:} thermal conductivity
\item{$\mu$:} absolute viscosity
\item{$\lambda$:} second coefficient of viscosity ($= - \frac{2}{3}\mu$ upon
invoking Stokes assumption)
\end{description}
\vspace{10pt}

\noindent A perfect gas is assumed, in which case,

\begin{equation}
e = E - \frac{u_i\,u_i}{2} = C_v T \Rightarrow T = \frac{1}{C_v}
\left\lbrack E - \frac{u_i\,u_i}{2} \right\rbrack\ .
\end{equation}

\noindent Furthermore, since $p = \rho R T$, then

\begin{equation}
p = \rho R T = \rho \frac{R}{C_v} \left\lbrack E - \frac{u_i\,u_i}{2}
\right\rbrack = (\gamma -1) \rho \left\lbrack E - \frac{u_i\,u_i}{2} \right\rbrack
\ .
\end{equation}
\vspace{10pt}

To simplify the subsequent development of the Beam and Warming Algorithm,
consider a two-dimensional flow in cartesian coordinates.  Expanding the
Navier-Stokes equations gives

\[
\begin{array}{rl}
\rho_t\ + & (\rho u)_x + (\rho v)_y = 0\ , \\
& \\
(\rho u)_t\ + & \left( \rho u^2 \right)_x + (\rho uv)_y + p_x - \left\lbrack
(2\mu + \lambda ) u_x + \lambda v_y \right\rbrack_x - \left\lbrack \mu
\left( u_y + v_x \right) \right\rbrack_y = \rho F_x\ , \\
& \\
(\rho v)_t\ + & (\rho uv)_x + \left(\rho v^2 \right)_y + p_y - \left\lbrack
\mu \left( u_y + v_x \right) \right\rbrack_x - \left\lbrack (2\mu + \lambda
) v_y + \lambda u_x \right\rbrack_y = \rho F_y\ , \\
& \\
(\rho E)_t\ + & (\rho u H)_x + (\rho v H)_y - \left\lbrack 2\mu u u_x + \mu v
\left( u_y + v_x \right) + \lambda u \left( u_x + v_y \right) + kT_x
\right\rbrack_x \\
& \\
- & \left\lbrack 2\mu v v_y + \mu u \left( u_y + v_x \right) + \lambda v
\left( u_x + v_y \right) + k T_y \right\rbrack_y = \rho u F_x + \rho v F_y +
\dot{q}\ .
\end{array}
\]
\vspace{10pt}

\noindent Rearrange the above equations according to
\begin{eqnarray}
\rho_t + (\rho u)_x + (\rho v)_y & = & 0\ , \nonumber \\
&& \nonumber \\
(\rho u)_t + \left( \rho u^2 + p\right)_x + (\rho u v)_y & - & \left\lbrack
(2\mu + \lambda ) u_x \right\rbrack_x - \left\lbrack \lambda v_y
\right\rbrack_x \nonumber \\
& - & \left\lbrack \mu v_x \right\rbrack_y - \left\lbrack \mu u_y
\right\rbrack_y = \rho F_x\ , \nonumber \\
&& \nonumber \\
(\rho v)_t + (\rho u v)_x + \left( \rho v^2 + p \right)_y & - & \left\lbrack
\mu v_x \right\rbrack_x - \left\lbrack \mu u_y \right\rbrack_x \nonumber \\
& - & \left\lbrack \lambda u_x \right\rbrack_y - \left\lbrack (2\mu +
\lambda ) v_y \right\rbrack_y = \rho F_y\ , \nonumber \\
&& \nonumber \\
(\rho E)_t + (\rho u H)_x + (\rho v H)_y & - & \left\lbrack 2\mu u u_x + \mu
v v_x + \lambda u u_x + kT_x \right\rbrack_x \nonumber \\
& - & \left\lbrack \mu v u_y + \lambda u v_y \right\rbrack_x \nonumber \\
& - & \left\lbrack \mu u v_x + \lambda v u_x \right\rbrack_y \nonumber \\
& - & \left\lbrack 2\mu v v_y + \mu u u_y + \lambda v v_y + kT_y \right\rbrack_y
\nonumber \\
& = & \rho u F_x + \rho v F_y + \dot{q}\ .
\end{eqnarray}

\noindent Next, define the following vector functions:

\[
{\bf q} = \left\lbrack \begin{array}{l} \rho \\
\rho u \\
\rho v \\
\rho E \\ \end{array}
\right\rbrack
; \qquad {\bf F} = \left\lbrack \begin{array}{l} \rho u \\
\rho u^2 + p \\
\rho u v \\
\rho u H \\ \end{array}
\right\rbrack
; \qquad {\bf G} = \left\lbrack \begin{array}{l} \rho v \\
\rho u v \\
\rho v^2 + p \\
\rho v H \\ \end{array}
\right\rbrack\ ;
\]

\[
{\bf U}_I = \left\lbrack \begin{array}{l} 0 \\
(2\mu + \lambda ) u_x \\
\mu v_x \\
2\mu u u_x + \mu v v_x + \lambda u u_x + k T_x \\ \end{array}
\right\rbrack
; \qquad {\bf U}_E = \left\lbrack \begin{array}{l} 0 \\
\lambda v_y \\
\mu u_y \\
\mu v u_y + \lambda u v_y \\ \end{array}
\right\rbrack\ ;
\]

\[
{\bf V}_I = \left\lbrack \begin{array}{l} 0 \\
\mu u_y \\
(2\mu + \lambda )v_y \\
2\mu v v_y + \mu u u_y + \lambda v v_y + k T_y \\ \end{array}
\right\rbrack
; \qquad {\bf V}_E = \left\lbrack \begin{array}{l}0 \\
\mu v_x \\
\lambda u_x \\
\mu u v_x + \lambda v u_x \\ \end{array}
\right\rbrack\ ;
\]

\noindent and, finally,

\begin{equation}
{\bf H} = \left\lbrack \begin{array}{l} 0 \\
\rho F_x \\
\rho F_y \\
\rho u F_x + \rho v F_y + \dot{q} \\ \end{array}
\right\rbrack\ .
\end{equation}
\vspace{10pt}

\noindent With the above definitions, the Navier-Stokes equations may be
written as

\begin{equation}
\frac{\partial{\bf q}}{\partial t} + 
\frac{\partial{\bf F}}{\partial x} +
\frac{\partial{\bf G}}{\partial y} = 
\frac{\partial{\bf U}_I}{\partial x} +
\frac{\partial{\bf U}_E}{\partial x} +
\frac{\partial{\bf V}_I}{\partial y} +
\frac{\partial{\bf V}_E}{\partial y} + {\bf H}
\end{equation}

\noindent{\bf Note:}
$\frac{\partial{\bf U}_I}{\partial x}, 
\frac{\partial{\bf V}_I}{\partial y}$ do not contain cross-derivatives, and
$\frac{\partial{\bf U}_E}{\partial x},
\frac{\partial{\bf V}_E}{\partial y}$ do contain cross-derivatives.
\vspace{10pt}

From equation~(6), it follows that

\[
\begin{array}{ll}
{\bf F} = {\bf F}({\bf q})\ , & {\bf G} = {\bf G}({\bf q})\ , \\
{\bf U}_I = {\bf U}_I \left( {\bf q},\, {\bf q}_x \right)\ , & {\bf U}_E = {\bf U}_E 
\left( {\bf q},\, {\bf q}_y \right)\ , \\
{\bf V}_I = {\bf V}_I \left( {\bf q},\,
{\bf q}_y \right)\ , & {\bf V}_E = {\bf V}_E \left( {\bf q},\, {\bf q}_x \right)\
. \\
\end{array}
\]
\clearpage

\noindent{\bf II.\ \ Second-order Temporal Discretization}
\vspace{10pt}

Let the equations be discretized at $t = (n+1)\Delta t$.  A Taylor
series expansion about $t+\Delta t$ (or $(n+1)\Delta t$) yields

\[
\hspace{3em} {\bf q}^{n\phantom{-1}}  = {\bf q}^{n+1} - \Delta t\,
\frac{\partial {\bf q}}{\partial t} \bigg|_{n+1} +
\frac{\Delta t^2}{2}\, \frac{\partial^2 {\bf q}}{\partial t^2} \bigg|_{n+1} +
O \left( \Delta t^3 \right)\ , \hspace{5.8em} {\rm (8a)} 
\]
\[
\hspace{3em} {\bf q}^{n-1} = {\bf q}^{n+1} - 2\Delta t\,
\frac{\partial{\bf q}}{\partial t} \bigg|_{n+1} + 2\Delta t^2\,
\frac{\partial^2 {\bf q}}{\partial t^2} \bigg|_{n+1} + O (\Delta t)^3\ .  
\hspace{5.3em} {\rm (8b)}
\]

\setcounter{equation}{8}

\noindent Subtracting equation~8(b) from four times equation 8(a) yields

\[
4{\bf q}^n - {\bf q}^{n-1} = 3{\bf q}^{n+1} -
2\Delta t\, \frac{\partial{\bf q}}{\partial t} \bigg|_{n+1} + O
\left( \Delta t^3 \right)\ , 
\]

\noindent or, upon rearranging the above,

\begin{equation}
{\bf q}^{n+1} - {\bf q}^n = \frac{2}{3} \Delta t
\frac{\partial{\bf q}}{\partial t} \bigg|_{n+1} + \frac{1}{3}
\left\lbrack {\bf q}^n - {\bf q}^{n-1} \right\rbrack +
O\left( \Delta t^3 \right)\,. 
\end{equation}

\noindent Denote

\[
\Delta {\bf q}^n  =\ {\bf q}^{n+1} - {\bf q}^n\ , 
\]
\[
\frac{\partial}{\partial t} \left( \Delta {\bf q}^n \right)  =\
\frac{\partial {\bf q}}{\partial t} \bigg|_{n+1} -
\frac{\partial {\bf q}}{\partial t} \bigg|_n\ . 
\]

\noindent Upon substituting the above into equation~(9),

\begin{equation}
\Delta {\bf q}^n = \frac{2}{3} \Delta t\, \frac{\partial}{\partial
t} \Delta {\bf q}^n + \frac{2}{3} \Delta t\,
\frac{\partial {\bf q}}{\partial t} \bigg|_n + \frac{1}{3}
\Delta {\bf q}^{n-1} + O \left( \Delta t^3 \right)\ . 
\end{equation}
\vspace{10pt}

\noindent {\bf Note:}\ (1) Equation~(9) indicates that $\frac{\partial
{\bf q}}{\partial t} \Big|_{n+1}$ is second order accurate in $t$. (2)
To achieve the above second-order accuracy, the time steps $\Delta t$
must be uniform.
\vspace{10pt}

A general form of the temporal discretization which encompasses
equation~(10) is given by

\begin{eqnarray}
\Delta {\bf q}^n & = & \frac{\theta\Delta t}{1+\alpha}\,
\frac{\partial}{\partial t} \left( \Delta {\bf q}^n \right) +
\frac{\Delta t}{1+\alpha}\, \frac{\partial {\bf q}}{\partial t}
\bigg|_n + \frac{\alpha}{1+\alpha} \Delta {\bf q}^{n-1} 
\nonumber \\
& + & O \left( \Big( \theta - \frac{1}{2} - \alpha \Big) \Delta t^2,\,
\Delta t^3 \right)\ .
\end{eqnarray}

\noindent The choice $\theta =1$, $\alpha = 1/2$ recovers the
second-order accurate equation~(10).
\vspace{10pt}

We now work with the general form equation~(11).  Upon substitution of
equation~(7) into equation~(11) yields

\begin{eqnarray}
\Delta {\bf q}^n & = & \frac{\theta\Delta t}{1+\alpha}
\left\lbrack \frac{\partial}{\partial x} \Bigl\lbrace -
\Delta {\bf F}^n + \Delta {\bf U}_I^n +
\Delta {\bf U}_E^n \Bigr\rbrace \nonumber \right. \\
& & \hspace{2em}+\  \frac{\partial}{\partial y}
\Bigl\lbrace - \Delta {\bf G}^n + \Delta {\bf V}_I^n +
\Delta {\bf V}_E^n \Bigr\rbrace + \Delta {\bf H}^n \left.
\right\rbrack \nonumber \\
& + & \frac{\Delta t}{1+\alpha} \left\lbrack \frac{\partial}{\partial
x} \Bigl\lbrace - {\bf F}^n + {\bf U}_I^n +
{\bf U}_E^n \Bigr\rbrace + \frac{\partial}{\partial y}
\Bigl\lbrace - {\bf G}^n + {\bf V}_I^n + {\bf V}_E^n
\Bigr\rbrace + {\bf H}^n \right\rbrack \nonumber \\
& + & \frac{\alpha}{1+\alpha} \Delta {\bf q}^{n-1} + O \left(
\Big( \theta - \frac{1}{2} - \alpha \Big) \Delta t^2 ,\, \Delta t^3
\right), 
\end{eqnarray} 

\noindent where the following notation is adopted:

\begin{eqnarray}
{\bf F}^n & = & {\bf F} \left( {\bf q}^n \right)\ ,
\nonumber \\
\Delta {\bf F}^n & = & {\bf F} \left( {\bf q}^{n+1}
\right) - {\bf F} \left( {\bf q}^n \right)\ . 
\end{eqnarray}

\noindent Similarly, the same notation is used for the other vector functions.
\vspace{10pt}

\noindent{\bf Note:} Equation~(12) is a nonlinear equation in
$\Delta {\bf q}^n$.  Once an appropriate spatial differencing
scheme is employed, the equations could be solved by using a
Newton-Raphson method\@; however, this is an expensive iterative procedure.

Three simplifications to equation~(12) are now made to allow for an
efficient solution procedure.
\vspace{10pt}

\noindent{\bf Step 1:}

\noindent From equation~(12), it follows that $\Delta {\bf q}^n \sim O (\Delta
t)$.  Also,

\begin{equation}
\left\lbrack \Delta {\bf q}^n \right\rbrack_x = \left\lbrack
{\bf q}^{n+1} - {\bf q}^n \right\rbrack_x =
{\bf q}^{n+1}_x - {\bf q}^n_x = \left\lbrack
\Delta {\bf q}_x \right\rbrack^n \sim O (\Delta t)\ .
\end{equation}

\noindent Similarly, for the $y$-derivative.  Let

\[
\begin{array}{ll}
\Delta {\bf q}^n_x & = \left\lbrack \Delta {\bf q}^n
\right\rbrack_x = \left\lbrack \Delta {\bf q}_x \right\rbrack^n\ , \\
& \\
\Delta {\bf q}^n_y & = \left\lbrack \Delta {\bf q}^n
\right\rbrack_y = \left\lbrack \Delta {\bf q}_y \right\rbrack^n\ ,
\end{array}
\]

\noindent and expand the vector functions {\bf F} and {\bf G} about
${\bf q}^n$.  Thus,

\[
{\bf F}^{n+1} = {\bf F}^n + \Delta {\bf q}^n\, 
\frac{\partial {\bf F}}{\partial {\bf q}}\bigg|_n +
O \left(\Delta {\bf q}^n \right)^2\ , 
\]

\[
{\bf G}^{n+1} = {\bf G}^n + \Delta {\bf q}^n\,
\frac{\partial {\bf G}}{\partial {\bf q}}\bigg|_n +
O \left(\Delta {\bf q}^n \right)^2\ . 
\]

\noindent Using the notation in equation~(13) and $\Delta {\bf q}^n \sim
O (\Delta t)$,

\[
\begin{array}{lcl}
\hspace{8em} \Delta {\bf F}^n & = & {\bf A}^n
\Delta {\bf q}^n + O \left( \Delta t^2 \right)\ , \hspace{9em} {\rm (15a)} \\
&& \\
\hspace{8em} \Delta {\bf G}^n & = & {\bf B}^n
\Delta {\bf q}^n + O \left( \Delta t^2 \right)\ , \hspace{9em} {\rm (15b)} 
\end{array}
\]

\setcounter{equation}{15}

\noindent where

\begin{equation}
{\bf A}^n =
\frac{\partial {\bf F}}{\partial {\bf q}} \bigg|_n \ ;\hspace{1em}
{\bf B}^n = 
\frac{\partial {\bf G}}{\partial {\bf q}} \bigg|_n\ .  
\end{equation}

\noindent Similarly, the functions ${\bf U}_I$ and ${\bf V}_I$ are expanded
according to

\begin{eqnarray}
{\bf U}_I^{n+1} & = & {\bf U}_I^n + \Delta {\bf q}^n\,
\frac{\partial {\bf U}_I}{\partial {\bf q}}\bigg|_n  + 
\Delta {\bf q}_x^n\, 
\frac{\partial {\bf U}_I}{\partial {\bf q}_x}\bigg|_n  + 
O \left(\Delta t^2\, \right)\ , \nonumber \\
{\bf V}_I^{n+1} & = & {\bf V}_I^n + \Delta {\bf q}^n\,
\frac{\partial {\bf V}_I}{\partial {\bf q}}\bigg|_n  + 
\Delta {\bf q}_y^n\, 
\frac{\partial {\bf V}_I}{\partial {\bf q}_y}\bigg|_n  + 
O \left(\Delta t^2 \right)\ , 
\end{eqnarray}

\noindent or,

\begin{eqnarray}
\Delta {\bf U}_I^n & = & {\bf P} \Delta {\bf q}^n
+ {\bf R} \Delta {\bf q}^n_x + O \left( \Delta t^2\right)\ ,
\nonumber \\
\Delta {\bf V}_I^n & = & {\bf Q} \Delta {\bf q}^n
+ {\bf S} \Delta {\bf q}^n_y + O \left( \Delta t^2\right) \ ,
\end{eqnarray}

\noindent where

\begin{eqnarray}
{\bf P} =
\frac{\partial {\bf U}_I}{\partial {\bf q}}\bigg|_n & ; &
{\bf Q} =
\frac{\partial {\bf V}_I}{\partial {\bf q}}\bigg|_n\ , \nonumber \\
&& \nonumber \\
{\bf R} =
\frac{\partial {\bf U}_I}{\partial {\bf q}_x}\bigg|_n & ; &
{\bf S} =
\frac{\partial {\bf V}_I}{\partial {\bf q}_y}\bigg|_n .
\end{eqnarray}

\noindent Substituting equations~(15)--(19) into equation~(12) gives

\begin{eqnarray}
\Delta {\bf q}^n & = & \frac{\theta\Delta t}{1+\alpha} \left\lbrack
\frac{\partial}{\partial x}\left\lbrace - {\bf A}^n
\Delta {\bf q}^n + {\bf P} \Delta {\bf q}^n +
{\bf R} \Delta {\bf q}^n_x + \Delta {\bf U}^n_E
\right\rbrace \right. \nonumber \\
& + & \frac{\partial}{\partial y} \left\lbrace - {\bf B}^n
\Delta {\bf q}^n + {\bf Q} \Delta {\bf q}^n +
{\bf S} \Delta {\bf q}^n_y + \Delta {\bf V}^n_E
\right\rbrace + \Delta {\bf H}^n \bigg] \nonumber \\
& + & \frac{\Delta t}{1+\alpha} \left\lbrack \frac{\partial}{\partial x}
\left\lbrace - {\bf F}^n + {\bf U}_I^n + {\bf U}_E^n \right\rbrace +
\frac{\partial}{\partial y} \left\lbrace - {\bf G}^n +
{\bf V}_I^n + {\bf V}_E^n \right\rbrace + {\bf H}^n
\right\rbrack \nonumber \\
& + & \frac{\alpha}{1+\alpha}\, \Delta {\bf q}^{n-1} + 0\left(
\Big(\theta - \frac{1}{2} - \alpha \Big) \Delta t^2 \,,\, \Delta t^3
\right)\ .
\end{eqnarray}

\noindent{\bf Note:} $\Delta {\bf U}_E^n$ and $\Delta {\bf V}_E^n$
are left unchanged and are not linearized like the other functions for the
following reason:

Suppose indeed that $\Delta {\bf U}_E^n$ and $\Delta {\bf V}_E^n$ were
linearlized in the same fashion as in equation~(18).  Then, upon
transferring the first term, i.e., $\frac{\theta\Delta t}{1+\alpha}
\left\lbrack \hspace{1em} \right\rbrack$, from the right-hand side to
the left-hand side, an implicit system of equations for the unknown
$\Delta {\bf q}^n$ would be set up.  Furthermore, suppose a simple
second-order central difference scheme is employed to discretize the
system of equations.  Then, for all terms other than
$\frac{\partial}{\partial x} \left\lbrack \Delta {\bf U}_E
\right\rbrack^n$ and $\frac{\partial}{\partial y} \left\lbrack \Delta
{\bf V}_E \right\rbrack^n$, a five-point stencil would result; however,
since the two terms above involve cross-derivatives, a nine-point
stencil is created.  (In three dimensions, it would be a nineteen-point
stencil).  The solution procedure would then be quite expensive.  In
order to circumvent this difficulty, the terms $\Delta {\bf U}_E^n$ and
$\Delta {\bf V}_E^n$ are treated in the following manner:
\vspace{10pt}

\noindent{\bf Step 2:}

A Taylor series expansion about $t$ gives

\[
{\bf U}_E^{n+1} = {\bf U}_E^n + \Delta t\, \frac{\partial 
{\bf U}_E}{\partial t}\bigg|_n + O \left( \Delta t^2 \right)\ ,
\]

\[ 
{\bf U}_E^{n-1} = {\bf U}_E^n - \Delta t\, \frac{\partial
{\bf U}_E}{\partial t}\bigg|_n + O \left( \Delta t^2 \right)\ . 
\]

\noindent Hence, for a uniform time step,

\[
\Delta {\bf U}_E^n = \Delta {\bf U}_E^{n-1} + O\left( \Delta
t^2 \right)\ ,
\]

\noindent and, similarly,

\begin{equation} 
\Delta {\bf V}_E^n = \Delta {\bf V}_E^{n-1} + O\left( \Delta
t^2 \right)\ . 
\end{equation}

\noindent  In effect, the cross-derivative terms are handled ``explicitly.''
With the simplification in equation~(21), equation~(20) now takes the following
form:

\[
\Delta {\bf q}^n\  -
\frac{\theta\Delta t}{1+\alpha} \left\lbrack \frac{\partial}{\partial x}
\biggl\lbrace - {\bf A}^n \Delta {\bf q}^n +
{\bf P} \Delta {\bf q}^n + {\bf R}
\Delta {\bf q}_x^n \biggr\rbrace \right.
\] 
\[
 +\  \frac{\partial}{\partial y} \left\lbrace - {\bf B}^n 
\Delta {\bf q}^n + {\bf Q} \Delta {\bf q}^n + 
{\bf S} \Delta {\bf q}_y^n \right\rbrace \bigg]\ \ 
\]
\[ 
\hspace{8em} =\ {\bf H}^{*n} +\  O\left( \Big( \theta - \frac{1}{2} - \alpha \Big)
\Delta t^2 \,,\, \Delta t^3 \right)\ , \hspace{5.0em} {\rm (22a)}
\]

\noindent where

\[
{\bf H}^{*n} = \frac{\theta\Delta t}{1+\alpha} \left\lbrack
\frac{\partial}{\partial x} \left( \Delta {\bf U}_E^{n-1} \right) +
\frac{\partial}{\partial y} \left( \Delta {\bf V}_E^{n-1} \right) +
\Delta {\bf H}^n \right\rbrack
\]
\[
\hspace{6em}\ +\ \frac{\Delta t}{1+\alpha}\, {\bf H}^n + \frac{\alpha}{1+\alpha} 
\Delta {\bf q}^{n-1}\ . \hspace{11em} {\rm (22b)}
\]
\vspace{10pt}

\noindent{\bf Note:}  The forcing function ${\bf H}^{*n}$ of the temporally
discretized equations is not the same as the physical forcing function
${\bf H}^n$.
\vspace{10pt}

An expansion of equation~(22a) gives

\begin{eqnarray*}
\Delta {\bf q}^n & - & \frac{\theta\Delta t}{1+\alpha}
\bigg[ - {\bf A}_x^n \Delta {\bf q}^n - {\bf A}^n \Delta {\bf q}_x +
{\bf P}_x^n \Delta {\bf q}^n + {\bf P}^n \Delta {\bf q}_x \\
& + & {\bf R}_x^n \Delta {\bf q}_x^n + {\bf R}^n \Delta {\bf q}_{xx}^n 
- {\bf B}_y^n \Delta {\bf q}^n - {\bf B}^n \Delta {\bf q}^n_y \\
& + & {\bf Q}_y^n \Delta {\bf q}^n + {\bf Q}^n \Delta {\bf q}_y +
{\bf S}_y^n \Delta {\bf q}_y^n + {\bf S}^n \Delta {\bf q}_{yy}^n \bigg] = 
{\bf H}^{*n}\ .
\end{eqnarray*}

\setcounter{equation}{22}

\noindent Denote
\begin{description}
\item{$\delta_y ,\, \delta_x$ :} generic first-order derivative 
\item{$\delta_{yy} ,\, \delta_{xx}$ :} generic second-order derivative 
\end{description}

\noindent With the above notation, the discretized equation may be written
compactly according to 

\begin{eqnarray}
\bigg[ {\bf I} & + & \frac{\theta\Delta t}{1+\alpha} \Big[ \left( {\bf A}_x^n -
{\bf P}_x^n \right) + \left( {\bf A}^n - {\bf P}^n - {\bf R}_x^n \right) \delta_x -
{\bf R}^n \delta_{xx} \nonumber \\
& + & \left( {\bf B}_y^n - {\bf Q}_y^n \right) + \left( {\bf B}^n -
{\bf Q}^n - {\bf S}_y^n \right) \delta_y -
{\bf S}^n \delta_{yy} \Big] \bigg] \Delta {\bf q}^n = {\bf H}^{*n}\ .
\end{eqnarray}

\noindent If three-point differences are used, the above discretization would
lead to a block pentadiagonal system of equations.
\vspace{20pt}

\noindent{\bf Step 3:} Approximate Factorization

Denote the operators ${\bf \alpha}_x$ and
${\bf \alpha}_y$ according to

\[
{\bf \alpha}_x = \left( {\bf A}_x^n -
{\bf P}_x^n \right) + \left( {\bf A}^n -
{\bf P}^n - {\bf R}_x^n \right) \delta_x -
{\bf R}^n \delta_{xx}
\]

\[ 
{\bf \alpha}_y = \left( {\bf B}_y^n -
{\bf Q}_y^n \right) + \left( {\bf B}^n -
{\bf Q}^n - {\bf S}_y^n \right) \delta_y -
{\bf S}^n \delta_{yy} 
\]

\noindent Then, equation~(23) may be written as

\begin{equation}
\left\lbrace {\bf I} + \frac{\theta\Delta t}{1+\alpha} \left( {\bf \alpha}_x +
{\bf \alpha}_y \right) \right\rbrace \Delta {\bf q}^n = 
{\bf H}^{*n} + O \left( \Big( \theta - \frac{1}{2} - \alpha \Big)
\Delta t^2 \,,\, \Delta t^3 \right)\ . 
\end{equation}

\noindent But,

\[
\left\lbrace {\bf I} + \frac{\theta\Delta t}{1+\alpha}\,
{\bf \alpha}_x \right\rbrace \left\lbrace I +
\frac{\theta\Delta t}{1+\alpha}\, {\bf \alpha}_y
\right\rbrace = \left\lbrace {\bf I} + \frac{\theta\Delta t}{1+\alpha}
\left( {\bf \alpha}_x + {\bf \alpha}_y
\right) \right\rbrace + O \left( \Delta t^2 \right)\ .
\]

\noindent Since $\Delta {\bf q}^n \sim O (\Delta t)$, equation~(24) may be
written as

\begin{equation}
\left\lbrace {\bf I} + \frac{\theta\Delta t}{1+\alpha}\,
{\bf \alpha}_x \right\rbrace \left\lbrace I +
\frac{\theta\Delta t}{1+\alpha}\, {\bf \alpha}_y
\right\rbrace \Delta {\bf q}^n = {\bf H}^{*n} + 
O \left( \Big( \theta - \frac{1}{2} - \alpha \Big)
\Delta t^2 \,,\, \Delta t^3 \right)\ . 
\end{equation}

\noindent Equations~(24) and (25) have the same formal temporary accuracy.
\vspace{10pt}

The factorization in equation~(25) permits an efficient solution procedure
consisting of a set of two block-tridiagonal system which is solved very
efficiently by the Block Thomas algorithm.  Specifically, solve

\[
\left\lbrace {\bf I} + \frac{\theta\Delta t}{1+\alpha}\, 
{\bf \alpha}_x \right\rbrace \Delta {\bf q}^* = {\bf H}^{*n}\ ,
\]

\noindent followed by

\[
\left\lbrace {\bf I} + \frac{\theta\Delta t}{1+\alpha}\, 
{\bf \alpha}_y \right\rbrace \Delta {\bf q}^n = \Delta {\bf q}^*\ .
\]
\newpage

{\bf Comment:}
\begin{itemize}
\item Cross-derivative terms are expensive to handle implicitly.
\item If possible, it is best to handle them ``explicitly'' by moving them
over to the right-hand side as demonstrated by the Beam and Warming
Algorithm.
\item A linearized (frozen coefficient assumption) stability analysis has
shown the above algorithm to be stable for $\alpha = 1/2$.  However, stable
calculations for nonlinear equations is not guaranteed.  This is usually a
problem when nonorthogonal grids are used where the cross-derivative terms
may not be small.  Furthermore, stability requirements are far more stringent
in three dimensions. 
\item The calculations at the first time step cannot be $O \left( \Delta t^2
\right)$.
\end{itemize}


\end{document}