%------------------------------------------------------------------------------
% Beginning of journal.tex
%------------------------------------------------------------------------------
%
% AMS-LaTeX 1.2 sample file for journals, based on amsart.cls.
%
% Replace amsart by the documentclass for the target journal, e.g. tran-l.
%
\documentclass{mcom-l}

\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}

\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{xca}[theorem]{Exercise}

\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}

\numberwithin{equation}{section}

%    Absolute value notation
\newcommand{\abs}[1]{\lvert#1\rvert}

%    Blank box placeholder for figures (to avoid requiring any
%    particular graphics capabilities for printing this document).
\newcommand{\blankbox}[2]{%
  \parbox{\columnwidth}{\centering
%    Set fboxsep to 0 so that the actual size of the box will match the
%    given measurements more closely.
    \setlength{\fboxsep}{0pt}%
    \fbox{\raisebox{0pt}[#2]{\hspace{#1}}}%
  }%
}

%%%%%%%%%%%%%%   MY  DEFINITIONS   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\pe}{\psi}
\def\d{\delta}
\def\loc{{\rm loc}}

\def\bfR{\mbox{\boldmath$R$}}
\def\bfZ{\mbox{\boldmath$Z$}}

\begin{document}

\title{Instability of numerical schemes for a convection-reaction equation}

%    Information for first author
\author{Yong Jung Kim}
%    Address of record for the research reported here
\address{School of Mathematics, University of Minnesota, Minneapolis
, Minnesota 55455-0436}
%    Current address
%\curraddr{School of Mathematics, University of Minnesota, Minneapolis, Minnesota 55455-0436}
\email{yjkim@math.umn.edu}
%    \thanks will become a 1st page footnote.
\thanks{The first author was supported in part by NSF Grant \#000000.}

%    General info
\subjclass[2000]{65M12, 65M06}
\keywords{instability of numerical schemes, convection, reaction, splitting method}


%\dedicatory{This paper is dedicated to our authors.}

\date{\bf Version : \today}
\begin{abstract}
In this report we observe several undesired behaviors of numerical
schemes for the Burgers equation with a linear source term.
Numerical solutions based on a splitting
method are compared to exact solutions of the problem.
method. From several examples we observe how some unexpected properties of
numerical schemes grow and how they can be neutralized.
\end{abstract}

\maketitle


%\pagestyle{myheadings} \thispagestyle{plain}
\markboth{Yong-Jung Kim}{Breaks down of numerical schemes}


%\vfil\eject%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Introduction}\hskip .001in
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





%\vfil\eject%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Exact Solutions}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We consider the inviscid Burgers equation with a linear source term,
\begin{equation}
u_t+({1\over2}u^2)_x=u,\quad u(x,0)=u_0(x),
                                                       \label{cr}\end{equation}
where $u,x\in\bfR,t>0$ and $u_0(\cdot)\in L^1(\bfR)$.
Under the change of variables
\begin{equation}
w=e^{-t/2}u,\quad\xi=e^{-t/2}x,
                                                    \label{var1}\end{equation}
the convection-reaction equation (\ref{cr}) is transformed to
\begin{equation}
w_t+{1\over 2}(w^2-\xi w)_\xi=0,\quad w(\xi,0)=u_0(\xi).
                                                \label{sburgers}\end{equation}
An advantage of considering (\ref{sburgers}) is that it has non-zero
steady states,
\begin{equation}
N_{p,q}(\xi)=
\left\{
\begin{array}{rcl}
\xi&,&\quad -\sqrt{2p}<\xi<\sqrt{2q},\\
0&,&\quad {\rm otherwise},\\
\end{array}
\right.
                                                   \label{nwave}\end{equation}
which represent the asymptotic behavior of the solution.
If an N-wave is given as the initial datum, $u_0(x)=N_{p,q}(x)$,
then the solution
of (\ref{sburgers}) is simply $w(\xi,t)=N_{p,q}(\xi)$.
So the solution of (\ref{cr}) with this initial data is obtained
using the change of variables (\ref{var1}), which is
\begin{equation}
u(x,t)=
\left\{
\begin{array}{ccc}
x&,& -\sqrt{2pe^t}<x<\sqrt{2qe^t},\\
0&,& {\rm  otherwise}.\\
\end{array}
\right.
                                                   \label{Nwave}\end{equation}

An interesting feature of the convection-reaction equation (\ref{cr}) is
the existence of unstable steady states.
For an example a roll wave,
\begin{equation}
{\bfR  }_1(x)=
\left\{
\begin{array}{ccc}
x+1&,& -1<x<0,\\
x-1&,& 0<x<1,\\
0&,& {\rm otherwise},\\
\end{array}
\right.
                                                 \label{rollwave}\end{equation}
is one of them. So $u(x,t)={\bfR  }_1(x)$ is a solution of (\ref{cr}).
But, since roll waves are unstable, there are several issues related
to computing those solutions. They are discussed in \cite{JinKim}.

It is also useful to consider variables,
\begin{equation}
v=e^{-t}u,\quad s=e^t-1.
                                                    \label{var2}\end{equation}
Then (\ref{cr}) is transformed to the inviscid Burgers equation
\begin{equation}
v_s+({1\over2}v^2)_x=0,\quad v(x,0)=u_0(x).
                                                  \label{burgers}\end{equation}
The linear source term is now disappeared and a pure convection equation
appears, which has been well studied analytically and numerically.
So people might think that the property of the solution of (\ref{cr}) can
be shown easily. But there are significant difference between these two.
The existence of unstable steady states is one of them. Numerically these
two have a huge difference in time scale (\ref{var2}) and, hence, numerical
solutions of these two represent different level of properties
of the same scheme.

%\vfil\eject%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Splitting Method}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We consider notations first. In a numerical scheme the mesh
width $h:=\Delta x$ is considered to be fixed.  We know that the wave
speed of solutions of (2.1) increases exponentially in time, and hence
the time step $k:=\Delta t$ should depend on the solution.
The time step $k=k_n$ and time mesh $t^n=\sum_{i=0}^{n-1} k_i$ are defined
inductively. Suppose that $k_i,i=0,...,n-1$ are given. Then a numerical scheme
provides numerical approximations $U(\cdot,t^n)$ and 
we may define numerical wave speed $\lambda(t^n)$,
\begin{equation}
\lambda(t^n)=2\sup(|U(\cdot,t^n)|).
                                                      \label{cfl}\end{equation}
Then the time step $k_n$ is given by $k_n=h/\lambda(t^n)$ and the next
time level is $t^{n+1}=t^n+k_n$.  It is standard to use the notation
$x_j=jh,j\in\bfZ$ together with $x_{j+1/2}=x_j+h/2$.
Sometimes we use the notation $k$ for the time step instead of $k_n$
for the convenience of the notation. Since the wave speed of the Burgers
equation is $|u|$, numerical speed (\ref{cfl}) implies that
the CFL number is fixed at $0.5$.
We view $U^n_j$ as an approximation to the cell average
of the true solution, given by
$$ U^n_j\sim {\bar u}^n_j:={1\over h}\int^{x_{j+1/2}}_{x_{j-1/2}} u(x,t^n)dx.
                                                                             $$
So the discretization of the initial data is always given by taking
cell averages, $U^0_j:={1\over h}\int^{x_{j+1/2}}_{x_{j-1/2}} u_0(x)dx$.


One way to handle the convection-reaction equation (\ref{cr})
is to use a fractional step splitting method, in which one alternates
between solving the conservation law
\begin{equation}
u_t+({1\over2}u^2)_x=0
                                                        \label{c}\end{equation}
and the ordinary differential equation
\begin{equation}
u_t=u,
                                                        \label{r}\end{equation}
in each step. For the convection part we consider schemes in a conservative form
\begin{equation}
\bar U^{n+1}_j=U^n_j-{k\over h}(F(U^n_{j-p},...,U^n_{j+q})
                                 -F(U^n_{j-p-1},...,U^n_{j+q-1})), 
                                                    \label{step1}\end{equation}
where numerical flux $F$ represents the real flux $f(u)=u^2/2$.
To simplify the writing we use a notation,
$$\Delta F(U^n_j)=F(U^n_{j-p},...,U^n_{j+q}) -F(U^n_{j-p-1},...,U^n_{j+q-1}). 
                                                                              $$
The ordinary differential equation (\ref{r}) is easily solved,
$u(t^{n+1})=e^ku(t^n)$, and hence
for each time step we get
\begin{equation}
U^{n+1}_j=e^k\bar U^{n+1}_j.
                                                    \label{step2}\end{equation}

So the corresponding splitting method for (\ref{cr}) is obtained by
repeating (\ref{step1}) and (\ref{step2}). In our case it can be written
in a single equation,
\begin{equation}
U^{n+1}_j=e^k\Big(U^n_j-{k\over h}\Delta F(U^n_j)\Big).
                                                \label{splitting}\end{equation}
In the next section we consider several numerical fluxes
$F(U^n_j,U^n_{j+1})$ and consider how the splitting method (\ref{splitting})
for the given numerical flux behaves.

In the rest of this section we consider a fact that shows the close
relationship between the splitting method (\ref{splitting}) and the change
of variables (\ref{var2}). The total variation of the numerical approximation
$U^n_j$ at the time level $t=t^n$ is defined by
$$TV(U^n)=\sum_{j=-\infty}^{\infty}|U^n_j-U^n_{j-1}|.                     $$
A numerical method is called {\it total variation diminishing}
(abbreviated TVD) if
$$TV(U^{n+1})\le TV(U^n),\qquad n\in\bfZ^+.$$

It is natural to ask if the properties of the numerical scheme (\ref{step1})
is transfered to the scheme (\ref{splitting}).  For example,
we may ask if (\ref{splitting}) is TVD when (\ref{step1}) is TVD.
Since the solutions of the convection-reaction equation (\ref{cr}) is not
TVD, numerical schemes for it can not be TVD. On the other hand
we have seen that (\ref{cr}) is transformed to the Burgers equation
under (\ref{var2}). So the right question is to ask
if the transformed sequence,
\begin{equation}
s_n=e^{t^n}-1,W^n_j=e^{-t}U^n_j,
                                                     \label{var3}\end{equation}
is TVD. It is true under a reasonable assumption on the numerical flux.

\begin{lemma} Let $U^n_j$ be the solution of numerical scheme (\ref{splitting}).
If the conservative scheme (\ref{step1}) is TVD and numerical flux satisfies
\begin{equation}
F(aU^n_{j-p},...,aU^n_{j+q})=a^2F(U^n_{j-p},...,U^n_{j+q}),
                                               \label{conditionF}\end{equation}
then $W^n_j$ generated by (\ref{var3}) is also TVD.
\end{lemma}
\begin{proof}
Consider
$$\Delta s=s_{n+1}-s_n=e^{t^n}(e^k-1),\quad \Delta x=h.
                                                                            $$
From (\ref{splitting}) we get a relation for variables in (\ref{var3}),
$$W^{n+1}_j=W^n_j-{k\over h}e^{-t^n}\Delta F(e^{t^n}W^n_j).
                                                                             $$
The condition (\ref{conditionF}) for the numerical flux $F$ implies
$$W^{n+1}_j=W^n_j-{\Delta s\over\Delta x}\Delta F(W^n_j) {k\over e^k-1}.
                                                                            $$
So we get
$$V^{n+1}_j=V^n_j-{\Delta s\over\Delta x}\Delta F(V^n_j),
                                                                            $$
where $V^n_j={k\over e^k-1}W^n_j$.
It shows that $V^n_j$ is the numerical solution of the convection
equation (\ref{c}) generated by the conservative numerical scheme
(\ref{step1}), which is TVD. Since
$$TV(W^{n+1}_j)={e^k-1\over k}TV(V^{n+1}_j)
              \le {e^k-1\over k}TV(V^n_j)=TV(W^n_j),              $$
we may conclude that $W^n_j$ is also TVD.
\end{proof}


Since ${k\over e^k-1}\to 1$ fast as $k\to0$, $W_j^n$ is a good approximation
for the numerical scheme (\ref{step1}).  It shows that splitting method is
closely related to the change of variables (\ref{var2}) and it keeps the
property of the conservative scheme (\ref{step1}) inside.
The main reason we present the lemma here is to show how well
the splitting method represents the change of variable (\ref{var2}).
Even though we see in the next section that numerical solutions
based on splitting method may blowup.

%\vfil\eject%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Blow-up of Numerical Solutions}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In this section we observe several unexpected behavior of numerical schemes
based on the splitting method (\ref{splitting}) for a convection-reaction
problem,
\begin{equation}
u_t+({1\over2}u^2)_x=u,\quad u(x,0)=N_{{1\over2},{1\over2}}(x).
                                                  \label{problem}\end{equation}
As we have observed in Section 3, the splitting method simply magnify
the structure of the conservative scheme (\ref{step1}) to the scale of
(\ref{var2}). So we can say that the origin of the unexpected behavior
we observe is the numerical flux.  We use the notation $f(u)$ to show
the dependence on the flux explicitly, even though we always consider
Burgers' flux $f(u)=u^2/2$.

\subsection{ Lax-Friedrichs Scheme}
First consider the numerical flux of Lax-Friedrichs,
\begin{equation}
F(U_j,U_{j+1})={h\over 2k}(U_j-U_{j+1})+{1\over2}(f(U_j)+f(U_{j+1})).
                                                             \end{equation}
Then the corresponding splitting method (\ref{splitting}) can be written as
\begin{equation}
U^{n+1}_j=e^k\Big( {1\over2}(U^n_{j+1}+U^n_{j-1})
                      -{k\over 2h}(f(U^n_{j+1})-f(U^n_{j-1}))\Big ).
                                                      \label{lf}\end{equation}


In Figure \ref{figlf1} a numerical solution of Lax-Friedrichs scheme for
(\ref{problem}) is displayed with the analytic solution.  We can see that
the numerical solution is behind the exact one. It is because of the
numerical viscosity in parts. Usually the viscosity
causes the size reduction of the humps for N-wave like solution
(see \cite{KimTzavaras}).
But there is another unexpected problem. The numerical solution
oscillates and eventually it is separated into two N-waves,
Figure \ref{figlf1} (b).


\begin{figure}[ht]
\hskip .8 in
\special{psfile="split1a.ps" hscale=35 vscale=15 hoffset=-220 voffset=-120 }
\special{psfile="split1b.ps" hscale=35 vscale=15 hoffset= -50  voffset=-120 }
\vskip 1.7 in
\centerline{\it \quad(a) Initial data discretization \hskip 0.6 in
(b) Lax-Friedrichs (\ref{lf}) at $t=10$}
\caption{Evolution of Lax-Friedrichs :
The mesh size is $\Delta x=1/5.5$ and the initial data of
(\ref{problem}) is discretized by taking the cell average, (a).
At time $t=10$ we can observe severe oscillation, or separation of two
different waves.  The exact solution of is displayed in a line.}
\label{figlf1}
\end{figure}

The reason for the separation is clear. In the scheme (\ref{lf}) only even
numbered grid points are used to decide an odd numbered one.
So odd numbered meshes and even numbered meshes represent different solutions
of (\ref{problem}) with initial data consist of even numbered or odd numbered
meshes. So if we modify the initial discretization and let even numbered meshes
and odd numbered meshes represent the same initial data, then the separation
disappears, see Figure \ref{figlf2}.

\begin{figure}[ht]
\hskip .8 in
\special{psfile="split2a.ps" hscale=35 vscale=15 hoffset=-220 voffset=-120 }
\special{psfile="split2b.ps" hscale=35 vscale=15 hoffset=-50 voffset=-120 }
\vskip 1.7 in
\centerline{\it \quad(a) Initial data discretization \hskip 0.5 in
(b) Lax-Friedrichs (\ref{lf}) at $t=10$}
\caption{Evolution of Lax-Friedrichs :
The mesh size is $\Delta x=1/11$. The initial data is not exactly the cell
average, (a).  The solution represents a single wave at $t=10$. Even though
it contains two identical waves, i.e., $U^n_{2j}=U^n_{2j+1}$.
The exact solution of is displayed in a line.}
\label{figlf2}
\end{figure}

There are two ways to cure the symptom. The first one is to modify the flux
to avoid such a behavior.  In the example, Figure \ref{figlf2},
we actually solves two identical solutions. So it is natural to consider
a way to compute with even numbered mesh points or odd numbered ones only.
We can easily check that following modified version of Lax-Friedrichs,
\begin{equation}
\begin{array}{c}
U^{n+1/2}_j=e^{k/2}\Big( {1\over2}(U^n_{j+1}+U^n_j)
                      -{k\over 2h}(f(U^n_{j+1})-f(U^n_j))\Big ),\\
U^{n+1}_j=e^{k/2}\Big( {1\over2}(U^{n+1/2}_j+U^{n+1/2}_{j-1})
                      -{k\over 2h}(f(U^{n+1/2}_j)-f(U^{n+1/2}_{j-1}))\Big ),\\
\end{array}
                                                    \label{newlf}\end{equation}
is the scheme if the mesh size and the time steps are doubled.
So the solutions in Figure \ref{figlf3} (a) and Figure \ref{figlf2} (b)
represent the same approximation.

\begin{figure}[ht]
\hskip .8 in
\special{psfile="split3a.ps" hscale=35 vscale=15 hoffset=-220 voffset=-120 }
\special{psfile="split3b.ps" hscale=35 vscale=15 hoffset=-50 voffset=-120 }
\vskip 1.7 in
\centerline{\it (a) Modified Lax-Friedrichs (\ref{newlf})\hskip 0.8 in
(b) Non-splitting method (\ref{lfdirect})}
\caption{Evolution of Lax-Friedrichs :
Solutions at $t=10$ with the mesh size $\Delta x=1/5.5$.
The initial discretization
is same as Figure \ref{figlf1} (a).
(a) represent Figure \ref{figlf2} (b).  The wave in (b) represents
the averaged one of Figure \ref{figlf1} (b).  Solutions are stable.
The exact solution is in a line.}
\label{figlf3}
\end{figure}

The second way is to handle the source in a different way.
The splitting method simply magnifies the properties of the numerical
schemes for the convection part. So we may handle the source term in a way
that the bad behavior of the Lax-Friedrichs can be neutralized
instead of being simply magnified. We may consider a direct method,
\begin{equation}
U^{n+1}_j={1\over2}(U^n_{j+1}+U^n_{j-1})-{k\over 2h}(f(U^n_{j+1})-f(U^n_{j-1}))
          +(e^k-1)U^n_j .
                                                  \label{lfdirect}\end{equation}
In this scheme the information of even numbered meshes and odd numbered ones
are mixed due to the source term and, hence, the separation between mesh points
does not appear, Figure \ref{figlf3} (b).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{McCormack's Scheme}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We can observe a blow-up of the splitting scheme (\ref{splitting}) with
the numerical flux of McCormack,
\begin{equation}
F(U_j,U_{j+1})={1\over2}\Big[f(U_{j+1})
                         +f\Big(U_j-{k\over h}(f(U_{j+1})-f(U_j))\Big)\Big],
                                                     \label{mc}\end{equation}
Figure \ref{figmc1}. The scheme gives the correct growth and wave speed
for small a amount of time, (a).  After certain amount of time
an inadmissible shock appears at the middle of the wave around, $t=5$,
and it grows breaking the solution down.

\begin{figure}[ht]
\hskip .8 in
\special{psfile="split4a.ps" hscale=35 vscale=15 hoffset=-220 voffset=-120 }
\special{psfile="split4b.ps" hscale=35 vscale=15 hoffset= -50 voffset=-120 }
\vskip 1.7 in
\centerline{\it\quad (a) Splitting with McCormack, $t=5$     \hskip 0.2 in
(b) Splitting with McCormack, $t=10$}
\caption{Evolution of McCormack :
The mesh size is $\Delta x=0.1$.  Until $t<5$ the splitting method
(\ref{splitting}) with McCormack approximate the exact solution (line) of
(\ref{problem}) well, (a).  A small inadmissible shock appears around at
$t=5$ and it grows to tear down the solution, (b).}
\label{figmc1}
\end{figure}

The reason for this behavior is related to the fact that the discrete function
$$U^n_j=
\left\{
\begin{array}{rcl}
-1&,\quad j\le0,\\
 1&,\quad j>0,\\
\end{array}
\right.                  $$
is a steady state of the conservative scheme (\ref{step1}) with the numerical
flux of McCormack.
This example shows that this situation can be developed from general initial
data for a convection-reaction problems and splitting method let the phenomena
grow together with the solution. So it is natural to modify the first step
(\ref{step1}) to avoid such an inadmissible shock. An easy to do that is to
mix the information at adjacent mesh points. Even if we mix them like
\begin{equation}
U^{n+1}_j=e^k\Big((U^n_{j-1}+8U^n_j+U^n_{j+1})/10-{k\over h}\Delta F(U^n_j)\Big),
                                                \label{splitting1}\end{equation}
we can avoid such a blow-up, Figure \ref{figmc2} (a).
One interesting thing is that the solution does not have the oscilation
near the shock, which is the typical property of McCormack's scheme.
So the property of the scheme has been changed.

Now consider
\begin{equation}
U^{n+1}_j=e^k\Big((U^n_{j-1}+9998U^n_j+U^n_{j+1})/10000-{k\over h}\Delta F(U^n_j)\Big),
                                                \label{splitting2}\end{equation}
if the mixing is impact is too tiny, for example,


we may still have such a blow-up, Figure \ref{figmc2} (b).

\begin{figure}[ht]
\hskip .8 in
\special{psfile="split9a.ps" hscale=35 vscale=15 hoffset=-220 voffset=-120 }
\special{psfile="split9b.ps" hscale=35 vscale=15 hoffset= -50 voffset=-120 }
\vskip 1.7 in
\centerline{\it (a) Splitting method (\ref{splitting1}), $t=7$ \hskip 0.5 in
(b) Splitting Method (\ref{splitting2}), $t=7$}
\caption{Evolution of McCormack :
The mesh size is $\Delta x=0.1$. It is possible to remove the defect of the
McCormack's scheme by mixing grids, (a). But the}
\label{figmc2}
\end{figure}
There should be a quantitative guideline for the mixing.

We can also consider non-splitting method. First we consider the same way
we did for the Lax-Friedrichs,
\begin{equation}
U^{n+1}_j=U^n_j-{k\over h}\Delta F(U^n_j)+(e^k-1)U^n_j.
                                                \label{mcdirect1}\end{equation}
Unfortunately it does not work, Figure \ref{figmc3} (a).
We get the almost identical solution as Figure \ref{figmc1} (b).
It seem like that mixing the adjacent mesh points is still needed. In fact
a direct method,
\begin{equation}
U^{n+1}_j=U^n_j-{k\over h}\Delta F(U^n_j)+(e^k-1){U^n_{j-1}+U^n_{j+1}\over2},
                                                \label{mcdirect2}\end{equation}
gives a correct solution, Figure \ref{figmc3} (b). It is natural that the
way to neutralize the defect of the scheme can be different to each schemes.

\begin{figure}[ht]
\hskip .8 in
\special{psfile="split5a.ps" hscale=35 vscale=15 hoffset=-220 voffset=-120 }
\special{psfile="split5b.ps" hscale=35 vscale=15 hoffset= -50 voffset=-120 }
\vskip 1.7 in
\centerline{\it (a) Non-splitting method (\ref{mcdirect1}), $t=10$\hskip 0.2 in
(b) Non-splitting method (\ref{mcdirect2}), $t=10$}
\caption{Evolution of McCormack :
In this example the mesh size is $\Delta x=0.1$ and the initial data of
(\ref{problem}) is discretized by taking the cell average.}
\label{figmc3}
\end{figure}

\subsection{Godunov Method}
We consider the Godunov scheme together with interface method
for the source. First we introduce the interface of the Godunov method,
\begin{equation}
U_{j+1/2}=
\left\{
\begin{array}{ccc}
     U_{j+1}&,&f'(U_j)+f'(U_{j+1})\le0,\;f'(U_{j+1})\le0,\\
       U_{j}&,&f'(U_j)+f'(U_{j+1})  >0,\;f'(U_j)>0,\\
           0&,&f'(U_j)<0,\quad f'(U_{j+1})>0.\\
\end{array}
\right.
                                                  \label{iface}\end{equation}
The flux of the Godunov scheme is simply
\begin{equation}
F(U_j,U_{j+1})=f(U_{j+1/2}).
                                                  \label{godunov}\end{equation}
Interface method is
\begin{equation}
U^{n+1}_j=U^n_j-{k\over h}\Delta F(U^n_j)
                             +(e^k-1){U^n_{j-1/2}+U^n_{j+1/2}\over2}.
                                                  \label{imethod}\end{equation}

Figure \ref{figgo1} shows that numerical solutions from these schemes are
stable.
\begin{figure}[ht]
\hskip .8 in
\special{psfile="split6a.ps" hscale=35 vscale=15 hoffset=-220 voffset=-120 }
\special{psfile="split6b.ps" hscale=35 vscale=15 hoffset= -50 voffset=-120 }
\vskip 1.7 in
\centerline{\it (a) Splitting method (\ref{godunov},\ref{splitting}) at $t=10$     \hskip 0.2 in
(b) Interface method (\ref{imethod}) at $t=10$}
\caption{Evolution of Godunov :
In this example the mesh size is $\Delta x=0.1$ and the initial data of
(\ref{problem}) is discretized by taking the cell average.}
\label{figgo1}
\end{figure}

A challenging problem is to find the unstable steady state (\ref{rollwave}).
Consider a periodic boundary value problem,
\begin{equation}
\begin{array}{c}
u_t+uu_x=u,\quad -1<x<1, 0<t\\
u(x,0)={\bfR  }_1(x),\; u(-1,t)=u(1,t).\\
\end{array}
                                                 \label{problem2}\end{equation}
Since the roll wave ${\bfR  }_1(x)$ is a steady state, the exact solution
of (\ref{problem2}) is $u(x,t)={\bfR  }_1(x)$. In Figure \ref{figgo2} numerical
solutions for the splitting method and interface methods are displayed.

\begin{figure}[ht]
\hskip .8 in
\special{psfile="split7a.ps" hscale=35 vscale=15 hoffset=-220 voffset=-120 }
\special{psfile="split7b.ps" hscale=35 vscale=15 hoffset= -50 voffset=-120 }
\vskip 1.7 in
\centerline{\it (a) Splitting method (\ref{godunov},\ref{splitting}) at $t=37$
   \hskip 0.2 in
(b) Interface method (\ref{imethod}) at $t=100$}
\caption{Evolution of Godunov :
In this example the mesh size is $\Delta x=0.05$ and the initial data of
(\ref{problem2}) is discretized by taking the cell average.}
\label{figgo2}
\end{figure}

For the numerical solutions of the splitting method the structure of
the roll wave breaks down around at $t=35$ and it represents a wrong one
at $t=37$, (a).  The interface method provide a stable
solution, (b). So if the objective of the computation is to get that kind
of unstable steady states, interface method is an excellent choice.

On the other hand it means that the scheme has a property to neutralize the
instability of the problem (\ref{problem2}) and, hence, the scheme possibly
gives incorrect evolution when initial data consists of a small perturbation
of roll waves (\ref{rollwave}). As an example we consider a periodic problem,
\begin{equation}
\begin{array}{c}
u_t+uu_x=u,\qquad \; u(-1,t)=u(1,t),\\
u(x,0)=
\left\{
\begin{array}{cc}
x+1,&-1<x<0,\\
c(x-1),&0<x<1.\\
\end{array}
\right.
\end{array}
                                                 \label{problem3}\end{equation}
For $c\ne1$ the solution diverges as $t\to\infty$. Figure \ref{figgo3} shows
that interface method method (\ref{imethod}) gives a roll wave type solution
for $c=1.04$. For $t=1.05$ the solution diverges.



\begin{figure}[ht]
\hskip .8 in
\special{psfile="split8a.ps" hscale=35 vscale=15 hoffset=-220 voffset=-120 }
\special{psfile="split8b.ps" hscale=35 vscale=15 hoffset= -50 voffset=-120 }
\vskip 1.7 in
\centerline{\it (a) I-method for (\ref{problem3}), $t=6.5,c=1.5$
   \hskip 0.2 in
(b) I-method for (\ref{problem3}), $t=100,c=1.4$ }
\caption{Evolution of Interface method :
In this example the mesh size is $\Delta x=0.05$ and the initial data of
(\ref{problem3}) is discretized by taking the cell average.}
\label{figgo3}
\end{figure}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The splitting method simply magnify the error caused by numerical schemes
for the convection term and it causes non-physical solutions. So to avoid
these non-physical solutions we can consider two ways to handle it.
First we can correct the numerical schemes for the convection term
which caused the problem. Fro example if we use the intermediate value
$u_j^n$ in the Lax-Friedrichs scheme, we can avoid the oscillation.

Another way is to handle the reaction part differently. Step (2.21) simply
magnify the error produced in Step (2.15). So we may consider other way to
absorb or distribute the error. For example we can consider some non-splitting
method.

It is very interesting.  If the equation is $u_t =
(u^m)_xx+(u^m)_yy+(u^m)_zz$ in 3-D, then the similarity space
variables are ${x\over\root{3m-1}\of{t}},{y\over\root{3m-1}\of{t}},
{z\over\root{3m-1}\of{t}}$ and $n-2 < nm < n-1$ implies
$0<3m-1<1$.  So the decay rate looks like $O(t^{-3\over 3m-1})$.
Probably it is not that simple.


\begin{thebibliography}{10}

\bibitem{JinKatsoulakis} {\sc S. Jin and M.~A. Katsoulakis}
{\em  Hyperbolic Systems with Supercharacteristic Relaxations and Roll Waves},
SIAM J. Appl. Math., 61 (2000), pp.~273--292.

\bibitem{JinKim} {\sc S. Jin and Y.-J. Kim},
{\em On the computations of roll waves},
Math. Model. Num. Anal., to appear.

\bibitem{KimTzavaras} {\sc Y.-J. Kim and A.~E. Tzavaras},
{\em Diffusive N-waves and Metastability in Burgers equation},
SIAM J. Math. Anal. to appear
\bibitem{LevequeYee} {\sc R.~J. LeVeque and H.~C. Yee}
{\em A study of numerical methods for hyperbolic conservation laws with stiff source terms},
J. Comput. Phys., 86 (1990), pp.~187--210. 



\end{thebibliography}

\end{document}



















~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\begin{abstract}
This paper is a sample prepared to illustrate the use of the American
Mathematical Society's \LaTeX{} document class \texttt{amsart} and
publication-specific variants of that class for AMS-\LaTeX{} version 1.2.
\end{abstract}

\maketitle

\section*{This is an unnumbered first-level section head}
This is an example of an unnumbered first-level heading.

%% The correct journal style for \specialsection is all uppercase; a known bug
%% in amsart.cls prevents this, so input must be uppercase until it is fixed.
%\specialsection*{This is a Special Section Head}
\specialsection*{THIS IS A SPECIAL SECTION HEAD}
This is an example of a special section head%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\footnote{Here is an example of a footnote. Notice that this footnote
text is running on so that it can stand as an example of how a footnote
with separate paragraphs should be written.
\par
And here is the beginning of the second paragraph.}%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
.

\section{This is a numbered first-level section head}
This is an example of a numbered first-level heading.

\subsection{This is a numbered second-level section head}
This is an example of a numbered second-level heading.

\subsection*{This is an unnumbered second-level section head}
This is an example of an unnumbered second-level heading.

\subsubsection{This is a numbered third-level section head}
This is an example of a numbered third-level heading.

\subsubsection*{This is an unnumbered third-level section head}
This is an example of an unnumbered third-level heading.

\begin{lemma}
Let $f, g\in  A(X)$ and let $E$, $F$ be cozero
sets in $X$.
\begin{enumerate}
\item If $f$ is $E$-regular and $F\subseteq E$, then $f$ is $F$-regular.

\item If $f$ is $E$-regular and $F$-regular, then $f$ is $E\cup
F$-regular.

\item If $f(x)\ge c>0$ for all $x\in E$, then $f$ is $E$-regular.

\end{enumerate}
\end{lemma}

The following is an example of a proof.

\begin{proof} Set $j(\nu)=\max(I\backslash a(\nu))-1$. Then we have
\[
\sum_{i\notin a(\nu)}t_i\sim t_{j(\nu)+1}
  =\prod^{j(\nu)}_{j=0}(t_{j+1}/t_j).
\]
Hence we have
\begin{equation}
\begin{split}
\prod_\nu\biggl(\sum_{i\notin
  a(\nu)}t_i\biggr)^{\abs{a(\nu-1)}-\abs{a(\nu)}}
&\sim\prod_\nu\prod^{j(\nu)}_{j=0}
  (t_{j+1}/t_j)^{\abs{a(\nu-1)}-\abs{a(\nu)}}\\
&=\prod_{j\ge 0}(t_{j+1}/t_j)^{
  \sum_{j(\nu)\ge j}(\abs{a(\nu-1)}-\abs{a(\nu)})}.
\end{split}
\end{equation}
By definition, we have $a(\nu(j))\supset c(j)$. Hence, $\abs{c(j)}=n-j$
implies (5.4). If $c(j)\notin a$, $a(\nu(j))c(j)$ and hence
we have (5.5).
\end{proof}

\begin{quotation}
This is an example of an `extract'. The magnetization $M_0$ of the Ising
model is related to the local state probability $P(a):M_0=P(1)-P(-1)$.
The equivalences are shown in Table~\ref{eqtable}.
\end{quotation}

\begin{table}[ht]
\caption{}\label{eqtable}
\renewcommand\arraystretch{1.5}
\noindent\[
\begin{array}{|c|c|c|}
\hline
&{-\infty}&{+\infty}\\
\hline
{f_+(x,k)}&e^{\sqrt{-1}kx}+s_{12}(k)e^{-\sqrt{-1}kx}&s_{11}(k)e^
{\sqrt{-1}kx}\\
\hline
{f_-(x,k)}&s_{22}(k)e^{-\sqrt{-1}kx}&e^{-\sqrt{-1}kx}+s_{21}(k)e^{\sqrt
{-1}kx}\\
\hline
\end{array}
\]
\end{table}

\begin{definition}
This is an example of a `definition' element.
For $f\in A(X)$, we define
\begin{equation}
\mathcal{Z} (f)=\{E\in Z[X]: \text{$f$ is $E^c$-regular}\}.
\end{equation}
\end{definition}

\begin{remark}
This is an example of a `remark' element.
For $f\in A(X)$, we define
\begin{equation}
\mathcal{Z} (f)=\{E\in Z[X]: \text{$f$ is $E^c$-regular}\}.
\end{equation}
\end{remark}

\begin{example}
This is an example of an `example' element.
For $f\in A(X)$, we define
\begin{equation}
\mathcal{Z} (f)=\{E\in Z[X]: \text{$f$ is $E^c$-regular}\}.
\end{equation}
\end{example}

\begin{xca}
This is an example of the \texttt{xca} environment. This environment is
used for exercises which occur within a section.
\end{xca}

The following is an example of a numbered list.

\begin{enumerate}
\item First item.
In the case where in $G$ there is a sequence of subgroups
\[
G = G_0, G_1, G_2, \dots, G_k = e
\]
such that each is an invariant subgroup of $G_i$.

\item Second item.
Its action on an arbitrary element $X = \lambda^\alpha X_\alpha$ has the
form
\begin{equation}\label{eq:action}
[e^\alpha X_\alpha, X] = e^\alpha \lambda^\beta
[X_\alpha X_\beta] = e^\alpha c^\gamma_{\alpha \beta}
 \lambda^\beta X_\gamma,
\end{equation}

\begin{enumerate}
\item First subitem.
\[
- 2\psi_2(e) =  c_{\alpha \gamma}^\delta c_{\beta \delta}^\gamma
e^\alpha e^\beta.
\]

\item Second subitem.
\begin{enumerate}
\item First subsubitem.
In the case where in $G$ there is a sequence of subgroups
\[
G = G_0, G_1, G_2, \ldots, G_k = e
\]
such that each subgroup $G_{i+1}$ is an invariant subgroup of $G_i$ and
each quotient group $G_{i+1}/G_{i}$ is abelian, the group $G$ is called
\textit{solvable}.

\item Second subsubitem.
\end{enumerate}
\item Third subitem.
\end{enumerate}
\item Third item.
\end{enumerate}

Here is an example of a cite. See \cite{A}.

\begin{theorem}
This is an example of a theorem.
\end{theorem}

\begin{theorem}[Marcus Theorem]
This is an example of a theorem with a parenthetical note in the
heading.
\end{theorem}

\begin{figure}[tb]
\blankbox{.6\columnwidth}{5pc}
\caption{This is an example of a figure caption with text.}
\label{firstfig}
\end{figure}

\begin{figure}[tb]
\blankbox{.75\columnwidth}{3pc}
\caption{}\label{otherfig}
\end{figure}

\section{Some more list types}
This is an example of a bulleted list.

\begin{itemize}
\item $\mathcal{J}_g$ of dimension $3g-3$;
\item $\mathcal{E}^2_g=\{$Pryms of double covers of $C=\openbox$ with
normalization of $C$ hyperelliptic of genus $g-1\}$ of dimension $2g$;
\item $\mathcal{E}^2_{1,g-1}=\{$Pryms of double covers of
$C=\openbox^H_{P^1}$ with $H$ hyperelliptic of genus $g-2\}$ of
dimension $2g-1$;
\item $\mathcal{P}^2_{t,g-t}$ for $2\le t\le g/2=\{$Pryms of double
covers of $C=\openbox^{C'}_{C''}$ with $g(C')=t-1$ and $g(C'')=g-t-1\}$
of dimension $3g-4$.
\end{itemize}

This is an example of a `description' list.

\begin{description}
\item[Zero case] $\rho(\Phi) = \{0\}$.

\item[Rational case] $\rho(\Phi) \ne \{0\}$ and $\rho(\Phi)$ is
contained in a line through $0$ with rational slope.

\item[Irrational case] $\rho(\Phi) \ne \{0\}$ and $\rho(\Phi)$ is
contained in a line through $0$ with irrational slope.
\end{description}


\bibliographystyle{amsplain}
\begin{thebibliography}{10}

\bibitem {A} T. Aoki, \textit{Calcul exponentiel des op\'erateurs
microdifferentiels d'ordre infini.} I, Ann. Inst. Fourier (Grenoble)
\textbf{33} (1983), 227--250.

\bibitem {B} R. Brown, \textit{On a conjecture of Dirichlet},
Amer. Math. Soc., Providence, RI, 1993.

\bibitem {D} R. A. DeVore, \textit{Approximation of functions},
Proc. Sympos. Appl. Math., vol. 36,
Amer. Math. Soc., Providence, RI, 1986, pp. 34--56.

\end{thebibliography}

\end{document}

%------------------------------------------------------------------------------
% End of journal.tex
%------------------------------------------------------------------------------