The cutoff method for the numerical computation of nonnegative solutions of parabolic PDEs with application to anisotropic diffusion and lubrication-type equations
Changna Lu\thanks{College of Math and Statistics
Nanjing
Jiangsu 210044
Weizhang Huang\thanks{ Department of Mathematics
Lawrence
KS 66045
Erik S. Van Vleck\thanks{ Department of Mathematics
The cutoff method, which cuts off the values of a function less than a given number, is studied for the numerical computation of nonnegative solutions of parabolic partial differential equations. A convergence analysis is given for a broad class of finite difference methods combined with cutoff for linear parabolic equations. Two applications are investigated, linear anisotropic diffusion problems satisfying the setting of the convergence analysis and nonlinear lubrication-type equations for which it is unclear if the convergence analysis applies. The numerical results are shown to be consistent with the theory and in good agreement with existing results in the literature. The convergence analysis and applications demonstrate that the cutoff method is an effective tool for use in the computation of nonnegative solutions. Cutoff can also be used with other discretization methods such as collocation, finite volume, finite element, and spectral methods and for the computation of positive solutions.
lhvv2012-revise.tex
\documentclass[DIV=14,letterpaper]{scrartcl}
%\documentclass[11pt]{article}
%\usepackage{fullpage}
%\textheight=9.0in
%\textwidth=6.5in
%\baselineskip=12pt
%%\oddsidemargin=16pt
%%\evensidemargin=16pt
%\oddsidemargin=0pt
%\evensidemargin=0pt
%\topmargin=-0.5in

\date{}

\usepackage{tikz}
\usetikzlibrary{calc}
\usetikzlibrary{matrix}
\usepackage{amsmath,amssymb,amsfonts,amsthm,subfig}
%\usepackage{amsmath,amssymb,amsfonts} % Typical maths resource packages
%\usepackage{color}                    % For creating coloured text and background
%\usepackage{hyperref}                 % For creating hyperlinks in cross references
\usepackage[pdftex]{hyperref} 

\usepackage{graphicx}
%\usepackage{showlabels}
%\usepackage{epsf}

\renewcommand{\baselinestretch}{1.1}

\newcommand{\V}[1]{\mbox{\boldmath $ #1 $}}
\newcommand{\mpar}[1]{\marginpar{\tiny\sloppy #1}}
%\newcommand{\mpar}[1]{\marginpar{}}

\def \lapprox{\stackrel{<}{\sim}}
\def \p{\partial}

\newcommand{\bey}{\begin{eqnarray}}
\newcommand{\eey}{\end{eqnarray}}
\newcommand{\nn}{\nonumber}
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newtheorem{thm}{\hspace{6mm}Theorem}[section]
\newtheorem{pro}{\hspace{6mm}Proposition}[section]
\newtheorem{lem}{\hspace{6mm}Lemma}[section]
\newtheorem{co}{\hspace{6mm}Corollary}[section]
\newtheorem{de}{\hspace{6mm}Definition}[section]
\newtheorem{exam}{\hspace{6mm}Example}[section]
\newtheorem{rem}{\hspace{6mm}Remark}[section]
\newcommand{\proofend}{\mbox{ }\hfill \raisebox{.4ex}{\framebox[1ex]{}}}
%

\title{The cutoff method for the numerical computation of nonnegative solutions of
parabolic PDEs with application to anisotropic diffusion and lubrication-type equations}

\vskip 4mm

\author{Changna Lu\thanks{College of Math and Statistics,
Nanjing University of Information Science and Technology,
Nanjing, Jiangsu 210044, China. ({\tt luchangna@nuist.edu.cn})}
\and Weizhang Huang\thanks{
Department of Mathematics, the University of Kansas, Lawrence, KS 66045,
U.S.A. ({\tt huang@math.ku.edu})}
\and Erik S. Van Vleck\thanks{
Department of Mathematics, the University of Kansas, Lawrence, KS 66045,
U.S.A. ({\tt evanvleck@math.ku.edu})}
}

\begin{document}
\vskip 1cm
\maketitle

\begin{abstract}
The cutoff method,  which cuts off the values of a function less than a given number,
is studied for the numerical computation of nonnegative solutions
of parabolic partial differential equations. A convergence analysis is given for a broad class of
finite difference methods combined with cutoff for linear parabolic equations.
Two applications are investigated, linear anisotropic diffusion problems satisfying the setting
of the convergence analysis and nonlinear lubrication-type equations
for which it is unclear if the convergence analysis applies.
The numerical results are shown to be consistent with the theory
and in good agreement with existing results in the literature.
The convergence analysis and applications demonstrate that the cutoff method is an effective
tool for use in the computation of nonnegative solutions. Cutoff can also be used with other
discretization methods such as collocation, finite volume, finite element, and spectral
methods and for the computation of positive solutions.
\end{abstract}

\noindent{\bf AMS 2010 Mathematics Subject Classification.}
65M06, 35K51

\noindent{\bf Key Words.}
cutoff, projection, nonnegative solution, positive solution, error analysis,
finite difference, anisotropic diffusion, lubrication-type equation 

\noindent{\bf Abbreviated title.}
Cutoff method for computation of nonnegative solutions

%----------------------- section 1 ---------------------------
\section{Introduction}

Many physical phenomena involve variables such as the density and concentration of a material
that take only nonnegative values. The mathematical reflection of this property is that
the partial differential equations (PDEs) modeling the phenomena admit nonnegative
solutions. For the numerical solution of those PDEs, it is crucial that numerical schemes
preserve the solution nonnegativity and produce physically meaningful numerical solutions.

A closely related concept is the maximum principle. Preserving the maximum principle
is equivalent to preserving the solution nonnegativity for linear problems \cite{TW08} but generally
the former is more difficult than the latter. It is known (e.g., see \cite{CR73}) that a conventional numerical method,
such as a finite difference (FD), finite volume (FV), or finite element (FE) method,
generally does not preserve the maximum principle and can produce negative undershoot in the solution
for diffusion problems, especially those with heterogeneous anisotropic diffusion coefficients.
Considerable effort has been made on developing numerical schemes satisfying the maximum principle; 
e.g., see \cite{BKK08, BE04, Cia70, CR73, KK09, KKK07, KL95, Let92, Sto86, SF73, WaZh11, XZ99}
for steady-state isotropic diffusion problems and
\cite{DDS04, GL09, GYK05, Hua10, KSS09, LH10, LSS07, LSSV07, LS08, LuHuQi2012, MD06, SH07}
for steady-state anisotropic diffusion problems. Particularly, Ciarlet and Raviart \cite{CR73} show that the linear
FE approximation of an elliptic isotropic diffusion problem satisfies the maximum principle
when the mesh is simplicial and nonobtuse. The result is generalized to anisotropic diffusion problems
by Huang and his coworkers in \cite{Hua10,LH10, LuHuQi2012}.
On the other hand, less progress has been made for time dependent problems.
For example, Fujii \cite{Fuj73} shows that the linear FE approximation of the heat equation satisfies the maximum
principle when the mesh is simplicial and acute and the time step size is bounded above
by a bound proportional to the squared maximal element size and below by a bound
proportional to the squared minimal element in-diameter.
He also shows that when the lumped mass matrix is used, the maximum principle holds without requiring the time step
size to be bounded from below.
Fujii's results are extended to more general isotropic diffusion problems by Farag{\'o} et al. \cite{FH2001,FHK05,FHK2011} and to anisotropic diffusion problems by Li and Huang \cite{LiHu2012}.
Farag{\'o} and Horv{\'a}th \cite{FH06} study the relations between the maximum principle, nonnegativity
preservation, and maximum norm contractivity for linear parabolic equations and their finite difference
and Galerkin finite element discretizations.
Thom{\'e}e and Wahlbin \cite{TW08} consider more general parabolic
PDEs and show that
the maximum principle cannot hold for the conventional semidiscrete FE problem
and Fujii's conditions on the mesh are essentially sharp for the lumped mass matrix.
Le Potier \cite{LePot05,LePot09b} (also see Lipnikov et al. \cite{LSSV07})
proposes two nonlinear FV schemes for linear anisotropic diffusion problems and shows that they are
second order in space and satisfy monotonicity \cite{LePot05} or the maximum principle \cite{LePot09b}.

We consider initial-boundary value problems of parabolic PDEs which are well posed and have
a unique, nonnegative solution. Instead of employing a maximum-principle preserving scheme,
we propose to use the cutoff method that cuts off the negative values in the computed solution
at each time step and then continues the time integration with the corrected solution;
see Fig.~\ref{f-1} for a sketch of the solution procedure.
The cutoff method shares a similar idea with many projection and correction methods.
For example, projection methods \cite{Chorin1968,Teman1968}
are an effective means of enforcing the incompressibility condition in the numerical solution of time dependent
incompressible fluid flow problems. 
Projection is also used by some researchers in the numerical solution
of Hamiltonian systems to preserve the energy; e.g., see \cite{Hai2006}.
A FE/implicit Euler discretization of a Cahn-Hilliard equation with logarithmic
nonlinearity is considered in \cite{CE92} and iterative methods to solve the resulting
nonlinear systems of equations are developed to keep the solution within its physically
relevant range.
Solution compression is used in the design of maximum-principle-satisfying high order schemes
for scalar conservation laws by Zhang and Shu \cite{ZhSh2010}.
It should be emphasized that the cutoff method provides several advantages over
many positivity-preserving or maximum-principle-preserving schemes.
Its implementation is simple and requires no significant changes in the
existing code. Moreover, the cutoff procedure does not impose any constraint on the mesh and time step
which maximum-principle-preserving schemes often impose. 
However, unlike projection methods, there do not seem to exist published theoretical or numerical
studies on the cutoff method.

The objective of this paper is to provide a theoretical and numerical study of the cutoff method
for use in the computation of nonnegative solutions of parabolic PDEs. 
We shall first prove that a broad class of finite difference methods are convergent when they are incorporated with
the cutoff method for linear parabolic PDEs. Two applications are then investigated.
The first is a linear anisotropic diffusion problems which satisfies the setting of the convergence analysis.
It is known that a conventional finite difference or finite element method with a uniform or quasi-uniform mesh
does not preserve the maximum principle (nor the solution nonnegativity) and
typically produces negative undershoot in
the computed solution. The cutoff method removes the unphysical negative values in the solution
while keeping the same convergence order of the underlying scheme.
The second application is nonlinear lubrication-type equations. It is unclear if the convergence analysis
(which can be extended to some nonlinear equations; see Remark~\ref{rem2.0})
applies to those equations. A distinct feature of lubrication-type equations is that a positive solution can
develop a finite time singularity of the form $u \to 0$ and become identically zero on some spatial interval
for a period of time.  A conventional finite difference or finite element method can
produce negative values in the computed solution during this development of singularity.
Negative solution values are not only unphysical but also cease the computation
for typical situations where the nonlinear diffusion coefficients are undefined for negative solution values.
The application of the cutoff method avoids this difficulty. Moreover, it is shown that
no regularization is needed when cutoff is used. Numerical results are
shown to be in good agreement with existing ones in the literature.

The outline of the paper is as follows. The cutoff method is described and some of its properties are given in
Section \ref{SEC:cutoff}. The convergence of a broad class of finite difference methods is proved also in
Section \ref{SEC:cutoff} when they are incorporated with the cutoff method for the computation of
nonnegative solutions of linear parabolic PDEs. Application of finite difference methods with
cutoff to a linear anisotropic diffusion problem and a lubrication-type PDE is studied in Sections~\ref{SEC:aniso-diffusion}
and \ref{SEC:lubrication}, respectively. Finally, Section~\ref{SEC:conclusion} contains conclusions and 
further comments.

% f-1
\begin{figure}[tb]
\centering
\begin{tikzpicture}[scale = 1]
\draw [thick, ->] (0,1.5) -- (8,1.5);
\draw [thick, ->] (0,1.5) -- (0,6.5);
\draw[right] (8.1,1.5) node {$t$};
\draw[above] (0,6.6) node {$u$};

\draw[thick] (1,3.5) .. controls (3,5.5) and (6, 4.5) .. (7.5, 5.5);
\draw[right] (7.53, 5.55) node {$u=u(t)$};
\fill (canvas cs:x=1cm,y=3.5cm) circle (1.5pt);
\draw[above] (1,3.5) node {$u^0$};
\fill (canvas cs:x=3cm,y=4.63cm) circle (1.5pt);
\draw[above] (3,4.7) node {$u^n$};
\fill (canvas cs:x=4.5cm,y=4.88cm) circle (1.5pt);
\draw[above] (4.5,4.9) node {$u^{n+1}$};

%\draw[thick] (3,3) .. controls (4.5,3.5) and (6, 3.3) .. (7.2, 4);
%\draw[right] (7.25, 4.05) node {$u=\widetilde{u}(t)$};
\fill (canvas cs:x=3cm,y=3cm) circle (1.5pt);
\draw[left] (3.8,2.6) node {${U}^n$};

\draw[thick, dashed] (2, 2.7) -- (3,3);
\draw[thick, dashed] (3, 3.5) -- (3,3);
\draw[above] (3,3.5) node {$(U^n)^{+}$};
\fill (canvas cs:x=3cm,y=3.5cm) circle (1.5pt);

\fill (canvas cs:x=4.5cm,y=3.6cm) circle (1.5pt);
%\draw[help lines,->] (5.5, 4.5) -- (4.7,3.85);
%\draw[right] (5.5, 4.8) node {$({U}^{n+1})^{+}$};
\draw[below] (4.5, 3.4) node {${U}^{n+1}$};
\draw[thick, dashed] (4.5, 4) -- (5.5,4.2);

\fill (canvas cs:x=4.5cm,y=4cm) circle (1.5pt);
\draw[above] (4.5, 4) node {$(U^{n+1})^{+}$};
\draw [thick, dashed] (3, 3.5) -- (4.5, 3.65);
\draw [thick, dashed] (4.5,4) -- (4.5, 3.65);

\fill (canvas cs:x=1cm,y=1.5cm) circle (1.5pt);
\draw[below] (1,1.5) node {$t_0$};
\fill (canvas cs:x=3cm,y=1.5cm) circle (1.5pt);
\draw[below] (3,1.5) node {$t_n$};
\fill (canvas cs:x=4.5cm,y=1.5cm) circle (1.5pt);
\draw[below] (4.5,1.5) node {$t_{n+1}$};

\end{tikzpicture}
\caption{An illustration of the time integration with cutoff (the dashed line).
Here, $U^n$ is the numerical solution, $({U}^{n})^{+}$ is the corrected numerical solution and
$u(t)$ is the exact solution of the IBVP starting from $u^0$ at $t=t_0$.
%and $\widetilde{u}(t)$ is the exact solution of the IBVP starting from $({U}^{n})^{+}$ at $t=t_n$.
}
\label{f-1}
\end{figure}



% section 2
\section{The cutoff method and convergence analysis}
\label{SEC:cutoff}

In this section, we first describe the cutoff method and some of its properties.
We then give a convergence analysis of the method for a class of finite difference methods
applied to linear parabolic problems. Possible generalizations are also discussed.

% section 2.1
\subsection{The cutoff method}

Consider a continuous function $f = f(\V{x})$ defined on a domain $\Omega \subset \mathbb{R}^d$ ($d \ge 1$).
For any given cutoff parameter $\delta \in \mathbb{R}$, we define the $\delta$-cutoff function as
\beq
f_\delta^{+}(\V{x}) = 
\begin{cases}
f(\V{x}),& \text{ for } f(\V{x}) \ge \delta  \\
\delta, & \text{ for } f(\V{x}) \le \delta 
\end{cases}
\quad \forall \V{x}\,   \in \Omega .
\label{corect-2}
\eeq
Note that $f^{+}(\V{x}) \equiv f_0^{+}(\V{x})$ is the nonnegative part of $f$, i.e.,
\beq
f^{+}(\V{x}) = \frac{1}{2} (|f(\V{x})| + f(\V{x})) = 
\begin{cases}
f(\V{x}),& \text{ for } f(\V{x}) \ge 0 \\
0, & \text{ for } f(\V{x}) \le 0
\end{cases}
\quad \forall \V{x}\,   \in \Omega .
\label{corect-1}
\eeq

The following three lemmas give some properties of cutoff functions.
Although these properties are easy to prove (and almost obvious), they play a central role
in the convergence analysis for a class of finite difference methods in the next subsection.

\begin{lem}
\label{lem2.1}
For a given function $f$ and any nonnegative continuous function $u$ defined on $\Omega$, we have
\beq
|f^{+}(\V{x}) - u(\V{x})| \le |f(\V{x}) - u(\V{x})|,\quad \forall \, \V{x} \in \Omega.
\label{lem2.1-1}
\eeq
\end{lem} 

{\bf Proof.}
Obviously, (\ref{lem2.1-1}) holds when $f(\V{x}) \ge 0$. When $f(\V{x}) < 0$, we have $f^{+}(\V{x}) = 0$.
By assumption, $u(\V{x}) \ge 0$. Thus, we have
\[
| f^{+}(\V{x}) - u(\V{x}) | = | u(\V{x}) | \le  |f(\V{x})| + | u(\V{x}) | = |f(\V{x}) - u(\V{x})|.
\]
\proofend

\begin{lem}
\label{lem2.2}
For a given function $f$ and any nonnegative continuous function $u$ defined on $\Omega$, we have
\beq
|f^{+}(\V{x}) - f(\V{x})| \le |u(\V{x}) - f(\V{x}) |,\quad \forall \, \V{x} \in \Omega.
\label{lem2.2-1}
\eeq
\end{lem} 

{\bf Proof.} The proof is similar to that of Lemma~\ref{lem2.1}.
\proofend

\begin{lem}
\label{lem2.3}
Given a cutoff parameter $\delta \ge 0$, for any function $f$ and any nonnegative continuous function $u$
defined on $\Omega$, we have
\bey
&&  |f_{\delta}^{+}(\V{x}) - f^{+}(\V{x})| \le \delta, \quad \forall \, \V{x} \in \Omega
\label{lem2.3-1}
\\
&& |f_{\delta}^{+}(\V{x}) - u(\V{x})| \le |f(\V{x}) - u(\V{x})| + \delta,\quad \forall \, \V{x} \in \Omega.
\label{lem2.3-2}
\eey
\end{lem} 

{\bf Proof.} Inequality (\ref{lem2.3-1}) follows from the definitions of $f^{+}$ and $f_{\delta}^{+}$.
Inequality (\ref{lem2.3-2}) follows from the triangle inequality $|f_{\delta}^{+}(\V{x}) - u(\V{x})| 
\le |f^{+}(\V{x}) - u(\V{x})| + |f_{\delta}^{+}(\V{x}) - f^{+}(\V{x})|$ and Lemma~\ref{lem2.1}.
\proofend

% section 2.2
\subsection{Error analysis for a linear IBVP problem}
\label{SEC:2.2}

We now use the results in the previous subsection to analyze the convergence of 
a class of finite difference schemes for a general initial-boundary value problem (IBVP) in the form
\beq
\begin{cases}
\frac{\p u}{\p t} = L(u),\quad \mbox{in} \quad \Omega\times (t_0, T] \\
u = g,\quad \mbox{on} \quad \p\Omega \\
u = u^0, \quad \mbox{on} \quad \Omega \times\{ t = t_0\}
\end{cases}
\label{ibvp-1}
\eeq
where $\Omega \subset \mathbb{R}^d$ ($d \ge 1$) is the physical domain, $L$ is a linear elliptic
differential operator and $g$ and $u^0$ are given sufficiently smooth functions. We assume that the IBVP
is well posed and admits a unique nonnegative continuous solution, $u=u(t,\V{x})$. We also assume that
$L$ and $g$ do not contain $t$ explicitly for notational simplicity.

We consider a class of finite difference schemes for (\ref{ibvp-1})  in the matrix form as
\beq
B_1 U^{n+1} = B_0 U^n + F^n,
\label{fd-1}
\eeq
where $B_0$ and $B_1$ are matrices independent of $n$ and $U^n$ is an approximation of
the exact solution at $t=t_n$, i.e., $U^n = \{ U^n_j\},\; u^n = \{u_j^n\},\; 
U^n_j  \approx u_j^n \equiv u(t_n, \V{x}_j)$ for each mesh vertex $\V{x}_j$.
We assume that (\ref{fd-1}) satisfies
\bey
&& \| B_1^{-1} \| \le K_1,
\label{fd-2}
\\
&& \| B_1^{-1}B_0 \| \le (1 + K\Delta t),
\label{fd-3}
\\
&& \frac{1}{\Delta t} \left [ B_1 u^{n+1} - B_0 u^n - F^n \right ] \longrightarrow \frac{\p u}{\p t} - L(u), \quad
	\mbox{ as} \quad \Delta t(h) \to 0 
\label{fd-4}
\eey
where $\Delta t$ and $h$ are the maximal time step size and the maximal element size, respectively,
$K$ and $K_1$ are positive constants, and $\|\cdot \|$ is a proper matrix norm.
Condition (\ref{fd-2}) requires that (\ref{fd-1}) produce a bounded solution while
(\ref{fd-3}) and (\ref{fd-4}) are the stability and consistency conditions, respectively.
The local truncation error of this scheme is defined as
\beq
{\tau}^n = B_1 {u}^{n+1} - B_0 {u}^n - F^n .
\label{lte-1}
\eeq
Assume that scheme (\ref{fd-1}) is $(p,q)$-order for some positive integers $p$ and $q$.
Then, there exists a constant $C_{lte}(u)$ (depending only on
the exact solution) such that
\beq
\| {\tau}^n \| \le C_{lte} (u) \, \Delta t \, (\Delta t^p + h^q) .
\label{lte-2}
\eeq

It is remarked that there exist many schemes satisfying assumptions (\ref{fd-2}) -- (\ref{fd-4}).
Examples include those employing central finite differences
for spatial discretization and the $\theta$-method for temporal discretization;
e.g., see Morton~and~Mayers~\cite{MM05}.

\begin{thm}
\label{thm2.1}
Assume that IBVP (\ref{ibvp-1}) is well posed and admits a unique nonnegative continuous exact solution $u=u(t,\V{x})$.
We also assume that scheme (\ref{fd-1}) satisfies (\ref{fd-2}) -- (\ref{fd-4}) and (\ref{lte-2}).
Then, the error for the cutoff solution procedure shown in Fig.~\ref{f-1} with scheme (\ref{fd-1}) is bounded by
\beq
\|(U^n)^{+} - u^n\| \le \| U^{n} - u^n \| \le \frac{K_1}{K}\, e^{K t_n} \, C_{lte}(u)\, (\Delta t^p + h^q) .
\label{thm2.1-1}
\eeq
\end{thm}

{\bf Proof.} Let $e^n = U^n - u^n$. Notice that the cutoff solution procedure shown in Fig.~\ref{f-1}
with (\ref{fd-1}) satisfies
\beq
B_1 U^{n+1} = B_0 ({U}^n)^{+} + F^n.
\label{thm2.1-2}
\eeq
Combining this with (\ref{lte-1}), we obtain
\[
B_1 e^{n+1} = B_0 (({U}^n)^{+} - u^n) - {\tau}^n .
\]
It follows that
\[
\| e^{n+1}\| \le \|B_1^{-1} B_0\| \cdot \| ({U}^n)^{+} - u^n \|  +   \| B_1^{-1}\| \cdot \| {\tau}^n\| .
\]
Combining this with (\ref{fd-2}) and (\ref{fd-3}) leads to
\beq
\| e^{n+1}\| \le (1+K\Delta t) \, \| ({U}^n)^{+} - u^n \|  +   K_1 \, \| {\tau}^n\| .
\label{thm2.1-5}
\eeq
Lemma~\ref{lem2.1} implies
\beq
 \| ({U}^n)^{+} - u^n \|  \le \| {U}^n - u^n \| .
 \label{thm2.1-7}
\eeq
Thus, we get
\beq
\| e^{n+1}\| \le (1+K\Delta t) \, \| {e}^n \|  +   K_1 \, \| {\tau}^n\| .
\label{thm2.1-6}
\eeq
Then it is standard to show that (\ref{thm2.1-1}) follows from (\ref{lte-2}) and (\ref{thm2.1-6}).
\proofend

\vspace{10pt}

\begin{rem} {\em
\label{rem2.0}
From the above proof we can see that the key condition is the convergence of
the original scheme (without cutoff). If it (without cutoff) is convergent,
using Lemma~\ref{lem2.1} we can readily show that the scheme with cutoff is also convergent.
This observation implies that the convergence analysis of the cutoff method can be extended to
more general linear or nonlinear IBVPs.
}\proofend
\end{rem}

\vspace{10pt}

When a strictly positive solution is desired, we can replace $(U^n)^{+}$ with $(U^n)^{+}_\delta$ for some positive
constant $\delta$. The following theorem can be proven
using Lemma~\ref{lem2.3} in a similar way as for Theorem~\ref{thm2.1}.

\begin{thm}
\label{thm2.2}
Assume that IBVP (\ref{ibvp-1}) is well posed and admits a unique nonnegative continuous exact solution $u=u(t,\V{x})$.
We also assume that scheme (\ref{fd-1}) satisfies (\ref{fd-2}) -- (\ref{fd-4}) and and (\ref{lte-2}).
For the cutoff solution procedure shown in Fig.~\ref{f-1} with scheme (\ref{fd-1}) and with $(U^n)^{+}$ being replaced
with $(U^n)^{+}_\delta$ for some positive $\delta$, the error is bounded by
\beq
\|(U^n)^{+}_\delta  - u^n\| \le \frac{K_1}{K}\, e^{K t_n} \, C_{lte}(u)\, (\Delta t^p + h^q)
+ \left ( \frac{(1+K \Delta t)}{K \Delta t} \, e^{K t_n} \, + 1\right ) \delta .
\label{thm2.2-1}
\eeq
\end{thm}

\begin{rem}{\em 
\label{rem2.1}
Inequality (\ref{thm2.2-1}) shows that $\delta$ should be chosen proportional to $\Delta t$
for the error bound to stay bounded as $\Delta t \to 0$. Ideally, $\delta$ should be
at the same level as the local truncation error, i.e.,
\beq
\delta = \mathcal{O}( \Delta t (\Delta t^p + h^q)) .
\label{thm2.2-2}
\eeq
This way, the terms on the right-hand side of (\ref{thm2.2-1}) have the same convergence order
as $\Delta t (h) \to 0$.
}\proofend
\end{rem}

%% section 2.3
%\subsection{Generalizations}

%
%It should also be pointed out that the setting in the previous subsection is basically for finite difference and
%finite volume methods. The procedure can be extended to finite element methods without major difficulty.

%% section 2.4
%\subsection{A conjecture}
%We conjecture that the difference between the the corrected numerical solution and the numerical
%solution is at the level of the local truncation error, i.e.,
%\beq
%\| (U^n)^{+} - U^n \| = \mathcal{O} (\| \tau^n \|) = \mathcal{O}(\Delta t (\Delta t^p + h^q) ).
%\label{diff-1}
%\eeq

%To see this, we can $\widetilde{u}(t,\V{x})$ (see Fig.~\ref{f-1}) for $u(\V{x})$ in Lemma~\ref{lem2.2} and get
%\[
%\| (U^{n+1})^{+} - U^{n+1} \| \le \| U^{n+1} - \widetilde{u}^{n+1} \|.
%\]
%Notice that $U^{n+1} - \widetilde{u}^{n+1}$ is essentially the local truncation error evaluated using
%$\widetilde{u}(t,\V{x})$. Thus, we have
%\[
%\| U^{n+1} - \widetilde{u}^{n+1} \| = \mathcal{O}( C_{lte}(\widetilde{u}) \Delta t (\Delta t^p + h^q) ).
%\]
%By the well-posedness (more specifically, the continuous dependence of the solution on the initial solution)
%of the continuous problem, we may expect
%\[
%C_{lte}(\widetilde{u}) = \mathcal{O} (C_{lte}({u})).
%\]

%{\color{red} The above argument is valid only if the numerical solution $(U^n)^{+}$ (and thus $\widetilde{u}(t,\V{x})$)
%are sufficiently smooth in space (which is not necessarily true). Any idea for proving (\ref{diff-1}) will be great!}

% section 3
\section{An anisotropic diffusion problem}
\label{SEC:aniso-diffusion}

In this section we present numerical results obtained for a 2D linear anisotropic diffusion problem
by the cutoff method described in the previous section.  The problem takes the form of 
IBVP (\ref{ibvp-1}) with
\beq
L(u) = \nabla \cdot (\mathbb{D} \, \nabla u) + f, \quad \Omega = [0,1]\times [0,1],\quad
\mathbb{D} =  \left ( \begin{array}{cc} 500.5 & 480 \\ 480 & 500.5  \end{array} \right ),
\label{ibvp-2}
\eeq
and $f$, $u^0$, and $g$ are chosen such that the exact solution of the IBVP is given by
\beq
u=\frac{1}{2} \mbox{exp}(-t) ( \tanh(-15( x-y )) +1 ).
\label{exam3.1-soln}
\eeq
The problem satisfies the maximum principle and the solution stays between 0 and 1.

It is worth pointing out that unlike lubrication-type equations which we shall consider
in the next section, for this problem the computation can continue when negative values
occur in the computed solution. From this point of view, it is unnecessary to remove negative
values at each time step. Nevertheless, this problem satisfies the setting of the convergence analysis
in the previous section and provides a good example for verifying the theory and testing
the effectiveness and accuracy of the cutoff method.

We use Cartesian grids of size $J\times J$ for the physical domain,
central finite differences for spatial discretization, and a third-order
singly diagonally implicit Runge-Kutta  (SDIRK) method \cite{Ale77,Cas79}
for temporal discretization of the underlying PDE. The discretization is standard and can be shown
to be convergent with order $(3,2)$ (third order in time and second order in space).
It is also known (and is shown below) that this standard scheme does not preserve the maximum principle
and produces spurious undershoot
and overshoot in the numerical solution. The numerical results presented below are obtained with
a fixed time step size $\Delta t = 10^{-2}$, which was found to be small enough so that the temporal
discretization error is ignorable compared to the spatial discretization error. 

Fig.~\ref{Exa3.1-1} shows the surface plot of a numerical solution (after cutoff) at $t=1$.
The contours of the numerical solutions before and after cutoff are shown in Fig.~\ref{Exa3.1-2} (a) and (b).
One may notice that the underlying finite difference scheme does not preserve the nonnegativity
and the numerical solution at $t=1$ contains negative values (with $U_{min}^N=-1.073\times 10^{-3}$
where $N$ is the last time step) before cutoff.

The $L^2$ error, $\| (U^N)^{+}-u^N \|_{L^2(\Omega)}$, and the maximal undershoot, $-u_{min}$,
are shown in Fig.~\ref{Exa3.1-3} as functions of the number of subintervals in $x$ (or $y$) direction. 
We can see that the $L^2$ error decreases at a rate of $O(J^{-2})$ (second order
in terms of element size) as $J$ increases. This is consistent with the theoretical prediction given
in Theorem~\ref{thm2.1}. On the other hand, the maximal undershoot, which is equal to the cutoff error, i.e.,
$-u_{min} = \| U^N - (U^N)^+\|_\infty$, decreases at a faster rate. This is somewhat surprising since we expect
the cutoff error to be at the level of the local truncation error, which is second order in space.

Recall that the cutoff strategy can be used to preserve the positivity of the solution. For example,
when the cutoff parameter is taken $\delta =  \Delta t ( \Delta x)^2$,
we have $(U^n_j)^{+}_{\delta} \ge \Delta t ( \Delta x)^2$. Theorem~\ref{thm2.2} implies that
the $L^2$ error, $\| (U^N)^{+}_{\delta}-u^N \|_{L^2(\Omega)}$, decreases at a rate of second order in space.
For the current example, the numerical results obtained with this cutoff parameter are indistinguishable
from those shown in Figs.~\ref{Exa3.1-2} and \ref{Exa3.1-3}. For this reason and to save space, we omit
those results here.

It is worth mentioning that computations were also performed for  
a similar problem with a convection term,
\beq
L(u) = \nabla \cdot (\mathbb{D} \, \nabla u) - \V { b } \cdot \nabla u + f, 
\label{ibvp-3}
\eeq
where $\mathbb{D}$ is given in (\ref{ibvp-2}), $\V{b} = [1000,1000]^T$, and the functions
$f$, $g$, and $u^0$ are chosen such that the exact solution is given by (\ref{exam3.1-soln}).
The obtained results (not shown) are very similar to those shown in Figs. \ref{Exa3.1-2} and \ref{Exa3.1-3}.

%Fig.~\ref{Exa3.2-2} shows the contours of the numerical solutions at $t=1$ before and after cutoff.
%It can be seen that the solution before the cutoff contains negative values (with
%$u_{min}=-1.102\times 10^{-3}$). Fig.~\ref{Exa3.2-3} shows the second order decrease rate
%of the $L^2$ norm of the error, which is consistent with the theoretical prediction (cf. Theorem~\ref{thm2.1}).
%Once again, the undershoot or the cutoff error decreases at a rate of higher than second order as
%$J=K$ increases.


%Exa3.1-1
\begin{figure}[thb]
\centering
%\centerline{(b): A numerical solution}
\includegraphics[height=3.5in]{ex1-solutions81-color.pdf}
\caption{Example~(\ref{ibvp-2}). The surface plot of a numerical solution (after cutoff) at $t=1$ obtained
with an $80 \times 80$ Cartesian grid.}
\label{Exa3.1-1}
\end{figure}

%Exa3.1-2
\begin{figure}[thb]
\centering
\hbox{
\begin{minipage}[b]{3in}
\centerline{(a): $U^N$}
\includegraphics[width=3.5in]{ex1-contour81-before-color.pdf}
\end{minipage}
\begin{minipage}[b]{3in}
\centerline{(b): $(U^N)^{+}$}
\includegraphics[width=3.5in]{ex1-contour81-after-color.pdf}
\end{minipage}
}
\caption{Example~(\ref{ibvp-2}). Contours of the numerical solutions at $t=1$ before and after cutoff.
The solutions are obtained with an $80 \times 80$ Cartesian grid.}
\label{Exa3.1-2}
\end{figure}

%Exa3.1-3
\begin{figure}[thb]
\centering
\includegraphics[scale=0.8]{ex1-error-umin-color.pdf}
\caption{Example~(\ref{ibvp-2}).
The $L^2$ error, $\| (U^N)^{+}-u^N \|_{L^2(\Omega)}$, and the maximal undershoot, $-u_{min}$,
are shown as functions of the number of subintervals in $x$ (or $y$) direction.
}
\label{Exa3.1-3}
\end{figure}

%%Exa3.1-4
%\begin{figure}[thb]
%\centering
%\includegraphics[scale=1]{ex1-error-umin-positive-color.pdf}
%\caption{Example~\ref{exam3.1}.
%The $L^2$  of the positive solution error and undershoot $-u_{min}$
%are shown as functions of the number of elements, the minimum positive solution is $\Delta t (\Delta x)^2$. 
%The solutions are obtained with a uniform mesh with $J \times J$ cells.}
%\label{Exa3.1-4}
%\end{figure}



%%Exa3.2-1
%\begin{figure}[thb]
%\centering
%\hbox{
%\hspace{2.5mm}
%\begin{minipage}[b]{3in}
%\centerline{(a): $N=2350$}
%\includegraphics[height=2in]{ex2-mesh2350.pdf}
%\end{minipage}
%\begin{minipage}[b]{3in}
%\centerline{(b): Numerical solutions}
%\includegraphics[height=2in]{ex2-solution2350.pdf}
%\end{minipage}
%}
%\caption{Example 3.2. Numerical solutions at the end time $t=1.0s$ obtained with  uniform meshes (a).}
%\label{Exa3.2-1}
%\end{figure}

%%Exa3.2-2
%\begin{figure}[thb]
%\centering
%\hbox{
%\begin{minipage}[b]{3in}
%\centerline{(a): $U^N$}
%\includegraphics[width=3.5in]{ex2-contour81-before.pdf}
%\end{minipage}
%\begin{minipage}[b]{3in}
%\centerline{(b): $(U^N)^{+}$}
%\includegraphics[width=3.5in]{ex2-contour81-after.pdf}
%\end{minipage}
%}
%\caption{Example~\ref{exam3.2}.
%Contours of the numerical solutions at $t=1$ before and after cutoff.
%The solutions are obtained with an $80 \times 80$ Cartesian grid.}
%\label{Exa3.2-2}
%\end{figure}

%%Exa3.2-3
%\begin{figure}[thb]
%\centering
%\includegraphics[scale=1]{ex2-error-umin-color.pdf}
%\caption{Example~\ref{exam3.2}.
%The $L^2$ error $\| (U^N)^{+}-u^N \|_{L^2(\Omega)}$ and the maximal undershoot $-u_{min}$
%are shown as functions of the number of subintervals in $x$ (or $y$) direction.
%}
%\label{Exa3.2-3}
%\end{figure}


% section 4
\section{Application to lubrication-type equations}
\label{SEC:lubrication}

In this section we consider the application of the cutoff method to a lubrication-type equation, which was first derived
from lubrication theory by Greenspan  \cite{Gre1978} for describing the movement of thin viscous films and
spreading droplets. Lubrication theory consists of a depth-averaged equation of mass conservation and a simplified
form of the Navier-Stokes equations that is appropriate for a thin layer of very viscous fluid.
We consider a lubrication-type equation in the general form
\beq
u_t + \nabla \cdot (f(u) \nabla \Delta u ) = 0,\quad f(u)\sim u^n, \;\; n \in [0, \infty )
\label{lub-1}
\eeq
where $u$ is the thickness of the viscous droplet and $n$ is a physical parameter which
has different values for different boundary conditions on the liquid solid interface.
Lubrication-type equations also appear in several other applications, including a thin neck model
in the Hele-Shaw cell \cite{CDGKS1993},  Cahn-Hilliard models  with degenerate mobility \cite{EllGar1996},
and problems in population dynamics \cite{Lew1994} and plasticity \cite{Grun1995}. 
It is noted that the solution $u$ must stay nonnegative to be physically meaningful in all of those applications.

It is known theoretically \cite{BeFr1990,BBDK1994} that in one dimension, (\ref{lub-1}) 
preserves the positivity of the solution for $n \ge 3.5$ and has a nonnegative weak solution
in a sense of distributions for $3/8 < n < 3$ (and in a weaker sense for $0 < n < 3/8$).
Such a weak solution can be obtained \cite{BeFr1990} as the limit as $\epsilon \to 0$ of the solution
of the regularized problem
\beq
u_t^\epsilon + \nabla \cdot (f^\epsilon(u^\epsilon) \nabla \Delta u^\epsilon ) = 0,\quad
f^\epsilon(u^\epsilon) = \frac{(u^\epsilon)^4 f(u^\epsilon)}{\epsilon f(u^\epsilon) + (u^\epsilon)^4} ,
\label{lub-2}
\eeq
where $\epsilon > 0$ is the regularization parameter. Moreover, numerical computation (such as see
\cite{Ber1995, Ber1996, BBDK1994}) shows that for small values of $n > 0$, a positive solution
can first develop a finite time singularity of the form $u \to 0$ at some point (which physically describes
the rupture of the liquid film), then becomes identically zero on an interval of time dependent length for a period of time,
becomes positive again at a later time, and eventually decays to the mean of the initial solution.
It is challenging to simulate this singularity development
since conventional numerical methods do not preserve the nonnegativity of the solution and requires
a huge number of mesh points to provide the necessary resolution to avoid spurious, negative undershoot.
Moreover, when $n$ is an even denominator rational number,  $f(u)$ is undefined for negative $u$.
In this situation, the computation typically stops around the initial formation of the singularity when negative values
start to appear in the numerical solution. For this reason, the simulation of the singularity development in (\ref{lub-1})
is a good test for the cutoff method although it is unclear if the convergence analysis described
in \S~\ref{SEC:2.2} applies to (\ref{lub-1}).

Great effort has been made in the numerical solution of (\ref{lub-1}).
In addition to the above mentioned references \cite{Ber1995, Ber1996, BBDK1994}, 
Zhornitskaya and Bertozzi \cite{ZhoBer2000} propose a positivity preserving finite difference scheme
for the regularized equation (\ref{lub-2}). Russell et al. \cite{RuWiXu2007} solve the same regularized
equation using a moving collocation method.
Sun et al. \cite{SunRusXu2007} solve (\ref{lub-1}) in two dimensions using an adaptive finite element method
and show that proper mesh adaptation can provide accurate resolution and there is no need 
to use regularization in the numerical simulation of the singularity development in (\ref{lub-1}).

We first consider an IBVP of the one-dimensional lubrication-type equation as
\beq
\begin{cases}
& u_t +  (f(u) u_{xxx})_x = 0, \;\; x \in \Omega \equiv (-1,1),\;\; f(u) = u^{\frac{1}{2}},
\\
& u(t,\pm 1) = 0, \; u_{xxx} (t, \pm 1) = 0,
\\
& u(0,x) = 0.8 -\cos (\pi x) + 0.25 \cos (2\pi x). 
\end{cases}
\label{ibvp-4}
\eeq
Notice that the initial condition is strictly positive (with minimum value 0.05 and mean value 0.8)
and $f(u)$ is undefined for negative $u$.
We use a uniform mesh of $J$ cells for $\Omega$, central finite
differences for spatial discretization, a third-order SDIRK method for the temporal discretization
of the underlying PDE. (A small fixed time step $\Delta t = 10^{-6}$ is used in our computation.)
This discretization is basically the same as that for the example in the previous section except
that a special treatment is needed for the diffusion coefficient for the current problem
since it is nonlinear. Recall that the scheme does not preserve the nonnegativity of the solution.
Thus, iterates can have negative values at some point during the Newton iterative solution
of the nonlinear algebraic equations resulting from the SDIRK/finite difference discretization 
of the underlying PDE.
Once this occurs, the computation stops because $f(u)$ is undefined for those values.
One way to avoid this difficulty is to use $f((U^{n+1})^{+})$ instead of $f(U^{n+1})$ in the scheme.
However, the nonsmoothness of $(U^{n+1})^{+}$ can cause difficulty in the 
convergence of the Newton iteration.
We use here the ``lagged diffusivity" method, i.e., the diffusion coefficient is computed using
the (corrected) numerical solution at the previous time step. An iterative method based on
the lagged diffusivity is used and proved to be convergent by Dobson and Vogel \cite{DoVo1997}
for a total variation denoising problem which also has a nonlinear diffusion coefficient.

%Here, the ``lagged diffusivity" method is used to discretize the time as described in \cite{DoVo1997}, i.e.
%\[
%\frac{u^{n+1}-u^n}{\Delta t} + \partial _x (f(u^n) \partial_x^3 u^{n+1}) = 0.
%\]
%The center finite difference method is used for spatial discretization on uniform mesh.

%The solution of this problem develops singularity in finite time. In \cite{ZhoBer2000}, the regularization
%\beq
%u_{\epsilon t} + \partial _x (f_{\epsilon}(u_{\epsilon}) \partial_x^3 u_{\epsilon}) = 0, \;\; f_{\epsilon}(u_{\epsilon}) = \frac{u_{\epsilon}^4 f(u_{\epsilon})}{\epsilon f(u_{\epsilon}) + u_{\epsilon}^4},
%\label{regularization}
%\eeq
%is used to guarantee nonnegative, where $\epsilon$ is a small positive constant.
%In this paper, the technique in section 2 is used for the FDM of the Lubrication-type equations
%to assure the nonnagetive numerical solution.


Fig.~\ref{Exa4.1-1} shows a numerical solution at various time instants obtained with the cutoff method
(without using regularization for $f(u)$) on a relatively coarse uniform mesh of 129 points. The result is almost identical
(in the eyeball norm) to those in \cite{Ber1995,SunRusXu2007}. To further verify the result, the computation
is done with a uniform mesh of 1001 points and the obtained solution is indistinguishable from the one shown
in Fig.~\ref{Exa4.1-1} in the plotting precision.

A uniform mesh of 1001 points is used for the numerical study of the singularity development of the solution.
Fig.~\ref{Exa4.1-2} shows the close views of the solution during the development. It can be seen that
the singularity is developed at around $t= 7.30\times 10^{-4}$, which is consistent with existing
numerical simulations,
including those in \cite{ZhoBer2000} with a positivity-preserving scheme.
Then, the solution becomes identically zero on an interval of time dependent length for a period
of time. The length of the interval increases from zero, attains its maximum (about $0.12$), and decreases to
zero. Afterward, the solution becomes positive again, and eventually decays to the mean of the initial solution.
It is known (e.g., see \cite{Ber1995}) that the 
development of the singularity is characterized by 
the third order derivative of the solution becoming infinite at certain points. 
This can be seen in Fig.~\ref{Exa4.1-3} where the third order derivative of the numerical solution is shown
at various time instants. In particular, the third order derivative becomes discontinuous when 
the singularity
occurs, then has jump discontinuities, and then becomes continuous again when the solution is positive.

As we have seen in the above, there is no need to use any regularization for $f(u)$ (cf. (\ref{lub-2})) with the cutoff
method. On the other hand, it is interesting to see the effects of the regularization of $f(u)$ on the solution.
We first point out that the regularization does not guarantee the nonnegativity of the solution for
the central finite difference discretization on a uniform mesh. Thus, we also need to use
the cutoff method for the regularized equation (\ref{lub-2}).
A solution obtained for (\ref{lub-2}) (with $\epsilon = 10^{-14}$) with the cutoff method is shown
in Fig.~\ref{Exa4.1-4}. By a direct comparison with Fig.~\ref{Exa4.1-3}, one can see that there
is only a slight difference (at the level of $\mathcal{O}(10^{-6})$) in the onset time and disappearance
time of the singularity. Moreover, the solution of the regularized problem does not become identically
zero on an interval for a period of time. Instead, there is a bump at the central part of the interval and
the height of the bump maintains constant almost for the whole appearance period of singularity.
Once again, the length of the interval first increases from zero, reaches its maximum, and decreases to zero
when the solution becomes positive again. 

Fig.~\ref{Exa4.1-6} shows comparison of numerical solutions at onset of singularity and at $t= 0.001$.
It can be seen that the regularization changes the onset pattern: the solution to the non-regularized equation
touches the $x$-axis at one point whereas those to the regularized one touch at two points simultaneously.
Moreover, the length of the touching interval and the height of the bump at onset depend on the regularization
parameter $\epsilon$. The smaller $\epsilon$, the smaller the touching interval and the lower the bump.
This suggests $u^\epsilon \to u$ as $\epsilon \to 0$. It is interesting to point out that all solutions, for regularization
or non-regularization, have almost the same length of the touching interval at $t=0.001$.

Finally, we present some numerical results for an IBVP of the two-dimensional lubrication-type equation as
\beq
\begin{cases}
& u_t + \nabla \cdot (f(u) \nabla \Delta u) = 0, \quad  (x,y) \in  \Omega \equiv [-1,1] \times [-1,1], \;\; f(u) = u^{\frac{1}{2}},
\\
& u = 0, \; \; \frac{\p \Delta u}{\p n}  = 0, \quad (x,y) \in  \partial \Omega
\\
& u(0, x,y) =( 0.8 -\cos (\pi x) + 0.25 \cos (2\pi x) ) ( 0.8 -\cos (\pi y) + 0.25 \cos (2\pi y) ), \quad  (x,y) \in  \Omega .
\end{cases}
\label{ibvp-5}
\eeq
Fig.~\ref{Exa4.2-1} (a) shows a numerical solution at $t = 0.001$ obtained with a uniform mesh of size $81 \times 81$
while its profiles along the diagonal line $x=y$ at various times are shown in Fig.~\ref{Exa4.2-1} (b).
One can see that the development of singularity has a similar behavior as in one dimension (cf. Fig.~\ref{Exa4.1-1}).
One may also see that the onset time for the current example is slightly earlier. This is because the current
initial solution has a smaller minimum value 0.025.
The contours of the numerical solutions at $t=0.001$ before and after cutoff are shown
in Fig.~\ref{Exa4.2-2}. The undershoot (with $u_{min}=-3.23\times 10^{-4}$) is visible
in the solution before cutoff. These results demonstrate that the cutoff method can be used
for the simulation of singularity development in the two-dimensional lubrication-type problem
without need of any type of regularization.

%Exa4.1-1
\begin{figure}[thb]
\centering
\includegraphics[width=4in]{ex3-soln-J129-dt1e-6-t5e-2-color.pdf}
%\hbox{
%\begin{minipage}[b]{3.0in}
%\centerline{(a) solution obtained with 129 grid points}
%\includegraphics[width=4in]{ex3-soln-J129-dt1e-6-t5e-2-color.pdf}
%\end{minipage}
%\begin{minipage}[b]{3.0in}
%\centerline{(b) solution obtained with 1001 grid points}
%\includegraphics[width=4in]{ex3-soln-J1001-dt1e-6-t5e-2-color.pdf}
%\end{minipage}
%}
\caption{The one-dimensional lubrication-type problem. The solution obtained with 129 uniform grid points is
shown at various time instants.}
\label{Exa4.1-1}
\end{figure}

%Exa4.1-2
\begin{figure}[thb]
\centering
\hbox{
\begin{minipage}[b]{2.0in}
\centerline{(a): t=7.29e-4}
\includegraphics[width=2in]{ex3-soln-J1001-dt1e-6-a.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(b): t=7.5e-4}
\includegraphics[width=2in]{ex3-soln-J1001-dt1e-6-b.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(c): t=8.0e-4}
\includegraphics[width=2in]{ex3-soln-J1001-dt1e-6-c.pdf}
\end{minipage}
}
\hbox{
\begin{minipage}[b]{2.0in}
\centerline{(d): t=8.5e-4}
\includegraphics[width=2in]{ex3-soln-J1001-dt1e-6-d.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(e): t=1.0e-3}
\includegraphics[width=2in]{ex3-soln-J1001-dt1e-6-f.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(f): t=1.5e-3}
\includegraphics[width=2in]{ex3-soln-J1001-dt1e-6-g.pdf}
\end{minipage}
}
\hbox{
\begin{minipage}[b]{2.0in}
\centerline{(g): t=2.0e-3}
\includegraphics[width=2in]{ex3-soln-J1001-dt1e-6-h.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(h): t=2.3e-3}
\includegraphics[width=2in]{ex3-soln-J1001-dt1e-6-i.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(i): t=2.34e-3}
\includegraphics[width=2in]{ex3-soln-J1001-dt1e-6-j.pdf}
\end{minipage}
}
\caption{The one-dimensional lubrication-type problem. The close views of the numerical solution
are shown at various time instants during the singularity development. The numerical solution is obtained
by the cutoff method (without using regularization for $f(u)$) on a uniform mesh of 1001 points
with $\Delta t = 10^{-6}$.}
\label{Exa4.1-2}
\end{figure}

%Exa4.1-3
\begin{figure}[thb]
\centering
\hbox{
\begin{minipage}[b]{2.0in}
\centerline{(a): t=7.29e-4}
\includegraphics[width=2in]{ex3-d-soln-J1001-dt1e-6-a.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(b): t=7.5e-4}
\includegraphics[width=2in]{ex3-d-soln-J1001-dt1e-6-b.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(c): t=8.0e-4}
\includegraphics[width=2in]{ex3-d-soln-J1001-dt1e-6-c.pdf}
\end{minipage}
}
\hbox{
\begin{minipage}[b]{2.0in}
\centerline{(d): t=8.5e-4}
\includegraphics[width=2in]{ex3-d-soln-J1001-dt1e-6-d.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(e): t=1.0e-3}
\includegraphics[width=2in]{ex3-d-soln-J1001-dt1e-6-f.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(f): t=1.5e-3}
\includegraphics[width=2in]{ex3-d-soln-J1001-dt1e-6-g.pdf}
\end{minipage}
}
\hbox{
\begin{minipage}[b]{2.0in}
\centerline{(g): t=2.0e-3}
\includegraphics[width=2in]{ex3-d-soln-J1001-dt1e-6-h.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(h): t=2.3e-3}
\includegraphics[width=2in]{ex3-d-soln-J1001-dt1e-6-i.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(i): t=2.34e-3}
\includegraphics[width=2in]{ex3-d-soln-J1001-dt1e-6-j.pdf}
\end{minipage}
}
\caption{The one-dimensional lubrication-type problem. The close views of the third order derivative
of the numerical solution are shown at various time instants during the singularity development.
The numerical solution is obtained by the cutoff method (without using regularization for $f(u)$)
on a uniform mesh of 1001 points with $\Delta t = 10^{-6}$.}
\label{Exa4.1-3}
\end{figure}


%Exa4.1-4
\begin{figure}[thb]
\centering
\hbox{
\begin{minipage}[b]{2.0in}
\centerline{(a): t=7.3e-4}
\includegraphics[width=2in]{ex3-reg-soln-J1001-dt1e-6-a.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(b): t=7.5e-4}
\includegraphics[width=2in]{ex3-reg-soln-J1001-dt1e-6-b.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(c): t=8.0e-4}
\includegraphics[width=2in]{ex3-reg-soln-J1001-dt1e-6-c.pdf}
\end{minipage}
}
\hbox{
\begin{minipage}[b]{2.0in}
\centerline{(d): t=8.5e-4}
\includegraphics[width=2in]{ex3-reg-soln-J1001-dt1e-6-d.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(e): t=1.0e-3}
\includegraphics[width=2in]{ex3-reg-soln-J1001-dt1e-6-f.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(f): t=1.5e-3}
\includegraphics[width=2in]{ex3-reg-soln-J1001-dt1e-6-g.pdf}
\end{minipage}
}
\hbox{
\begin{minipage}[b]{2.0in}
\centerline{(g): t=2.0e-3}
\includegraphics[width=2in]{ex3-reg-soln-J1001-dt1e-6-h.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(h): t=2.3e-3}
\includegraphics[width=2in]{ex3-reg-soln-J1001-dt1e-6-i.pdf}
\end{minipage}
\begin{minipage}[b]{2.0in}
\centerline{(i): t=2.34e-3}
\includegraphics[width=2in]{ex3-reg-soln-J1001-dt1e-6-j.pdf}
\end{minipage}
}
\caption{The one-dimensional regularized lubrication-type problem (with $\epsilon = 10^{-14}$).
The close views of the numerical solution
are shown at various time instants during the singularity development. The numerical solution is obtained
by the cutoff method with regularization for $f(u)$ on a uniform mesh of 1001 points with
$\Delta t = 10^{-6}$.}
\label{Exa4.1-4}
\end{figure}

%%Exa4.1-5
%\begin{figure}[thb]
%\centering
%\hbox{
%\begin{minipage}[b]{2.0in}
%\centerline{(a): t=7.29e-4}
%\includegraphics[width=2in]{ex3-reg-d-soln-J1001-dt1e-6-a.pdf}
%\end{minipage}
%\begin{minipage}[b]{2.0in}
%\centerline{(b): t=7.5e-4}
%\includegraphics[width=2in]{ex3-reg-d-soln-J1001-dt1e-6-b.pdf}
%\end{minipage}
%\begin{minipage}[b]{2.0in}
%\centerline{(c): t=8.0e-4}
%\includegraphics[width=2in]{ex3-reg-d-soln-J1001-dt1e-6-c.pdf}
%\end{minipage}
%}
%\hbox{
%\begin{minipage}[b]{2.0in}
%\centerline{(d): t=8.5e-4}
%\includegraphics[width=2in]{ex3-reg-d-soln-J1001-dt1e-6-d.pdf}
%\end{minipage}
%\begin{minipage}[b]{2.0in}
%\centerline{(e): t=1.0e-3}
%\includegraphics[width=2in]{ex3-reg-d-soln-J1001-dt1e-6-f.pdf}
%\end{minipage}
%\begin{minipage}[b]{2.0in}
%\centerline{(f): t=1.5e-3}
%\includegraphics[width=2in]{ex3-reg-d-soln-J1001-dt1e-6-g.pdf}
%\end{minipage}
%}
%\hbox{
%\begin{minipage}[b]{2.0in}
%\centerline{(g): t=2.0e-3}
%\includegraphics[width=2in]{ex3-reg-d-soln-J1001-dt1e-6-h.pdf}
%\end{minipage}
%\begin{minipage}[b]{2.0in}
%\centerline{(h): t=2.3e-3}
%\includegraphics[width=2in]{ex3-reg-d-soln-J1001-dt1e-6-i.pdf}
%\end{minipage}
%\begin{minipage}[b]{2.0in}
%\centerline{(i): t=2.34e-3}
%\includegraphics[width=2in]{ex3-reg-d-soln-J1001-dt1e-6-j.pdf}
%\end{minipage}
%}
%\caption{Example of one dimensional Lubirication problem. Solutions of the third derivation at different times near to singularity.}
%\label{Exa4.1-5}
%\end{figure}


%Exa4.1-6
\begin{figure}[thb]
\centering
\hbox{
\begin{minipage}[b]{3.0in}
\centerline{(a): Solutions at onset of singularity}
\includegraphics[width=3.0in]{ex3-soln-t0-J1001-c.pdf}
\end{minipage}
\begin{minipage}[b]{3.0in}
\centerline{(b): Solutions at $t=0.001$}
\includegraphics[width=3.0in]{ex3-soln-t1e-3-J1001-c.pdf}
\end{minipage}
}
\caption{The one-dimensional lubrication-type problem with/without regularization. Numerical solutions
at onset of singularity or at $t=0.001$.}
\label{Exa4.1-6}
\end{figure}

%Exa4.2-1
\begin{figure}[thb]
\centering
\hbox{
\begin{minipage}[b]{3.0in}
\centerline{(a): Numerical solution at $t=0.001$}
\includegraphics[height=2.5in]{ex3-2d-solution1e-3.pdf}
\end{minipage}
\begin{minipage}[b]{3.0in}
\centerline{(b): Solution along $x=y$ at various times}
\includegraphics[height=2.5in]{ex3-2d-soln-J81-dt1e-6-t1e-3-c.pdf}
\end{minipage}
}
\caption{The two-dimensional lubrication problem. The solution is obtained with the cutoff method with
a uniform mesh of size $81\times 81$.}
\label{Exa4.2-1}
\end{figure}


%Exa4.2-2
\begin{figure}[thb]
\centering
\hbox{
\begin{minipage}[b]{3.0in}
\centerline{(a): $U^N$}
\includegraphics[width=3.5in]{ex3-2D-contour81-t1e-3-before.pdf}
\end{minipage}
\begin{minipage}[b]{3.0in}
\centerline{(b): $(U^N)^+$}
\includegraphics[width=3.5in]{ex3-2D-contour81-t1e-3-after.pdf}
\end{minipage}
}
\caption{The two-dimensional lubrication problem. Contours of the numerical solutions at $t=0.001$ before and after cutoff.}
\label{Exa4.2-2}
\end{figure}



% section 5
\section{Conclusions and further comments}
\label{SEC:conclusion}

In the previous section we have studied the cutoff method for the numerical computation of
nonnegative solutions of parabolic partial differential equations. Several properties
of the cutoff method are given in Lemmas~\ref{lem2.1} -- \ref{lem2.3}. Convergence of a class of
finite difference methods is proved (Theorems \ref{thm2.1} and \ref{thm2.2}) when they are incorporated
with the cutoff method for linear parabolic PDEs. The method is investigated for two applications, linear anisotropic
diffusion problems and nonlinear lubrication-type PDEs. The numerical results are consistent with theoretical
predictions and in good agreement with existing results in the literature.

We have considered finite difference methods in this work.
But it is worth pointing out that the cutoff method can also be used with other discretization methods, e.g.
collocation, finite volume, finite element, or spectral methods.
As an example, we show in Figs.~\ref{Exa3.1-6} and \ref{Exa3.1-7} results obtained with the standard linear
finite element method on Delaunay meshes for anisotropic diffusion problem (\ref{ibvp-2}).
They are comparable with those in Figs.~\ref{Exa3.1-2} and \ref{Exa3.1-3}.
Theoretical analysis for finite element methods with the cutoff strategy is currently underway.
%(with $u_{min} = -4.4017\time 10^{-4}$)


%%Exa3.1-5
%\begin{figure}[thb]
%\centering
%\includegraphics[width=3.5in]{ex1-fem-mesh2360.pdf}
%\caption{Example~\ref{exam3.1}. A Delaunay mesh of $N=2388$ cells.}
%\label{Exa3.1-5}
%\end{figure}

%Exa3.1-6
\begin{figure}[thb]
\centering
\hbox{
\begin{minipage}[b]{3in}
\centerline{(a): $U^N$}
\includegraphics[width=3.5in]{ex1-fem-contour-N2350-before.pdf}
\end{minipage}
\begin{minipage}[b]{3in}
\centerline{(b): $(U^N)^+$}
\includegraphics[width=3.5in]{ex1-fem-contour-N2350-after.pdf}
\end{minipage}
}
\caption{The standard linear finite element method for Example~(\ref{ibvp-2}).
Contours of the numerical solutions at $t = 1$ before and after cutoff obtained 
with a Delaunay mesh of $N=2350$ elements.}
\label{Exa3.1-6}
\end{figure}

%Exa3.1-7
\begin{figure}[thb]
\centering
\includegraphics[scale=0.8]{ex1-fem-error-umin-color.pdf}
\caption{The standard linear finite element method for Example~(\ref{ibvp-2}).
The $L^2$ error, $\| (U^N)^{+}-u^N \|_{L^2(\Omega)}$, and the maximal undershoot, $-u_{min}$,
are shown as functions of $\sqrt{N}$ with $N$ being the number of elements.
}
\label{Exa3.1-7}
\end{figure}

\vspace{20pt}

{\bf Acknowledgment.} 
The work was supported in part by the National
Science Foundation (U.S.A.) under Grants DMS-1115118 and DMS-1115408
and the Natural Science Foundation of China under Grant 40906048.


%---------------------------------------------------
%\bibliographystyle{plain}
%\bibliographystyle{abbrv}
%\bibliography{/Users/huang/tex/bib/mmesh}
%\bibliography{/Users/changnalu/tex/bib/paper}

\begin{thebibliography}{10}

\bibitem{Ale77}
R.~Alexander.
\newblock Diagonally implicit {Runge-Kutta} methods for stiff {O.D.E.}'s.
\newblock {\em SIAM J. Numer. Anal.}, 14:1006 -- 1021, 1977.

\bibitem{BeFr1990}
F.~Bernis and A.~Friedman.
\newblock Higher order nonlinear degenerate parabolic equations.
\newblock {\em J. Differential Equations}, 83:179--206, 1990.

\bibitem{Ber1995}
A.~L. Bertozzi.
\newblock Loss and gain of regularity in a lubrication equation for thin
  viscous films.
\newblock In {\em Free boundary problems: theory and applications ({T}oledo,
  1993)}, volume 323 of {\em Pitman Res. Notes Math. Ser.}, pages 72--85.
  Longman Sci. Tech., Harlow, 1995.

\bibitem{Ber1996}
A.~L. Bertozzi.
\newblock Symmetric singularity formation in lubrication-type equations for
  interface motion.
\newblock {\em SIAM J. Appl. Math.}, 56:681--714, 1996.

\bibitem{BBDK1994}
A.~L. Bertozzi, M.~P. Brenner, T.~F. Dupont, and L.~P. Kadanoff.
\newblock Singularities and similarities in interface flows.
\newblock In {\em Trends and perspectives in applied mathematics}, volume 100
  of {\em Appl. Math. Sci.}, pages 155--208. Springer, New York, 1994.

\bibitem{BKK08}
J.~Brandts, S.~Korotov, and M.~K\v{r}\'i\v{z}ek.
\newblock The discrete maximum principle for linear simplicial finite element
  approximations of a reaction-diffusion problem.
\newblock {\em Lin. Alg. Appl.}, 429:2344--2357, 2008.

\bibitem{BE04}
E.~Burman and A.~Ern.
\newblock Discrete maximum principle for {G}alerkin approximations of the
  {L}aplace operator on arbitrary meshes.
\newblock {\em C. R. Acad. Sci. Paris}, Ser.I 338:641--646, 2004.

\bibitem{Cas79}
J.~R. Cash.
\newblock Diagonally implicit runge-kutta formulae with error estimates.
\newblock {\em J. Inst. Math. Appl.}, 24:293--301, 1979.

\bibitem{Chorin1968}
A.~J. Chorin.
\newblock Numerical solution of the navier-stokes equations.
\newblock {\em Math. Comp.}, 22:745 -- 762, 1968.

\bibitem{Cia70}
P.~G. Ciarlet.
\newblock Discrete maximum principle for finite difference operators.
\newblock {\em Aequationes Math.}, 4:338--352, 1970.

\bibitem{CR73}
P.~G. Ciarlet and P.-A. Raviart.
\newblock Maximum principle and uniform convergence for thefinite element
  method.
\newblock {\em Comput. Meth. Appl. Mech. Engrg.}, 2:17--31, 1973.

\bibitem{CDGKS1993}
P.~Constantin, T.~F. Dupont, R.~E. Goldstein, L.~P. Kadanoff, M.~J. Shelley,
  and S.-M. Zhou.
\newblock Droplet breakup in a model of the {H}ele-{S}haw cell.
\newblock {\em Phys. Rev. E (3)}, 47:4169--4181, 1993.

\bibitem{CE92}
M.~I.~M. Copetti and C.~M. Elliott.
\newblock Numerical analysis of the {C}ahn-{H}illiard equation with a
  logarithmic free energy.
\newblock {\em Numer. Math.}, 63:39--65, 1992.

\bibitem{DoVo1997}
D.~C. Dobson and C.~R. Vogel.
\newblock Convergence of an iterative method for total variation denoising.
\newblock {\em SIAM J. Numer. Anal.}, 34:1779--1791, 1997.

\bibitem{DDS04}
A.~Dr\v{a}g\v{a}nescu, T.~F. Dupont, and L.~R. Scott.
\newblock Failure of the discrete maximum principle for an elliptic finite
  element problem.
\newblock {\em Math. Comp.}, 74:1--23, 2004.

\bibitem{EllGar1996}
C.~M. Elliott and H.~Garcke.
\newblock On the {C}ahn-{H}illiard equation with degenerate mobility.
\newblock {\em SIAM J. Math. Anal.}, 27:404--423, 1996.

\bibitem{FH2001}
I.~Farag{\'o} and R.~Horv{\'a}th.
\newblock On the nonnegativity conservation of finite element solutions of
  parabolic problems.
\newblock In {\em Finite element methods ({J}yv\"askyl\"a, 2000)}, volume~15 of
  {\em GAKUTO Internat. Ser. Math. Sci. Appl.}, pages 76--84. Gakk\=otosho,
  Tokyo, 2001.

\bibitem{FH06}
I.~Farag{\'o} and R.~Horv{\'a}th.
\newblock Discrete maximum principle and adequate discretizations of linear
  parabolic problems.
\newblock {\em SIAM J. Sci. Comput.}, 28:2313--2336, 2006.

\bibitem{FHK05}
I.~Farag{\'o}, R.~Horv{\'a}th, and S.~Korotov.
\newblock Discrete maximum principle for linear parabolic problems solved on
  hybrid meshes.
\newblock {\em Appl. Numer. Math.}, 53:249--264, 2005.

\bibitem{FHK2011}
I.~Farag{\'o}, R.~Horv{\'a}th, and S.~Korotov.
\newblock Discrete maximum principles for {FE} solutions of nonstationary
  diffusion-reaction problems with mixed boundary conditions.
\newblock {\em Numer. Methods Partial Differential Equations}, 27:702--720,
  2011.

\bibitem{Fuj73}
H.~Fujii.
\newblock Some remarks on finite element analysis of time-dependent field
  problems.
\newblock In {\em Theory and Proactice in Finite Element Structural Analysis},
  pages 91--106. University of Tokyo, Tokyo, 1973.

\bibitem{GL09}
S.~G\H{u}nter and K.~Lackner.
\newblock A mixed implicit-explicit finite difference scheme for heat transport
  in magnetised plasmas.
\newblock {\em J. Comput. Phys.}, 228:282--293, 2009.

\bibitem{GYK05}
S.~G\H{u}nter, Q.~Yu, J.~Kruger, and K.~Lackner.
\newblock Modelling of heat transport in magnetised plasmas using non-aligned
  coordinates.
\newblock {\em J. Comput. Phys.}, 209:354--370, 2005.

\bibitem{Gre1978}
H.~P. Greenspan.
\newblock On the motion of a small viscous droplet that wets a surface.
\newblock {\em J. Fluid Mech.}, 84:125--143, 1978.

\bibitem{Grun1995}
G.~Gr{\"u}n.
\newblock Degenerate parabolic differential equations of fourth order and a
  plasticity model with non-local hardening.
\newblock {\em Z. Anal. Anwendungen}, 14:541--574, 1995.

\bibitem{Hai2006}
E.~Hairer.
\newblock Long-time energy conservation of numerical integrators.
\newblock In {\em Foundations of computational mathematics, {S}antander 2005},
  volume 331 of {\em London Math. Soc. Lecture Note Ser.}, pages 162--180.
  Cambridge Univ. Press, Cambridge, 2006.

\bibitem{Hua10}
W.~Huang.
\newblock Discrete maximum principle and a delaunay-type mesh condition for
  linear finite element approximations of two-dimensional anisotropic diffusion
  problems.
\newblock {\em Numer. Math. Theory Meth. Appl.}, 4:319--334, 2011.
\newblock (arXiv:1008.0562).

\bibitem{KK09}
J.~Kar{\'a}tson and S.~Korotov.
\newblock An algebraic discrete maximum principle in {H}ilbert space with
  applications to nonlinear cooperative elliptic systems.
\newblock {\em SIAM J. Numer. Anal.}, 47:2518--2549, 2009.

\bibitem{KKK07}
J.~Kar\'atson, S.~Korotov, and M.~K\v{r}\'i\v{z}ek.
\newblock On discrete maximum principles for nonlinear elliptic problems.
\newblock {\em Math. Comput. Sim.}, 76:99--108, 2007.

\bibitem{KSS09}
D.~Kuzmin, M.~J. Shashkov, and D.~Svyatskiy.
\newblock A constrained finite element method satisfying the discrete maximum
  principle for anisotropic diffusion problems.
\newblock {\em J. Comput. Phys.}, 228:3448--3463, 2009.

\bibitem{KL95}
M.~K\v{r}\'i\v{z}ek and Q.~Lin.
\newblock On diagonal dominance of stiffness matrices in {3D}.
\newblock {\em East-West J. Numer. Math.}, 3:59--69, 1995.

\bibitem{LePot05}
C.~Le~Potier.
\newblock Sch\'ema volumes finis monotone pour des op\'erateurs de diffusion
  fortement anisotropes sur des maillages de triangles non structur\'es.
\newblock {\em C. R. Math. Acad. Sci. Paris}, 341:787--792, 2005.

\bibitem{LePot09b}
C.~Le~Potier.
\newblock Un sch\'ema lin\'eaire v\'erifiant le principe du maximum pour des
  op\'erateurs de diffusion tr\`es anisotropes sur des maillages d\'eform\'es.
\newblock {\em C. R. Math. Acad. Sci. Paris}, 347:105--110, 2009.

\bibitem{Let92}
F.~W. Letniowski.
\newblock Three-dimensional {D}elaunay triangulations for finite element
  approximations to a second-order diffusion operator.
\newblock {\em SIAM J. Sci. Stat. Comput.}, 13:765--770, 1992.

\bibitem{Lew1994}
M.~A. Lewis.
\newblock Spatial coupling of plant and herbivore dynamics: The contribution of
  herbivore dispersal to transient and persistent ``waves" of damage.
\newblock {\em Theoret. Population Biol.}, 45:277--312, 1994.

\bibitem{LH10}
X.~P. Li and W.~Huang.
\newblock An anisotropic mesh adaptation method for the finite element solution
  of heterogeneous anisotropic diffusion problems.
\newblock {\em J. Comput. Phys.}, 229:8072--8094, 2010 (arXiv:1003.4530v2).

\bibitem{LiHu2012}
X.~P. Li and W.~Huang.
\newblock Maximum principle for the finite element solution of time dependent
  anisotropic diffusion problems.
\newblock {\em Numer. Meth. P.D.E.}, 2012.
\newblock (arXiv:1209.5657).

\bibitem{LSS07}
X.~P. Li, D.~Svyatskiy, and M.~Shashkov.
\newblock Mesh adaptation and discrete maximum principle for {2D} anisotropic
  diffusion problems.
\newblock Technical Report LA-UR 10-01227, Los Alamos National Laboratory, Los
  Alamos, NM, 2007.

\bibitem{LSSV07}
K.~Lipnikov, M.~Shashkov, D.~Svyatskiy, and Y.~Vassilevski.
\newblock Monotone finite volume schemes for diffusion equations on
  unstructured triangular and shape-regular polygonal meshes.
\newblock {\em J. Comput. Phys.}, 227:492--512, 2007.

\bibitem{LS08}
R.~Liska and M.~Shashkov.
\newblock Enforcing the discrete maximum principle for linear finite element
  solutions of second-order elliptic problems.
\newblock {\em Comm. Comput. Phys.}, 3:852--877, 2008.

\bibitem{LuHuQi2012}
C.~Lu, W.~Huang, and J.~Qiu.
\newblock Maximum principle in linear finite element approximations of
  anisotropic diffusion-convection-reaction problems.
\newblock {\em Numer. Math.} in press, 2012 (arXiv:1201.3564v1).

\bibitem{MD06}
M.~J. Mlacnik and L.~J. Durlofsky.
\newblock Unstructured grid optimization for improved monotonicity of discrete
  solutions of elliptic equations with highly anisotropic coefficients.
\newblock {\em J. Comput. Phys.}, 216:337--361, 2006.

\bibitem{MM05}
K.~W. Morton and D.~F. Mayers.
\newblock {\em Numerical Solution of Partial Differential Equations}.
\newblock Cambridge University Press, Cambridge, U.K., 2005.
\newblock 2nd Edition.

\bibitem{RuWiXu2007}
R.~D. Russell, J.~F. Williams, and X.~Xu.
\newblock M{OVCOL}4: a moving mesh code for fourth-order time-dependent partial
  differential equations.
\newblock {\em SIAM J. Sci. Comput.}, 29:197--220 (electronic), 2007.

\bibitem{SH07}
P.~Sharma and G.~W. Hammett.
\newblock Preserving monotonicity in anisotropic diffusion.
\newblock {\em J. Comput. Phys.}, 227:123--142, 2007.

\bibitem{Sto86}
G.~Stoyan.
\newblock On maximum principles for monotone matrices.
\newblock {\em Lin. Alg. Appl.}, 78:147--161, 1986.

\bibitem{SF73}
G.~Strang and G.~J. Fix.
\newblock {\em An Analysis of the Finite Element Method}.
\newblock Prentice Hall, Englewood Cliffs, NJ, 1973.

\bibitem{SunRusXu2007}
P.~Sun, R.~D. Russell, and J.~Xu.
\newblock A new adaptive local mesh refinement algorithm and its application on
  fourth order thin film flow problem.
\newblock {\em J. Comput. Phys.}, 224:1021--1048, 2007.

\bibitem{Teman1968}
R.~Temam.
\newblock Une m\'ethode d'approximation de la solution des \'equations de
  {N}avier-{S}tokes.
\newblock {\em Bull. Soc. Math. France}, 96:115--152, 1968.

\bibitem{TW08}
V.~Thom{\'e}e and L.~B. Wahlbin.
\newblock On the existence of maximum principles in parabolic finite element
  equations.
\newblock {\em Math. Comput.}, 77:11--19, 2008.

\bibitem{WaZh11}
J.~Wang and R.~Zhang.
\newblock Maximum principle for {P1}-conforming finite element approximations
  of quasi-linear second order elliptic equations.
\newblock {\em SIAM J Numer. Anal.}, 50:626--642, 2012.
\newblock (arXiv:1105.1466).

\bibitem{XZ99}
J.~Xu and L.~Zikatanov.
\newblock A monotone finite element scheme for convection-diffusion equations.
\newblock {\em Math. Comput.}, 69:1429--1446, 1999.

\bibitem{ZhSh2010}
X.~Zhang and C.-W. Shu.
\newblock On maximum-principle-satisfying high order schemes for scalar
  conservation laws.
\newblock {\em J. Comput. Phys.}, 229:3091--3120, 2010.

\bibitem{ZhoBer2000}
L.~Zhornitskaya and A.~L. Bertozzi.
\newblock Positivity-preserving numerical schemes for lubrication-type
  equations.
\newblock {\em SIAM J. Numer. Anal.}, 37:523--555, 2000.

\end{thebibliography}

\end{document}

\begin{figure}[thb]
\centering
\hbox{
\begin{minipage}[t]{1.8in}
\includegraphics[width=2in, height=2in]{ade1e/mesh11b.pdf}
\end{minipage}
\hspace{10mm}
\begin{minipage}[t]{1.8in}
\includegraphics[width=2in, height=2in]{ade1e/mesh16b.pdf}
\end{minipage}
\hspace{10mm}
\begin{minipage}[t]{1.8in}
\includegraphics[width=2in, height=2in]{ade1e/mesh20b.pdf}
\end{minipage}
}
\hbox{
\begin{minipage}[t]{1.8in}
\includegraphics[width=2in, height=2in]{ade1e/mesh11c.pdf}
\end{minipage}
\hspace{10mm}
\begin{minipage}[t]{1.8in}
\includegraphics[width=2in, height=2in]{ade1e/mesh16c.pdf}
\end{minipage}
\hspace{10mm}
\begin{minipage}[t]{1.8in}
\includegraphics[width=2in, height=2in]{ade1e/mesh20c.pdf}
\end{minipage}
}
\caption{Example 3.1. Adaptive moving mesh (first row) and computed solution (second row) for ADRE.}
\label{ade1e}
\end{figure}

% table 0 
\begin{table}[htb]
\begin{center}
\caption{Quality measures and interpolation error for Shishkin-type and uniform meshes for Example (\ref{4.1}).
In the table, $|e|_{H^1(\Omega)} = (\sum_K |v - \Pi_{k,K} v|_{H^1(K)}^2 )^{\frac{1}{2}}$.
It can be shown that $Q_{soln,2} = \mbox{O}(\epsilon^{-3/2})$.}
\begin{tabular}{c|l|l|l|l|l|l}\hline \hline
 & \multicolumn{3}{|c}{Uniform Mesh} & \multicolumn{3}{|c}{Shishkin-type Mesh} \\ \hline\hline
 & $\Omega_0$ &  $\Omega_1$, $\Omega_2$ & $\Omega_3$
 & $\Omega_0$ &  $\Omega_1$, $\Omega_2$ & $\Omega_3$ \\ \cline{2-7}
 $Q_{geo}(x,y)$ & O(1) & O(1) & O(1) & O(1) & O$(\epsilon^{-1}|\ln\epsilon|^{-1})$ & O(1)\\
 $Q_{ali,2}(x,y)$ & O(1) & O($\epsilon^{-1}$) & O(1) & O(1) & O($|\ln\epsilon|$) & O(1) \\
 $Q_{adp,2}(x,y)$ & O(1) & O($\epsilon^{-1}$) & O($\epsilon^{-2}$) & O(1) & O($|\ln\epsilon|$) &
        O($|\ln\epsilon|^2$) \\ \hline\hline
$ Q_{mesh,2}$ & \multicolumn{3}{|c}{O($\epsilon^{-\frac 3 2}|\ln\epsilon|$)} &
       \multicolumn{3}{|c}{O($\epsilon^{-\frac 1 2}|\ln\epsilon|^{\frac 3 2}$)} \\ \hline
$Q_{mesh,2}/Q_{soln,2}$ & \multicolumn{3}{|c}{O(1)} & \multicolumn{3}{|c}{O($\epsilon$)} \\ \hline
$|e|_{H^1(\Omega)}$ & \multicolumn{3}{|c}{O$\left (M^{-1} \epsilon^{-\frac 3 2}|\ln \epsilon|\right )$} &
\multicolumn{3}{|c}{O$\left (M^{-1} \epsilon^{-\frac 1 2}|\ln\epsilon|^{\frac 3 2}\right )$} \\\hline \hline
\end{tabular}
\label{t.0}
\end{center}
\end{table}


% figure f-1
\begin{figure}
\begin{center}
{\footnotesize
\tikzset{my node/.code=\ifpgfmatrix\else\tikzset{matrix of nodes}\fi}
\begin{tikzpicture}[every node/.style={my node},scale=0.4] 
\draw[->,thick] (-2,0)--(0,0);
\draw[thick] (2,0) circle (2);
\node (node1) at (2,0) {$U^n$\\}; 

\draw[->,thick] (4,0)--(6,0);
\draw[thick] (6,-2) rectangle (12,2);
\node (node2) at (9,0) {cutoff\\ (removal of\\ negative values)\\};

\draw[->,thick] (12,0)--(14,0);
\draw[thick] (16,0) circle (2);
\node (node1) at (16,0) {$(U^{n})^{+}$\\}; 

\draw[->,thick] (18,0)--(20,0);
\draw[thick] (20,-2) rectangle (26,2);
\node (node2) at (23,0) {Integrate PDE\\};

\draw[->,thick] (26,0)--(28,0);
\draw[thick] (30,0) circle (2);
\node (node1) at (30,0) {$U^{n+1}$\\}; 
\draw[->,thick] (32,0)--(34,0);

\end{tikzpicture} 
}
\caption{A flow chart for the numerical solution of parabolic PDEs with cutoff.}
\label{f-1}
\end{center}
\end{figure}

We consider an initial-boundary value problem (IBVP) for a time interval $[t_0, T]$ and assume that
a numerical scheme of order $p \ge 1$ in time and order $p\ge 1$ in space is given
for its numerical solution. We make two assumptions on the IBVP.
\begin{itemize}
\item[H1.]  {\em Well posedness.} We assume that the IBVP is well posed. Specifically,
	 the IVBP has a unique solution for any given admissible initial solution $\phi$ at any given
	 time $s\in [t_0, T]$ (i.e., $u(s) = \phi$). Moreover, all the partial derivatives of the exact solution
	 of order up to $(p+1)$ in time and $(q+1)$ in space depend continuously on the initial solution $\phi$.
\item[H2.] {\em Positivity.} The exact solution of the IBVP starting from any nonnegative initial solution
	$\phi$ at any given time $s\in [t_0, T]$ will stay nonnegative for all time up to $T$.
\end{itemize}
Assumption H2 is a natural one since we are concerned with problems with nonnegative solutions.
Assumption H1 is essentially the standard well-posedness assumption except that we requires
that the continuous dependence hold also for (possible) high-order derivatives. This is because
the $(p,q)$-order numerical scheme will be used for solving the IBVP and the partial derivatives
of order up to $(p+1)$ in time and $(q+1)$ in space need to be bounded by the counterparts
of the unperturbed problem.


The time integration with cutoff is illustrated in Fig.~\ref{f-2}, where $U^n$ and $\tilde{U}^{n}$
denote the numerical solution and the corrected numerical solution at $t=t_n$ (the cutoff method
will be defined later), respectively and $u=u(t)$ (resp. $u=\widetilde{u}(t)$) is the exact solution of the IBVP starting
from the initial solution $u^0$ at $t=t_0$ (resp. $\tilde{U}^{n}$ at $t = t_n$).
We simply define the corrected solution $\widetilde{U}^{n+1}$ as the nonnegative part of $U^{n+1}$.
Specifically, denote by $U^{n+1}_j$
a nodal value of $U^{n+1}$ associated with the $j^{\text{th}}$ node. Then,
\beq
\widetilde{U}^{n+1}_j = \begin{cases} U^{n+1}_j, & \mbox{ for } U^{n+1}_j \ge 0 \\
0, & \mbox{ for } U^{n+1}_j < 0.
\end{cases}
\label{U-1}
\eeq
Bounds on the cutoff error are given in the next two lemmas in terms of the local truncation error.

\begin{lem}
\label{lem2.1}
Assume that the IBVP satisfies Assumption H2 and a numerical scheme is
employed for its numerical solution. Then,  the correction error is bounded by
\beq
\| \widetilde{U}^{n+1} - U^{n+1} \| \le \|\widetilde{\tau}^n\| ,
\label{err-4-2}
\eeq
where $\| \cdot \|$ is a proper norm, $\widetilde{\tau}^n$ is the local truncation
error defined in terms of the exact solution $\widetilde{u}$
of the IBVP with initial solution $\widetilde{U}^n$ (cf. (\ref{lte-1}) in Section~\ref{SEC:2.2}).
\end{lem}

{\bf Proof.} We first prove
\beq
\| \widetilde{U}^{n+1} - U^{n+1} \| \le \| U^{n+1} - \widetilde{u}^{n+1} \| .
\label{err-4}
\eeq
Indeed, the positivity assumption (H2) implies that  $\widetilde{u}^{n+1}$ is nonnegative, i.e.,
$\widetilde{u}^{n+1}_j \ge 0$ for all $j$. When $U^{n+1}_j \ge 0$, from (\ref{U-1}) we have
\[
| \widetilde{U}^{n+1}_j - U^{n+1}_j| = 0 \le |U^{n+1}_j - \widetilde{u}^{n+1}_j | .
\]
When $U^{n+1}_j \ge 0$, on the other hand, noticing $\widetilde{u}^{n+1}_j > 0$ we have
\[
| \widetilde{U}^{n+1}_j - U^{n+1}_j| = | U^{n+1} | \le |U^{n+1}_j - \widetilde{u}^{n+1}_j | .
\]
The above two inequalities lead to (\ref{err-4}).

Notice that the local truncation error is equal to the solution error incurred over one step.
Since both $U^{n+1}$ and $\widetilde{u}^{n+1}$ result from the same initial solution $\widetilde{U}^n$ at $t=t_n$,
we have $\widetilde{\tau}^n = U^{n+1} - \widetilde{u}^{n+1}$. Combining this with (\ref{err-4}) gives rise to (\ref{err-4-2}).
\proofend

\vspace{10pt}

It is emphasized that  the local truncation error in the above lemma depends on $\widetilde{u}$ which is inconvenient
to use. We would like to use the exact solution $u$ of the original IBVP (with initial solution $u^0$ at $t=t_0$).
This requires the continuous dependence of the exact solution on the initial solution and $\widetilde{U}^{n}$ is
close enough to $u^n$.

\begin{lem}
\label{lem2.2}
Assume that the IBVP satisfies Assumption H1 and a $(p,q)$-order numerical scheme is employed for its
numerical solution. If
\beq
\widetilde{U}^n - u^n = \mathcal{O}(\Delta t^p + h^q), \quad \mbox{ as } \Delta t (h) \to 0
\label{H3}
\eeq
then we have
\beq
\| \widetilde{\tau}^{n} \| \le C_{lte}(u)\; \Delta t\; (\Delta t^p + h^q) ,
\label{err-1}
\eeq
where $\Delta t$ and $h$ are the maximal time step size
and the maximal element size, respectively, $u$ is the exact solution of the IBVP with initial
solution $u^0$ at $t=t_0$, and $C_{lte}(u)$ is a constant depending
on partial derivatives of $u$ of order up to $(p+1)$ in time and $(q+1)$ in space.
\end{lem}

{\bf Proof.}
Since the numerical scheme is of order $p$ in time and $q$ in space,
there exists a constant $C_{lte}(\widetilde{u})$ depending
on partial derivatives of $\widetilde{u}$ of order up to $(p+1)$ in time and $(q+1)$ in space
such that the local truncation error is bounded by
\beq
\| \widetilde{\tau}^n \| \le \frac{1}{2} C_{lte}(\widetilde{u}) \; \Delta t \;(\Delta t^p + h^q).
\label{err-2}
\eeq
From (\ref{H3}) and the well-posedness assumption H1, the partial derivatives of $\widetilde{u}$
of order up to $(p+1)$ in time and $(q+1)$ in space are close to those of $u$ and therefore 
$C_{lte}(\widetilde{u}) \approx C_{lte}(u)$ (we use $ C_{lte}(\widetilde{u}) \le 2\; C_{lte}(u)$).
As a consequence, we can rewrite (\ref{err-2}) into (\ref{err-1}).
\proofend

{\bf Remarks.}
\begin{enumerate}
\item The error estimate (\ref{erreq-3}) is simply the standard error estimate plus the sum of the projection
	errors at different time levels.
\item From (\ref{erreq-3}) we can see that the projection error should be of order
	\beq
	\| \widetilde{U}^n - U^n \| = \mathcal{O} (\Delta t (\Delta t^p + h^q))
	\label{erreq-4}
	\eeq
	to maintain the same convergence order as the scheme without projection.
\item Sometimes we may define $R_h$ using interpolation (such as in mesh rezoning) and the projection error
	is typically of order
	\beq
	\| \widetilde{U}^n - U^n \| = \mathcal{O}(h^{q+1}) .
	\label{erreq-5}
	\eeq
	In this case, the scheme's convergence can deteriorate or even becomes divergent as $\Delta t \to 0$
	while $h$ is fixed. 
\end{enumerate}


The analysis can be generalized to more general situations including
\begin{itemize}
\item Neumann boundary conditions; Robin boundary conditions.
\item Nonlinear PDEs. The key is to derive a stability inequality.
\item Projection: interpolation (such as the mesh rezoning approach), inexact Newton iterations, project to a manifold
	(such as energy conservation)
\end{itemize}




1206.6312.extracted_bib
\begin{thebibliography}{10}

\bibitem{Ale77}
R.~Alexander.
\newblock Diagonally implicit {Runge-Kutta} methods for stiff {O.D.E.}'s.
\newblock {\em SIAM J. Numer. Anal.}, 14:1006 -- 1021, 1977.

\bibitem{BeFr1990}
F.~Bernis and A.~Friedman.
\newblock Higher order nonlinear degenerate parabolic equations.
\newblock {\em J. Differential Equations}, 83:179--206, 1990.

\bibitem{Ber1995}
A.~L. Bertozzi.
\newblock Loss and gain of regularity in a lubrication equation for thin
  viscous films.
\newblock In {\em Free boundary problems: theory and applications ({T}oledo,
  1993)}, volume 323 of {\em Pitman Res. Notes Math. Ser.}, pages 72--85.
  Longman Sci. Tech., Harlow, 1995.

\bibitem{Ber1996}
A.~L. Bertozzi.
\newblock Symmetric singularity formation in lubrication-type equations for
  interface motion.
\newblock {\em SIAM J. Appl. Math.}, 56:681--714, 1996.

\bibitem{BBDK1994}
A.~L. Bertozzi, M.~P. Brenner, T.~F. Dupont, and L.~P. Kadanoff.
\newblock Singularities and similarities in interface flows.
\newblock In {\em Trends and perspectives in applied mathematics}, volume 100
  of {\em Appl. Math. Sci.}, pages 155--208. Springer, New York, 1994.

\bibitem{BKK08}
J.~Brandts, S.~Korotov, and M.~K\v{r}\'i\v{z}ek.
\newblock The discrete maximum principle for linear simplicial finite element
  approximations of a reaction-diffusion problem.
\newblock {\em Lin. Alg. Appl.}, 429:2344--2357, 2008.

\bibitem{BE04}
E.~Burman and A.~Ern.
\newblock Discrete maximum principle for {G}alerkin approximations of the
  {L}aplace operator on arbitrary meshes.
\newblock {\em C. R. Acad. Sci. Paris}, Ser.I 338:641--646, 2004.

\bibitem{Cas79}
J.~R. Cash.
\newblock Diagonally implicit runge-kutta formulae with error estimates.
\newblock {\em J. Inst. Math. Appl.}, 24:293--301, 1979.

\bibitem{Chorin1968}
A.~J. Chorin.
\newblock Numerical solution of the navier-stokes equations.
\newblock {\em Math. Comp.}, 22:745 -- 762, 1968.

\bibitem{Cia70}
P.~G. Ciarlet.
\newblock Discrete maximum principle for finite difference operators.
\newblock {\em Aequationes Math.}, 4:338--352, 1970.

\bibitem{CR73}
P.~G. Ciarlet and P.-A. Raviart.
\newblock Maximum principle and uniform convergence for thefinite element
  method.
\newblock {\em Comput. Meth. Appl. Mech. Engrg.}, 2:17--31, 1973.

\bibitem{CDGKS1993}
P.~Constantin, T.~F. Dupont, R.~E. Goldstein, L.~P. Kadanoff, M.~J. Shelley,
  and S.-M. Zhou.
\newblock Droplet breakup in a model of the {H}ele-{S}haw cell.
\newblock {\em Phys. Rev. E (3)}, 47:4169--4181, 1993.

\bibitem{CE92}
M.~I.~M. Copetti and C.~M. Elliott.
\newblock Numerical analysis of the {C}ahn-{H}illiard equation with a
  logarithmic free energy.
\newblock {\em Numer. Math.}, 63:39--65, 1992.

\bibitem{DoVo1997}
D.~C. Dobson and C.~R. Vogel.
\newblock Convergence of an iterative method for total variation denoising.
\newblock {\em SIAM J. Numer. Anal.}, 34:1779--1791, 1997.

\bibitem{DDS04}
A.~Dr\v{a}g\v{a}nescu, T.~F. Dupont, and L.~R. Scott.
\newblock Failure of the discrete maximum principle for an elliptic finite
  element problem.
\newblock {\em Math. Comp.}, 74:1--23, 2004.

\bibitem{EllGar1996}
C.~M. Elliott and H.~Garcke.
\newblock On the {C}ahn-{H}illiard equation with degenerate mobility.
\newblock {\em SIAM J. Math. Anal.}, 27:404--423, 1996.

\bibitem{FH2001}
I.~Farag{\'o} and R.~Horv{\'a}th.
\newblock On the nonnegativity conservation of finite element solutions of
  parabolic problems.
\newblock In {\em Finite element methods ({J}yv\"askyl\"a, 2000)}, volume~15 of
  {\em GAKUTO Internat. Ser. Math. Sci. Appl.}, pages 76--84. Gakk\=otosho,
  Tokyo, 2001.

\bibitem{FH06}
I.~Farag{\'o} and R.~Horv{\'a}th.
\newblock Discrete maximum principle and adequate discretizations of linear
  parabolic problems.
\newblock {\em SIAM J. Sci. Comput.}, 28:2313--2336, 2006.

\bibitem{FHK05}
I.~Farag{\'o}, R.~Horv{\'a}th, and S.~Korotov.
\newblock Discrete maximum principle for linear parabolic problems solved on
  hybrid meshes.
\newblock {\em Appl. Numer. Math.}, 53:249--264, 2005.

\bibitem{FHK2011}
I.~Farag{\'o}, R.~Horv{\'a}th, and S.~Korotov.
\newblock Discrete maximum principles for {FE} solutions of nonstationary
  diffusion-reaction problems with mixed boundary conditions.
\newblock {\em Numer. Methods Partial Differential Equations}, 27:702--720,
  2011.

\bibitem{Fuj73}
H.~Fujii.
\newblock Some remarks on finite element analysis of time-dependent field
  problems.
\newblock In {\em Theory and Proactice in Finite Element Structural Analysis},
  pages 91--106. University of Tokyo, Tokyo, 1973.

\bibitem{GL09}
S.~G\H{u}nter and K.~Lackner.
\newblock A mixed implicit-explicit finite difference scheme for heat transport
  in magnetised plasmas.
\newblock {\em J. Comput. Phys.}, 228:282--293, 2009.

\bibitem{GYK05}
S.~G\H{u}nter, Q.~Yu, J.~Kruger, and K.~Lackner.
\newblock Modelling of heat transport in magnetised plasmas using non-aligned
  coordinates.
\newblock {\em J. Comput. Phys.}, 209:354--370, 2005.

\bibitem{Gre1978}
H.~P. Greenspan.
\newblock On the motion of a small viscous droplet that wets a surface.
\newblock {\em J. Fluid Mech.}, 84:125--143, 1978.

\bibitem{Grun1995}
G.~Gr{\"u}n.
\newblock Degenerate parabolic differential equations of fourth order and a
  plasticity model with non-local hardening.
\newblock {\em Z. Anal. Anwendungen}, 14:541--574, 1995.

\bibitem{Hai2006}
E.~Hairer.
\newblock Long-time energy conservation of numerical integrators.
\newblock In {\em Foundations of computational mathematics, {S}antander 2005},
  volume 331 of {\em London Math. Soc. Lecture Note Ser.}, pages 162--180.
  Cambridge Univ. Press, Cambridge, 2006.

\bibitem{Hua10}
W.~Huang.
\newblock Discrete maximum principle and a delaunay-type mesh condition for
  linear finite element approximations of two-dimensional anisotropic diffusion
  problems.
\newblock {\em Numer. Math. Theory Meth. Appl.}, 4:319--334, 2011.
\newblock (arXiv:1008.0562).

\bibitem{KK09}
J.~Kar{\'a}tson and S.~Korotov.
\newblock An algebraic discrete maximum principle in {H}ilbert space with
  applications to nonlinear cooperative elliptic systems.
\newblock {\em SIAM J. Numer. Anal.}, 47:2518--2549, 2009.

\bibitem{KKK07}
J.~Kar\'atson, S.~Korotov, and M.~K\v{r}\'i\v{z}ek.
\newblock On discrete maximum principles for nonlinear elliptic problems.
\newblock {\em Math. Comput. Sim.}, 76:99--108, 2007.

\bibitem{KSS09}
D.~Kuzmin, M.~J. Shashkov, and D.~Svyatskiy.
\newblock A constrained finite element method satisfying the discrete maximum
  principle for anisotropic diffusion problems.
\newblock {\em J. Comput. Phys.}, 228:3448--3463, 2009.

\bibitem{KL95}
M.~K\v{r}\'i\v{z}ek and Q.~Lin.
\newblock On diagonal dominance of stiffness matrices in {3D}.
\newblock {\em East-West J. Numer. Math.}, 3:59--69, 1995.

\bibitem{LePot05}
C.~Le~Potier.
\newblock Sch\'ema volumes finis monotone pour des op\'erateurs de diffusion
  fortement anisotropes sur des maillages de triangles non structur\'es.
\newblock {\em C. R. Math. Acad. Sci. Paris}, 341:787--792, 2005.

\bibitem{LePot09b}
C.~Le~Potier.
\newblock Un sch\'ema lin\'eaire v\'erifiant le principe du maximum pour des
  op\'erateurs de diffusion tr\`es anisotropes sur des maillages d\'eform\'es.
\newblock {\em C. R. Math. Acad. Sci. Paris}, 347:105--110, 2009.

\bibitem{Let92}
F.~W. Letniowski.
\newblock Three-dimensional {D}elaunay triangulations for finite element
  approximations to a second-order diffusion operator.
\newblock {\em SIAM J. Sci. Stat. Comput.}, 13:765--770, 1992.

\bibitem{Lew1994}
M.~A. Lewis.
\newblock Spatial coupling of plant and herbivore dynamics: The contribution of
  herbivore dispersal to transient and persistent ``waves" of damage.
\newblock {\em Theoret. Population Biol.}, 45:277--312, 1994.

\bibitem{LH10}
X.~P. Li and W.~Huang.
\newblock An anisotropic mesh adaptation method for the finite element solution
  of heterogeneous anisotropic diffusion problems.
\newblock {\em J. Comput. Phys.}, 229:8072--8094, 2010 (arXiv:1003.4530v2).

\bibitem{LiHu2012}
X.~P. Li and W.~Huang.
\newblock Maximum principle for the finite element solution of time dependent
  anisotropic diffusion problems.
\newblock {\em Numer. Meth. P.D.E.}, 2012.
\newblock (arXiv:1209.5657).

\bibitem{LSS07}
X.~P. Li, D.~Svyatskiy, and M.~Shashkov.
\newblock Mesh adaptation and discrete maximum principle for {2D} anisotropic
  diffusion problems.
\newblock Technical Report LA-UR 10-01227, Los Alamos National Laboratory, Los
  Alamos, NM, 2007.

\bibitem{LSSV07}
K.~Lipnikov, M.~Shashkov, D.~Svyatskiy, and Y.~Vassilevski.
\newblock Monotone finite volume schemes for diffusion equations on
  unstructured triangular and shape-regular polygonal meshes.
\newblock {\em J. Comput. Phys.}, 227:492--512, 2007.

\bibitem{LS08}
R.~Liska and M.~Shashkov.
\newblock Enforcing the discrete maximum principle for linear finite element
  solutions of second-order elliptic problems.
\newblock {\em Comm. Comput. Phys.}, 3:852--877, 2008.

\bibitem{LuHuQi2012}
C.~Lu, W.~Huang, and J.~Qiu.
\newblock Maximum principle in linear finite element approximations of
  anisotropic diffusion-convection-reaction problems.
\newblock {\em Numer. Math.} in press, 2012 (arXiv:1201.3564v1).

\bibitem{MD06}
M.~J. Mlacnik and L.~J. Durlofsky.
\newblock Unstructured grid optimization for improved monotonicity of discrete
  solutions of elliptic equations with highly anisotropic coefficients.
\newblock {\em J. Comput. Phys.}, 216:337--361, 2006.

\bibitem{MM05}
K.~W. Morton and D.~F. Mayers.
\newblock {\em Numerical Solution of Partial Differential Equations}.
\newblock Cambridge University Press, Cambridge, U.K., 2005.
\newblock 2nd Edition.

\bibitem{RuWiXu2007}
R.~D. Russell, J.~F. Williams, and X.~Xu.
\newblock M{OVCOL}4: a moving mesh code for fourth-order time-dependent partial
  differential equations.
\newblock {\em SIAM J. Sci. Comput.}, 29:197--220 (electronic), 2007.

\bibitem{SH07}
P.~Sharma and G.~W. Hammett.
\newblock Preserving monotonicity in anisotropic diffusion.
\newblock {\em J. Comput. Phys.}, 227:123--142, 2007.

\bibitem{Sto86}
G.~Stoyan.
\newblock On maximum principles for monotone matrices.
\newblock {\em Lin. Alg. Appl.}, 78:147--161, 1986.

\bibitem{SF73}
G.~Strang and G.~J. Fix.
\newblock {\em An Analysis of the Finite Element Method}.
\newblock Prentice Hall, Englewood Cliffs, NJ, 1973.

\bibitem{SunRusXu2007}
P.~Sun, R.~D. Russell, and J.~Xu.
\newblock A new adaptive local mesh refinement algorithm and its application on
  fourth order thin film flow problem.
\newblock {\em J. Comput. Phys.}, 224:1021--1048, 2007.

\bibitem{Teman1968}
R.~Temam.
\newblock Une m\'ethode d'approximation de la solution des \'equations de
  {N}avier-{S}tokes.
\newblock {\em Bull. Soc. Math. France}, 96:115--152, 1968.

\bibitem{TW08}
V.~Thom{\'e}e and L.~B. Wahlbin.
\newblock On the existence of maximum principles in parabolic finite element
  equations.
\newblock {\em Math. Comput.}, 77:11--19, 2008.

\bibitem{WaZh11}
J.~Wang and R.~Zhang.
\newblock Maximum principle for {P1}-conforming finite element approximations
  of quasi-linear second order elliptic equations.
\newblock {\em SIAM J Numer. Anal.}, 50:626--642, 2012.
\newblock (arXiv:1105.1466).

\bibitem{XZ99}
J.~Xu and L.~Zikatanov.
\newblock A monotone finite element scheme for convection-diffusion equations.
\newblock {\em Math. Comput.}, 69:1429--1446, 1999.

\bibitem{ZhSh2010}
X.~Zhang and C.-W. Shu.
\newblock On maximum-principle-satisfying high order schemes for scalar
  conservation laws.
\newblock {\em J. Comput. Phys.}, 229:3091--3120, 2010.

\bibitem{ZhoBer2000}
L.~Zhornitskaya and A.~L. Bertozzi.
\newblock Positivity-preserving numerical schemes for lubrication-type
  equations.
\newblock {\em SIAM J. Numer. Anal.}, 37:523--555, 2000.

\end{thebibliography}