A Computational Study of Residual KPP Front Speeds in Time-Periodic Cellular Flows in the Small Diffusion Limit
Penghe Zu
Long Chen
Jack Xin \thanks{Department of Mathematics
UC Irvine
Irvine
CA 92697
USA.\newline Corresponding author: Penghe Zu. Phone number: 9493784420.\newline \hspace{.2 in
The minimal speeds ($c^*$) of the Kolmogorov-Petrovsky-Piskunov (KPP) fronts at small diffusion ($\epsilon \ll 1$) in a class of time-periodic cellular flows with chaotic streamlines is investigated in this paper. The variational principle of $c^*$ reduces the computation to that of a principal eigenvalue problem on a periodic domain of a linear advection-diffusion operator with space-time periodic coefficients and small diffusion. To solve the advection dominated time-dependent eigenvalue problem efficiently over large time, a combination of finite element and spectral methods, as well as the associated fast solvers, are utilized to accelerate computation. In contrast to the scaling $c^*=\mathcal{O(\epsilon^{1/4)$ in steady cellular flows, a new relation $c^* = \mathcal{O(1)$ as $\epsilon \ll 1$ is revealed in the time-periodic cellular flows due to the presence of chaotic streamlines. Residual propagation speed emerges from the Lagrangian chaos which is quantified as a sub-diffusion process.
ZuCX.tex
\documentclass[10pt]{article}
\usepackage{colordvi,amsfonts,amsmath,epsfig,subfigure,cite}
%\usepackage[margin=3.68cm]{geometry}
%\usepackage[pdftex,colorlinks=true]{hyperref}
%  \newcounter{mnote}
%  \setcounter{mnote}{0}
%  \newcommand{\mnote}[1]{\addtocounter{mnote}{1}
%    \ensuremath{{}^{\bullet\arabic{mnote}}}
%    \marginpar{\footnotesize\em\color{red}\ensuremath{\bullet\arabic{mnote}}#1}}
%  \let\oldmarginpar\marginpar
%    \renewcommand\marginpar[1]{\-\oldmarginpar[\raggedleft\footnotesize #1]%
%    {\raggedright\footnotesize #1}}


%\setlength{\oddsidemargin}{-0.08in}
%\setlength{\textheight}{9.0in}
%\setlength{\textwidth}{6.5in}
%\setlength{\topmargin}{- in}
%\usepackage{amsmath,latexsym,amssymb,graphicx,epsf,amsfonts,amsthm}
%\usepackage{epsf}
\usepackage{color}
\newcommand{\PZ}[1]{\textcolor{red}{#1}}
\newcommand{\LC}[1]{\textcolor{cyan}{#1}}
%\usepackage{mathrsfs}
%\usepackage{epstopdf}
%\usepackage{epsfig}
%\usepackage{subfigure}
%\usepackage[toc,page]{appendix}
%\usepackage{framed}
\numberwithin{equation}{section} %\numberwithin{figure}{section}
\newcommand{\ra}{\rightarrow}
\newcommand{\wh}{\widehat}
\newcommand{\nit}{\noindent}
\newcommand{\no}{\nonumber}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\ba}{\begin{eqnarray}}
\newcommand{\ea}{\end{eqnarray}}
\newcommand{\bi}{\begin{itemize}}
\newcommand{\ei}{\end{itemize}}
\newcommand{\bc}{boundary condition}
\newcommand{\eq}{equation}
\newcommand{\dta}{\mbox{$\delta$}}
\newcommand{\fee}{\mbox{$\varphi$}}
\newcommand{\ar}{\mbox{$\alpha$}}
\newcommand{\lam}{\mbox{$\lambda$}}
\newcommand{\eps}{\mbox{$\epsilon$}}
\newcommand{\gam}{\mbox{$\gamma$}}
\newcommand{\rn}{\mbox{$R^{n}$}}
\newcommand{\bk}{\mbox{$k$}}
\newcommand{\DF}[2]{\frac{\D#1}{\D#2}}
\newtheorem{thm}{Theorem}
\newtheorem{theorem}{Theorem}
\newtheorem{pro}[thm]{Proposition}
\newtheorem{lemma}[thm]{Lemma}
\newtheorem{cor}[thm]{Corollary}
\newtheorem{corollary}[thm]{Corollary}
\newtheorem{alg}[thm]{Algorithm}
\newtheorem{algorithm}[thm]{Algorithm}
\newtheorem{ass}{Assumption}
\newtheorem{defi}[thm]{Definition}
\newtheorem{rem}[thm]{Remark}
%\textwidth  = 15.3truecm \textheight = 24truecm
\begin{document}
%\hoffset = -1truecm \voffset = -1truecm

\title{A Computational Study of Residual KPP Front Speeds in Time-Periodic Cellular Flows 
in the Small Diffusion Limit}
\author{Penghe Zu, Long Chen, Jack Xin \thanks{Department of Mathematics, 
UC Irvine, Irvine, CA 92697, USA.\newline Corresponding author: Penghe Zu. Phone number: 9493784420.\newline \hspace{.2 in} Emails: (pzu,chenlong,jxin)@math.uci.edu.}}
\date{}
\maketitle
%******************************************************************************************************
% abstract
%******************************************************************************************************
\begin{abstract}
The minimal speeds ($c^*$) of the Kolmogorov-Petrovsky-Piskunov (KPP) fronts  
at small diffusion ($\epsilon \ll 1$) in a class of time-periodic cellular flows with chaotic streamlines is investigated in this paper. 
The variational principle of $c^*$ reduces the computation to 
that of a principal eigenvalue problem on a periodic domain of 
a linear advection-diffusion operator with space-time periodic coefficients and small diffusion.
To solve the advection dominated time-dependent eigenvalue problem efficiently over large time, 
a combination of finite element and spectral methods, as well as the associated 
fast solvers, are utilized to accelerate computation. 
In contrast to the scaling $c^*=\mathcal{O}(\epsilon^{1/4})$ in steady cellular flows, 
a new relation $c^* = \mathcal{O}(1)$ as $\epsilon \ll 1$ is revealed in the time-periodic cellular flows due to the presence of chaotic streamlines. 
Residual propagation speed emerges from the Lagrangian chaos which is quantified as a sub-diffusion process.
\end{abstract}
\medskip

\providecommand{\keywords}[1]{\textbf{\textit{Keywords:}} #1}
\keywords{KPP front speeds, time periodic cellular flows, chaotic streamlines, sub-diffusion, residual front speeds.\\
\\
{\bf PACS}: 42.25.Dd, 02.50.Fz, 05.10.-a, 02.60.Cb}
\thispagestyle{empty}
%\newpage
\setcounter{equation}{0} \setcounter{page}{0}
\newpage
%******************************************************************************************************
% Introduction
%******************************************************************************************************
\section{Introduction}
Front propagation in complex fluid flows arises in many scientific 
areas such as combustion~\cite{R,W85}, population
growth of ecological communities (plankton) in the ocean~\cite{A}, and 
reactive chemical front in liquids~\cite{PS05,X09}. 
An interesting problem is the front speed enhancement 
in time dependent fluid flows with chaotic streamlines and random media~\cite{PS05,X09}. 

In this paper, we shall consider a two dimensional time dependent cellular flow:
\be
(\cos(2\pi y) + \theta \sin(2\pi y)\cos(t), \, \cos(2\pi x) + \theta \sin(2\pi x)\cos(t)), \; \; 
\theta \in (0,1], \label{tdcell}
\ee 
whose steady part $(\cos (2\pi y), \cos(2\pi x))$ is subject to 
time periodic perturbation that causes transition to Lagrangian chaos. 
Chaotic behavior of a similar flow field is qualitatively analyzed by formal dynamical system 
methods in~\cite{CW91}. The enhanced residual diffusion in (\ref{tdcell}) is observed numerically in~\cite{BCVV95}. 
The enhanced diffusion and propagation speeds in steady cellular flows with ordered streamlines and 
their extensions have been extensively studied~\cite{ACVV02,Aud,ChSo,FP,Nov11,H,NR05,RZ07,SZX13,Liu_13,Yu1,Yu2} among others.  
%
We shall take a statistical look at the 
chaotic streamlines of (\ref{tdcell}) in terms of the scaling of 
the mean square displacements from random initial data, and  
quantify the Lagrangian chaos of (\ref{tdcell}) as a sub-diffusion process. 

Consider then the advection-reaction-diffusion equation:
\begin{equation}\label{eq:nonlinear1}
\partial_t u = \epsilon \triangle u +\vec{B}(x,t)\cdot\nabla u + \frac{1}{\tau}f(u), \quad x\in \mathbb{R}^2, t>0, 
\end{equation}
where $f(u)=u(1-u)$ is the Kolmogorov-Petrovsky-Piskunov (KPP) nonlinearity,
$\epsilon$ is the molecular diffusion parameter,  
$\tau$ is the reaction rate and $\vec{B}(x,t)$ is a space-time periodic, mean zero, and incompressible 
flow field such as (\ref{tdcell}).  If the initial data for $u$ is nonnegative and compactly supported, 
the large time behavior of $u$ is an outward propagating front,  
with speed $c^*=c^*(\vec{e})$ in the direction $\vec{e} = \langle1,0\rangle$. The variational principle of $c^*$ is~\cite{NX05}:
\begin{equation}\label{eq:inf}
c^*=\inf_{\lambda>0} \frac{\mu(\lambda)}{\lambda},
\end{equation}
where $\mu(\lambda)$ is the principal eigenvalue of the periodic-parabolic operator:
\begin{equation}
L^{\lambda}\Phi = \epsilon\triangle\Phi + (\vec{B}+2\lambda \vec{e})\cdot \nabla\Phi + 
(\epsilon \lambda^2 + \lambda \vec{B}\cdot \vec{e} + \tau^{-1}\, f'(0))\Phi -\Phi_t, \label{varp1}
\end{equation}
on the space-time periodic cell $\Omega$.
The principal eigenvalue $\mu(\lambda)$ can be computed by solving the evolution problem:
\ba
&& w_t  =  \epsilon\triangle w +(2\epsilon\lambda\vec{e} + \vec{B}(x,t))\cdot\nabla w + 
 (-C_M +\epsilon\lambda^2 +\lambda\vec{e}\cdot \vec{B}(x,t) + \frac{1}{\tau}f'(0)) w \nonumber \\
&& w(x,0) = 1. \label{para}
\ea
Here the constant $C_M = \epsilon\lambda^2 +\lambda||\vec{B}||_{\infty} 
+ ||\frac{1}{\tau}f'(0)||_{\infty}$, 
so that $w$ remains positive and bounded by one.  The $\mu (\lambda)$ is then given by:
\begin{equation}\label{lyapunov}
\mu(\lambda) = C_M + \lim_{t\rightarrow\infty} \frac{1}{t} \ln\int_{\Omega} w(x,t)\, dx
\end{equation}
The number $\mu(\lambda)$ is also the principle 
Lyapunov exponent of the parabolic equation (\ref{para}), and the formula 
(\ref{lyapunov}) extends to the more general case when $\vec{B}$ is a stationary 
ergodic field\cite{NX09}. The limit then holds almost surely and 
$\mu$ is deterministic\cite{NX09}. KPP fronts are examples of the so called ``pulled fronts''~\cite{S03} 
because their speed is determined by the behavior of the solution far 
beyond the front interface, in the region where the solution is close to zero 
(the unstable equilibrium). The minimal speed of the planar wave solution of the linearized equation 
near the unstable equilibrium gives (\ref{eq:inf}), also known as the marginal stability criterion\cite{S03,X09}.  

Equation (\ref{para}) and formula (\ref{lyapunov}) 
will be discretized for approximating $\mu$ at sufficiently large time for a range of 
$\lambda $ values. The minimal point of $\mu (\lambda)/\lambda $ is searched by the 
golden section algorithm. At each new search of $\lambda$, 
the large time solution of (\ref{para}) is computed. Because $\epsilon$ is small, 
(\ref{para}) is advection dominated. To this end, 
upwinding type finite element methods (EAFE~\cite{XZ99}) and spectral method
(SM) are utilized. Spectral method has high accuracy but the computation is slower 
because it takes small time step due to the restriction of stability. On the other
hand, EAFE can be solved faster using relative larger time step with
implicit discretization. Without the compromise of accuracy, our strategy to obtain 
the minimal point is to narrow down the search interval by EAFE first, and solve
more accurately in a smaller interval by SM.

%Our first finding is that the chaotic streamlines of the time periodic cellular flow (\ref{tdcell}) 
%can be quantified as a sub-diffusive process with distinct $\theta$-dependent scaling exponents. 
%The resulting motion appears ergodic inside the invariant infinite channel domains. Our next finding is  
%the front speed phenomenon $c^*(\epsilon) = O(1)$ as $\epsilon \downarrow 0$. In contrast, 
%$c^{*} = O(\epsilon^{1/4})$ in steady cellular flows \cite{NR05}. 
%The streamlines of the steady cellular flows are ordered with closed orbits 
%leading to a much slower rate of enhancement for the transport. We also observe the presence of 
%layer and circular structures in the generalized eigenfunction $w$ at large time, as a 
%reflection of the advection-dominated transport in time periodic cellular flow (\ref{tdcell}).

The rest of the paper is organized as follows. In section 2, we illustrate the 
difference of streamline geometry and dynamical properties of 
steady and unsteady cellular flows. We find that 
the chaotic streamlines of the time periodic cellular flow (\ref{tdcell}) 
can be quantified as a sub-diffusive process with distinct $\theta$-dependent scaling exponents. 
The resulting motion appears ergodic inside the invariant infinite channel domains.
In section 3, we discuss numerical methods for advection-dominated problems 
such as EAFE, and semi-implicit SM, as well as the stability of 
these methods and time step constraints. 
In section 4, we show numerical results of KPP front speeds by 
a combination of these methods and the existence of  
residual front speed in the sense that $c^*(\epsilon) = O(1) > 0$ as $\epsilon \downarrow 0$. In contrast, 
$c^{*} = O(\epsilon^{1/4})$ in steady cellular flows \cite{NR05}. 
The streamlines of the steady cellular flows are ordered with closed orbits 
leading to a much slower rate of enhancement for the transport. We also observe the presence of 
layer and circular structures in the generalized eigenfunction $w$ at large time, a 
reflection of the advection-dominated transport in time periodic cellular flow (\ref{tdcell}). The 
corresponding energy spectrum of $w$ shows decay towards high frequencies with an intermediate scaling range. 
In section 5, we give concluding remarks about our findings.

%*******************************************************************************************************
% Section 2
%*******************************************************************************************************
\section{Properties of Cellular Flows}
 In this section, we illustrate the 
difference of streamline geometry and dynamical properties of 
steady and unsteady cellular flows. In particular, we quantify the chaotic 
behavior of the motion along streamlines of (\ref{tdcell}) using the empirical 
mean square distance inside infinite invariant channels 
in the direction $(1,1)$, or $E[|X(t)\cdot (1,1)|^2]$, where 
$X\in \mathbb R^2$ denotes the particle trajectories in the flow. 
We show computationally that it scales 
with time as $O(t^p)$, $p \in (0,1)$, for sampled values of $\theta \in (0,1]$.
\subsection{Steady Cellular Flow}
A phase portrait of the steady cellular flow 
\be
B_s=(\cos (2\pi y), \cos (2\pi x)), \label{steady}
\ee
with a cell square is presented in Fig.~\ref{stcell}. 
The phase portrait away from the square simply repeats. The streamlines are ordered. The 
closed orbits form elliptic part of the phase space. The 
saddles are located at half-integer points $(n,m)/2$, $n$,$m \in \mathbb{Z}$ 
with connecting separatrices forming hyperbolic part of the phase space. 
Any particle trajectory is either a closed orbit or a separatrix, hence the motion is 
bounded inside the square.
  
\begin{figure}
\includegraphics[width=1 \textwidth, height = 0.8\textwidth]{scell.eps}
\caption{Ordered and localized streamlines of steady cellular flow (\ref{steady}) with a cell square containing 
four vertices.}\label{stcell}
\end{figure}

\subsection{Unsteady Cellular Flow}
For the unsteady flow 
\be
B_u=(\cos(2\pi y) + \theta \, \sin(2\pi y)\cos(t), \cos(2\pi x) + \theta \, \sin(2\pi x)\cos(t)),\label{unsteady}
\ee 
the invariant manifold consists of lines $y=x+n$, $n \in \mathbb{Z}$. 
The flow trajectories are restricted in the channels bounded by two neighboring lines 
$y=x+n+1$ and $y =x+n$. At $\theta > 0$, the flow trajectory extends itself from one cell square 
to another. The Lagrangian particle undergoes chaotic motion, see Fig.~\ref{cntFlow} for an 
illustration. The computation is done by a fourth order symplectic scheme.
The local flow direction is colored red. The intensified red region in the middle indicates 
that the particle spends a lot of time wandering in and out of the cells there. 
Fig.~\ref{projTr} plots a projected trajectory in the direction $(1,1)/\sqrt{2}$ vs. time. 
The stochastic feature is visible.  

\begin{figure}[htdp]
\includegraphics[width=0.5\textwidth, height = 0.45\textwidth]{tdTrajFinTime20.eps}
\includegraphics[width=0.5\textwidth, height = 0.45\textwidth]{tdTrajFinTime200.eps}
\caption{Disordered and extended trajectories in time periodic cellular flow (\ref{unsteady})  
over time $[0,20]$ (left) and $[0,200]$ (right) in a channel bounded by invariant lines $y=x+1$ and $y=x$. The local flow direction is colored red.}\label{cntFlow}
\end{figure}

\begin{figure}[htdp]
\includegraphics[width=1 \textwidth, height = 0.8\textwidth]{projRandTraj.eps} 
\caption{A projected trajectory $X(t)\cdot (1,1)/\sqrt{2}$ in time periodic cellular flow (\ref{unsteady}), 
$\theta =1$, showing stochastic character.}\label{projTr}
\end{figure}

To quantify the disorder and draw a connection with diffusion process, 
we compute $10^4$ trajectories $X(t,\omega)$ over time interval $[0,1000]$ with uniformly distributed initial points on a base cell square, 
where $\omega$ denotes the random samples. We calculate the empirical mean square distance ($E[|X(t,\omega)|^2]$)
and the projected mean square distance 
$E[ (X(t,\omega)\cdot (1,1))^2 ]$, and plot them as a function of time on the logarithmic 
scale to recover scaling laws. For efficiency, the samples are obtained by a 500 node parallel implementation of 
the standard 4th order Runge-Kutta method. Fig.~\ref{subd1} shows that the mean square distances  
$E[|X(t,\omega)|^2] \sim O(t^p)$, $p\approx 0.483, 0.673$ at $\theta =0.1, 0.4$ respectively, hence belonging to the sub-diffusive regime. 
Similar sub-diffusive behavior can be observed in Fig.~\ref{subd2} where the projected mean square distance scales like  
$E[ (X(t,\omega)\cdot (1,1))^2 ]\sim O(t^q)$, $q \approx 0.605, 0.786$ at $\theta =0.1, 0.4$.

\begin{figure}[ht!]
\includegraphics[width=0.5\textwidth, height = 0.45\textwidth]{var_theta=1e-1.eps}
\includegraphics[width=0.5\textwidth, height = 0.45\textwidth]{var_theta=4e-1.eps}
\caption{Mean square distances $E[|X(t,\omega)|^2]$ vs. time on the logarithmic scale with robust fit showing sub-diffusive scaling $O(t^p)$,
$p=0.483, 0.673 $ at $\theta=0.1$ (left) and $0.4$ (right).}\label{subd1}
\end{figure}

\begin{figure}[ht!]
\includegraphics[width=0.5\textwidth, height = 0.45\textwidth]{dot_theta=1e-1.eps}
\includegraphics[width=0.5\textwidth, height = 0.45\textwidth]{dot_theta=4e-1.eps}
\caption{Projected mean square distances $E[ (X(t,\omega)\cdot (1,1))^2 ]$ vs. time on the logarithmic scale with robust fit 
showing sub-diffusive scaling $O(t^q)$, $q=0.605, 0.786$ at $\theta=0.1$ (left) and $0.4$ (right).}\label{subd2}
\end{figure}




%*******************************************************************************************************
% Section 3
%*******************************************************************************************************
\section{Numerical Methods}
Numerical approximation of advection-dominated problems is a 
topic of independent interest. We shall use the 
Edge-Averaged Finite Element (EAFE) method~\cite{XZ99} due to the nice  
discrete maximum principle obeyed by EAFE and the corresponding fast multigrid solvers. 
On the other hand, pseudo-spectral methods are a class of highly accurate numerical methods for 
solving partial differential equations. 
In practice, pseudo-spectral method has excellent error reduction properties 
with the so-called ``exponential convergence'' being the fastest possible 
as long as the solution is smooth. 
In this section, we briefly review these two methods used in our simulation.

%We have compared Finite Volume method~\cite{BBFS}, Edge-Averaged Finite Element (EAFE) method~\cite{XZ99} and 
%Streamline Diffusion (SD) method~\cite{ESW} for solving equation (\ref{para}). 
%The advantage of the first two schemes is that the discrete maximum principle  
%holds and consequently the positivity of numerical solutions is preserved. Also 
%since we are concerned with the long time behavior of the numerical solution, huge amount of 
%work on solving sparse matrix is inevitable. So we will use more advanced solver to accelerate this process.

%In this section, we discuss {\bf four} methods, namely the Box method~\cite{BBFS}, EAFE method~\cite{XZ99}, streamline diffusion method~\cite{ESW} and 
%pseudo-spectral method, to solve 
Recall that we are solving the following evolution 
problem over the unit square domain $\Omega = [0,1]\times[0,1]$ with 
periodic boundary condition and a constant initial condition:
\begin{equation}\label{eq:linear5}
\left\{
\begin{array}{l}
w_t = \epsilon\triangle w +(2\epsilon\lambda\vec{e} + B)\cdot\nabla w - (C_M - \epsilon\lambda^2 -\lambda\vec{e}\cdot B - \frac{1}{\tau}f'(0)) w,\\
w(0,y,t)=w(1,y,t); w(x,0,t) = w(x,1,t), \\
w(x,y,0)=1.
\end{array}
\right.
\end{equation}

%For the first three finite element schemes, we have tested them on this triangulation. The numerical result shows that among these three schemes, streamline diffusion method has the highest accuracy. {\bf EAFE method is slightly more accurate than Box method.} And the matrix of the linear system assembled by EAFE is M-matrix. In other words, the matrix in EAFE method has zero column sums, with nonnegative diagonal and nonpositive off diagonal entries. The solution of the equation will be nonnegative if the right hand side vector is nonnegative. This property is also called discrete maximum principle.
%
%For the above reason, we choose EAFE method instead of Box method or streamline diffusion method to compute the advection dominant problem. 



%The Box method was first proposed in~\cite{BBFS} and was formally 
%derived as a finite volume approximation of advection-dominated problem. 
%In~\cite{BBFS}, the authors pointed out that as long as the flow function $\vec{B}$ is divergence free, the linear equation formed by discretizing the system will be an M-matrix. In other words, the matrix in Box scheme has zero column sums, with nonnegative diagonal and nonpositive off diagonal entries. The solution of the equation will be nonnegative if the right hand side vector is nonnegative. This property is also called discrete maximum principle.

%Another discrete maximal principle preserving scheme is~\cite{XZ99}. 
%It requires that in $R^2$ the sum of angles opposite to any edge be less than or 
%equal to $\pi$. This condition implies the Delaunay triangulation.

%The streamline diffusion method~\cite{BH82} does not preserve 
%the discrete maximal principle. However, it has second order accuracy and 
%even high order if the stabilizer parameter is chosen properly. 
%The streamline diffusion scheme will be combined with the 
%pseudo-spectral method to compute $c^*$.  {\bf Clarify where we used Box or EAFE: We compare these finite element method before. After computing the bench mark example, the result shows that EAFE has higher accuracy than Box method. On the other hand, even though SD method has the highest accuracy among those three method, but it losses the discrete maximum principle. Based on the numerical result, we choose EAFE method as the preliminary method to narrow down the searching interval of Golden-section method}

%\subsection{Box Method}
%For discretization, we write the differential operator of (\ref{para}) in the form:
%\begin{equation}\label{eq:operator}
%\mathcal{L} u = -\nabla\cdot(\nabla u + \beta u).
%\end{equation}
%Based on the calculation in~\cite{BBFS} we have the following formula for stiff matrix:
%\begin{equation}
%\begin{array}{ll}
%A_b =& (1+\frac{l_1|\beta\cdot t_1|}{2})\frac{s_1}{l_1}
%\left[
%\begin{array}{rrr}
%0&0&0\\
%0&1&-1\\
%0&-1&1
%\end{array}
%\right]
%%%
%+(1+\frac{l_2|\beta\cdot t_2|}{2})\frac{s_2}{l_2}
%\left[
%\begin{array}{rrr}
%1&0&-1\\
%0&0&0\\
%-1&0&1
%\end{array}
%\right] \\
%%%
%&+(1+\frac{l_3|\beta\cdot t_3|}{2})\frac{s_3}{l_3}
%\left[
%\begin{array}{rrr}
%1&-1&0\\
%-1&1&0\\
%0&0&0
%\end{array}
%\right]
%%%
%+\frac{s_1\beta\cdot t_1}{2}
%\left[
%\begin{array}{rrr}
%0&0&0\\
%0&-1&-1\\
%0&1&1
%\end{array}
%\right] \\
%%%
%&+\frac{s_2\beta\cdot t_2}{2}
%\left[
%\begin{array}{rrr}
%1&0&1\\
%0&0&0\\
%-1&0&-1
%\end{array}
%\right]
%%%
%+\frac{s_3\beta\cdot t_3}{2}
%\left[
%\begin{array}{rrr}
%-1&-1&0\\
%1&1&0\\
%0&0&0
%\end{array}
%\right]
%\end{array}
%\end{equation}
%where $l_i$ represents length of $i$-th edge, $t_i$ is tangent vector of 
%each edge counterclockwise and $s_i$ is the length of box edge corresponding 
%to the perpendicular bisector of $i$-th edge.

\subsection{EAFE Method}
We first decompose the domain $\Omega$ into a triangulation $\mathcal T$ which is a set of triangles such that 
\begin{equation}
\cup_{\tau\in\mathcal{T}} \tau = \bar{\Omega}\quad \text{and}\quad \mathring{\tau}_i\cap\mathring{\tau}_j = \emptyset.
\end{equation} 
For the unit square, we simply use the uniform mesh obtained by setting length of each triangle $h$. 

We describe the EAFE method using a simple advection-diffusion equation
\begin{equation}\label{eq:CReq}
\left\{
\begin{array}{ll}
-\nabla\cdot(\nabla u + \beta(x) u) = f(x), &x\in\Omega\\
u=0 & x\in \partial\Omega.
\end{array}
\right.
\end{equation}
Associated with each $\mathcal{T}_h$, let $V_h\subset H^1_0(\Omega)$ be the piecewise 
linear finite element space. The space $H^1_0(\Omega)$ is defined as the subspace of $H^1(\Omega)$ 
with zero trace on $\partial\Omega$.

Given any edge $E$ in the triangulation, we introduce a function $\phi_E$ defined locally on $E$ 
(up to an arbitrary constant) by the relation:
\begin{equation}
\frac{\partial\phi_E}{\partial t_E} = \frac{1}{|t_E|} \epsilon^{-1}(\beta\cdot t_E),
\end{equation}
where $E$ is the edge connecting two vertices $q_i$ and $q_j$, and $t_E$ is the vector such that $t_E = q_i-q_j$.
	
The EAFE formulation of the problem (\ref{eq:CReq}) is: Find $u_h\in V_h$ such that
\begin{equation}
a_h(u_h,v_h) = f(v_h)\quad \text{for}\ \text{any}\ v_h\in V_h
\end{equation}
where
\begin{equation*}
a_h(u_h,v_h) = \sum_{T\in\mathcal{T}_h}
\left\{
\sum_{E\subset T} 
\frac{\omega^T_E}{|t_E|} \int_E e^{\phi_E}J(u)\cdot t_Eds
\ \delta_E v_h
\right\}
\end{equation*}
and
\begin{equation*}
\left\{
\begin{array}{l}
\omega_E = \frac{1}{2} \sum_{E\subset T}\cot\theta_E^T \geq 0\\
J(u) = \epsilon\nabla u + \beta u
\end{array}
\right.
\end{equation*}
where $\theta_E^T$ is 
the angle between the edges that share the common edge $E$ and $\delta_E$ is related to the tangential derivative along $E$.

The authors of~\cite{XZ99} show that if the triangulation is a so-called Delaunay triangulation, i.e., the summation of two angles of an interior edge is less than or equal to $\pi$, the matrix of the linear system generated by EAFE 
would be an M-matrix with nonnegative row sum. The detail of the scheme can be found in~\cite{XZ99}.

When applied to equation \eqref{eq:linear5}, notice that $\vec{B}$ is divergence free, the operator can be written in the form of \eqref{eq:CReq} and the reaction term is always positive. We apply the mass lumping method~\cite{TPB71} so that the discretization of the reaction term makes positive contribution to the diagonal. Therefore the M-matrix property and consequently the discrete maximum principle still holds. 

\begin{rem}
We have tested the Box method~\cite{BBFS}, 
EAFE method~\cite{XZ99} and streamline diffusion method~\cite{ESW} on the advection-dominated problem. 
The accuracy of EAFE scheme is better than Box method, but lower than the streamline diffusion method. 
We chose EAFE instead of streamline diffusion scheme because  
the matrix of linear system generated by EAFE is M-matrix. 
The M-matrix not only preserves the discrete maximum principle 
(so that the numerical solution stays between 0 and 1) but also benefits from the 
fast solvers.   
\end{rem}
%\subsection{Streamline diffusion method}
%The streamline diffusion method~\cite{BH82} gives:
%\begin{framed}
%Find $u_h\in V_h$ such that
%\begin{equation}
%a_{sd}(u_h, v_h) = l_{sd}(v_h)\quad \text{for all } v_h\in V_h^0
%\end{equation}
%where
%\begin{equation}\label{eq:streamline}
%\begin{array}{ll}
%a_{sd}(u,v) =& \epsilon\int_{\Omega}\nabla u\cdot\nabla v + \int_{\Omega}(\vec{B}\cdot\nabla u)v +\delta\int_{\Omega}(\vec{B}\cdot\nabla u)(\vec{B}\cdot\nabla v)\\
%l_{sd}(v) =& \int_{\Omega} fv+\delta\int_{\Omega}f\vec{B}\cdot\nabla v
%\end{array}
%\end{equation}
%\end{framed}

%To implement (\ref{eq:streamline}), we  
%choose the parameter $\delta$, or locally parameter $\delta_k$ as:
%\begin{equation}
%\delta_k=
%\left\{
%\begin{array}{ll}
%\frac{h_k}{2|\vec{B}|} (1-\frac{1}{Pe_k}) &\text{if } Pe_k >1,\\
%0 & \text{if } Pe_k\leq 1.
%\end{array}
%\right.
%\end{equation}
%Here $|\vec{B}|$ is the $l^2$ norm of the wind at the element centroid, $h_k$ is the measure of element length in the direction of the wind, and $Pe_k=|\vec{B}|h_k/(2\epsilon)$ is the element Peclet number.

To accelerate the speed of computation, Algebraic Multi-Grid (AMG) solver is
used to solve the linear algebraic equation arising from each time step. When the matrix is an M-matrix, the corresponding AMG solver is proven to be efficient. Specifically for the advection-diffusion equations, a multigrid method preserving the M-matrix property in coarse level is developed in \cite{KXZ03}. 
Among many AMG software packages, we use AGMG (aggregation-based algebraic multigrid method) package~\cite{N}.
%  
%AMG is a fast solver for large systems of linear equations~\cite{N}. 
%{\bf It has high performance} on large scale sparse linear systems arising 
%from upwind finite difference discretization of advection-diffusion equations. 
The numerical test shows that AMG solver is $8$ times faster than the direct solver.     


\subsection{Pseudo-Spectral Method with Semi-Explicit Scheme}
To describe the pseudo-spectral method, we consider the problem:
\ba
& & u_t = \epsilon\Delta u + \vec{B}(x,y,t)\cdot\nabla u + Cu \nonumber \\
& &  u(x,y,0)=1; u(0,y,t)=u(1,y,t), \\
&& u(x,0,t)=u(x,1,t).
\ea
The derivatives of solutions are $\nabla u = \mathcal{F}^{-1} (i\vec{k}\mathcal{F}(u))$ 
and $\triangle u = \mathcal{F}^{-1} ((i\vec{k})^2\mathcal{F}(u))$. Here $\mathcal{F}$ is the discrete Fourier transformation and $\vec{k}$ represents the wave numbers.

The pseudo-spectral scheme with semi-explicit scheme is that diffusion term uses implicit scheme and advection term uses half step lagged explicit scheme:
\begin{equation}\label{eq:semi}
\frac{\hat{U}^{n+1}-\hat{U}^{n}}{\Delta t} = -\epsilon\vec{k}^2 \hat{U}^{n+1} + \mathcal{F}(\vec{B}\cdot \mathcal{F}^{-1}(i\vec{k}\hat{U}^n)) + C \hat{U}^n
\end{equation}
%
Here the discretization on the spatial domain is $N\times N$ and $\hat{U}^{n}$ is the frequency of $u$ at time step $n$. 
Pseudo-spectral method has higher order accuracy than the finite element methods.  
However, its time step size $\Delta t$ depends on the small diffusion parameter $\epsilon$. Therefore it is more costly to reach the large time solution. 

\subsection{Stability Analysis}
The pseudo-spectral method is not so efficient 
when diffusion parameter $\epsilon$ goes to zero. Here 
we compare the upwinding finite element method and pseudo-spectral method. 
The analysis suggests a hybrid algorithm combining finite element method and pseudo-spectral method.

\subsubsection{Upwinding Finite element method}
The stability condition of EAFE method 
is equivalent to that of an upwinding scheme.  To simplify analysis, 
the advection term is taken as one-dimensional and the reaction term is ignored. 
The proto-type equation is
\begin{equation}\label{proto1}
u_t = \epsilon\, u_{xx} + b\, u_x
\end{equation}

The upwinding semi-explicit scheme is
\begin{equation}\label{eq:upwinding}
\left\{
\begin{array}{l}
\frac{U^{n+1}_j-U^n_j}{\Delta t} = \epsilon(\frac{U^n_{j+1}+U^n_{j-1} - 2U^n_j}{2} + \frac{U^{n+1}_{j+1}+U^{n+1}_{j-1} - 2U^{n+1}_j}{2})/\Delta x^2 + b\frac{U^n_{j+1}-U^n_j}{\Delta x}\quad\text{if}\  b\geq 0 \\
\frac{U^{n+1}_j-U^n_j}{\Delta t} = \epsilon(\frac{U^n_{j+1}+U^n_{j-1} - 2U^n_j}{2} + \frac{U^{n+1}_{j+1}+U^{n+1}_{j-1} - 2U^{n+1}_j}{2})/\Delta x^2 + b\frac{U^n_{j}-U^n_{j-1}}{\Delta x}\quad\text{if}\  b< 0
\end{array}
\right.
\end{equation}
Denoting $$ b^+=\max(b,0),\quad b^-=\min(b,0),$$ 
we carry out the von Neumann analysis to 
obtain the stability condition.
Upon substitution $U^n_j = \rho^n e^{ij\theta}$, (\ref{eq:upwinding}) becomes:
\begin{equation}\label{eq:upwinding2}
\frac{\rho-1}{\Delta t} = \frac{\epsilon(\rho+1)(\cos(\theta)-1)}{\Delta x^2} + \frac{|b|(\cos(\theta)-1+i\sin(\theta))}{\Delta x}
\end{equation}
Let $\mu = \frac{\epsilon\Delta t}{\Delta x^2}$ and 
CFL number $C = \frac{|b|\Delta t}{\Delta x}$. The 
Peclet number 
$$Pe = \frac{\text{advection}}{\text{diffusion}} = \frac{|b|\Delta x}{\epsilon} = C/\mu. $$
Equation (\ref{eq:upwinding2}) simplifies to:
\begin{equation}
\begin{array}{rl}
(1+\mu (1-\cos\theta))\rho =& (1-\mu(1-\cos\theta)) + Pe\, \mu(\cos\theta-1+i\sin\theta)\\
=& (1-\mu(1+\, Pe\, )(1-\cos\theta)) + i\, Pe\, \mu\sin\theta\\
\end{array}
\end{equation}

\noindent The stable condition is:
\begin{equation}
%\begin{array}{rl}
\rho \leq 1 
%\Longleftrightarrow
%&|1-\mu(1+Pe)(1-\cos\theta) + iPe\mu\sin\theta| \leq |1+\mu(1-\cos\theta)|\\
%\Longleftrightarrow
%&1-2\mu(1+Pe)(1-\cos\theta) + (\mu(1+Pe)(1-\cos\theta))^2 + (Pe\mu\sin\theta)^2\\
%&\leq 1+2\mu(1-\cos\theta) +(\mu(1-\cos\theta))^2 \\
%\Longleftrightarrow
%&((1+Pe)^2(1-\cos\theta)^2 + (Pe\sin\theta)^2)\mu\leq 2(1-\cos\theta)(2+Pe)\\
\Longleftrightarrow
 \mu \leq \frac{2(1-\cos\theta)(2+Pe) } {(1+Pe)^2(1-\cos\theta)^2 + (Pe\sin\theta)^2}
%\end{array}
\end{equation}
It suffices to impose:
\begin{equation*}
\begin{array}{rl}
\mu\leq &\min_{\theta\in[-\pi,\pi]}\frac{2(1-\cos\theta)(2+Pe) } {(1+Pe)^2(1-\cos\theta)^2 + (Pe\sin\theta)^2}\\
=&(2+Pe)/Pe^2\\
\rightarrow& 1/Pe \quad\text{while}\quad \epsilon\rightarrow 0
\end{array}
\end{equation*}
The stability condition of the upwinding scheme in the limit of $\epsilon \downarrow 0$ is:
\begin{equation}
\Delta t\leq\frac{\Delta x}{|b|}, 
\end{equation}
a CFL condition independent of $\epsilon$.

\subsubsection{Pseudo-Spectral Method}\
Consider pseudo-spectral method on the proto-type equation (\ref{proto1}) again:
we have shown before:
\begin{equation}\label{eq:semispectral}
\frac{U^{n+1}-U^n}{\triangle t} = -\epsilon \vec{k}^2 U^{n+1} + b \mathcal{F}^{-1}(i\vec{k}\mathcal{F}U^n).
\end{equation}

To prove the stability of pseudo-spectral method with Semi-Explicit scheme, we assume:
\begin{equation}\label{eq:basis}
u(x,t) = \sum^{\infty}_{k=-\infty} \hat{u}_k(t) e^{ikx}.
\end{equation}
and we denote the projection operator:
\begin{equation}\label{eq:projection}
(P_N)u(x,t) = \sum_{-\frac{N}{2}<k\leq\frac{N}{2}} \hat{u}_k(t) e^{ikx}.
\end{equation}

If we plug (\ref{eq:basis}) and (\ref{eq:projection}) into (\ref{eq:semispectral}), we have that the k-th component of the vector, which should be the corresponding coefficient of $e^{ikx}$:
\begin{equation}\label{eq:elementwise}
\frac{\hat{u}^{n+1}_k-\hat{u}^{n}_k}{\triangle t} = -\epsilon k^2 u^{n+1}_{k}  + ik b \hat{u}^n_k.
\end{equation}
From \eqref{eq:elementwise}, we obtain
\begin{equation}
|\hat{u}^{n+1}_k| = \left |\frac{1+ikb\triangle t}{1+\epsilon k^2\triangle t} \hat{u}^n_k \right | \leq \frac{|1+ikb\triangle t|}{|1+\epsilon k^2 \triangle t|} |\hat{u}^n_k|
\end{equation}
%
So $U^n$ in (\ref{eq:semispectral}) is stable if and only if:
\begin{equation}
\left |\frac{1+ikb\triangle t}{1+\epsilon k^2 \triangle t}\right |\leq 1 \quad \forall k \in \{ -\frac{N}{2}<k\leq \frac{N}{2}\}
\end{equation}
which is equivalent to
\begin{equation}
\begin{array}{rll}
%&|\frac{1+ikb\triangle t}{1+\epsilon k^2 \triangle t}|\leq 1 & \forall k \in \{ -\frac{N}{2}<k\leq \frac{N}{2}\} \\
%%\Longleftrightarrow
%%&|1+ikb\triangle t| \leq |1+\epsilon k^2 \triangle t| & \forall k \in \{ -\frac{N}{2}<k\leq \frac{N}{2}\}\\
%%\Longleftrightarrow
%%& 1+ k^2b^2\triangle t^2 \leq 1 + 2\epsilon k^2 \triangle t + \epsilon^2k^4\triangle t^2 & \forall k \in \{ -\frac{N}{2}<k\leq \frac{N}{2}\}\\
%%\Longleftrightarrow
%%&b^2\triangle t \leq 2\epsilon + \epsilon^2 k^2 \triangle t & \forall k \in \{ -\frac{N}{2}<k\leq \frac{N}{2}, k\neq 0\}\\
%\Longleftrightarrow 
& \triangle t \leq \frac{2\epsilon}{b^2-\epsilon^2 k^2}& \forall k \in \{ -\frac{N}{2}<k\leq \frac{N}{2}, k\neq 0, b^2-\epsilon^2 k^2 >0 \}
\end{array}
\end{equation}
So if the equation is advection dominant, then at least $k=1$ satisfies  
\begin{equation}
-\frac{N}{2}<k\leq \frac{N}{2},\quad k\neq 0,\quad b^2-\epsilon^2 k^2 >0
\end{equation}
and the stable condition would be $\triangle t\leq  \frac{2\epsilon}{b^2-\epsilon^2}$. When $ \epsilon \ll 1 $, the time step is in the order of $\mathcal O(\epsilon)$ which is very restrictive.


%\subsubsection{Pseudo-Spectral Method with nonconstant advection}\
%Consider pseudo-spectral method on advection-diffusion equation 
%\begin{equation}\label{proto2}
%u_t = \epsilon\, u_{xx} + B(x)\, u_x
%\end{equation}
%we have shown before:
%\begin{equation}\label{eq:semispectral2}
%\frac{U^{n+1}-U^n}{\triangle t} = -\epsilon \vec{k}^2 U^{n+1} + B(x)\cdot \mathcal{F}^{-1}(i\vec{k}\mathcal{F}U^n)
%\end{equation}

%To prove the stability of pseudo-spectral method with Semi-Explicit scheme, we assume:
%\ba
%u(x,t) &=& \sum^{\infty}_{k=-\infty} \hat{u}_k(t) e^{ikx}\\
%B(x) &=& \sum^{\infty}_{k=-\infty} \hat{B}_k e^{ikx}\label{eq:basis2}
%\ea
%and we denote the projection operation:
%\begin{equation}\label{eq:projection2}
%(P_N)u(x,t) = \sum_{-\frac{N}{2}<k\leq\frac{N}{2}} \hat{u}_k(t) e^{ikx}
%\end{equation}

%If we plug (\ref{eq:basis2}) and (\ref{eq:projection2}) into (\ref{eq:semispectral2}), we have that the i-th component of the vector, which should be the corresponding coefficient of $e^{ikx}$:
%\begin{equation}\label{eq:elementwise2}
%\frac{\hat{u}^{n+1}_k-\hat{u}^{n}_k}{\triangle t} = -\epsilon k^2 u^{n}_k + \sum_{
%\begin{array}{l}
%s+t=k,\forall s\\
%-\frac{N}{2}<t	\leq\frac{N}{2}
%\end{array}
%} it \hat{B}_s \hat{u}^n_t
%\end{equation}

%and the matrix form of (\ref{eq:semispectral2}) will be:
%\begin{equation}\label{eq:MatrixForm2}
%\hat{U}^{n+1} = D^{-1}\,(I + \triangle t\,\hat{B}\, T)\,\hat{U}^n
%\end{equation}
%where 

%\begin{equation}
%\hat{U}^n = (\hat{u}_{-\frac{N}{2}+1} \cdots \hat{u}_{\frac{N}{2}})^T
%\end{equation}

%\begin{equation}
%D = \text{diag}(1+\epsilon\triangle t(\frac{N-1}{2})^2, 1+\epsilon\triangle t(\frac{N-2}{2})^2, \cdots 1 \cdots, 1+\epsilon\triangle t(\frac{N}{2})^2) 
%\end{equation}

%\begin{equation}
%T = \text{diag}(-\frac{N}{2}+1, -\frac{N}{2}+2, \cdots, \frac{N}{2})
%\end{equation}

%and
%\begin{equation}
%\hat{B} = \left[
%\begin{array}{llll}
%\hat{B}_0, & \hat{B}_{-1}, & \cdots & \hat{B}_{-N+1} \\
%\hat{B}_{-1}, & \hat{B}_{-2}, & \cdots & \hat{B}_{-N+2} \\
%\cdots & \cdots & & \\
%\hat{B}_{N-1}, & \hat{B}_{N-2}, & \cdots & \hat{B}_0
%\end{array}
%\right]
%\end{equation}

%Consider the infinite norm of vector:
%\begin{equation}
%||\vec{U}||_{\infty} = \max_{k} |U_k|
%\end{equation}
%and the infinite norm of matrix:
%\begin{equation}
%||\hat{B}||_{\infty} = \max_{i}\sum_{j} |\hat{B}_{ij}|
%\end{equation}
%So from (\ref{eq:elementwise}) we have 

%\begin{equation}
%\begin{array}{ll}
%|\hat{u}^{n+1}_k| &\leq \frac{1}{1+\epsilon k^2\triangle t}|\hat{u}^n_k| + \frac{\frac{N}{2}\triangle t||\hat{B}||_{\infty} ||\hat{U}^n||_{\infty}}{1+\epsilon k^2\triangle t}\\
% & \leq \frac{1+\frac{N}{2}\triangle t ||\hat{B}||_{\infty}}{1+\epsilon k^2\triangle t} ||\hat{U}^n||_{\infty} 
%\end{array}
%\end{equation}

%{\bf 
%Last Friday, we made a mistake about the numerator and thought it's $k$ instead of $N$. If the estimator is  $\frac{1+\frac{k}{2}\triangle t ||\hat{B}||_{\infty}}{1+\epsilon k^2\triangle t}$, then the minimal point is when $k=[\frac{\sqrt{\epsilon^2-||\hat{B}||^2_{\infty}\epsilon\triangle t/4} - ||\hat{B}||_{\infty}\triangle t/2}{\epsilon ||\hat{B}||_{\infty} \triangle t}]$.

%However if the coefficient is $\frac{1+\frac{k}{2}\triangle t ||\hat{B}||_{\infty}}{1+\epsilon k^2\triangle t}$ as it supposed to be, the minimal value of $\frac{1+\frac{k}{2}\triangle t ||\hat{B}||_{\infty}}{1+\epsilon k^2\triangle t}$ is when $k=\frac{N}{2}$ and we can't get the stability condition because $\triangle 5$ will be canceled.
%}
%*******************************************************************************************************
% Section 5
%*******************************************************************************************************
\section{Numerical Results on Residual Speeds}
For each given $\lambda$, we compute the Lyapunov exponent of~\eqref{varp1} by finding the long time 
quasi-stationary solution $u$ and applying formula \eqref{lyapunov}. A quasi-stationary 
solution with $\lambda=1$ is shown in the left plot of Fig.~\ref{w_hat}. The layer structure of the numerical solution 
moves along the $y=x$ direction and oscillates between the two neighboring invariant manifolds (lines of unit slope). 
The right plot of Fig.~\ref{w_hat} shows the power spectrum of the auxiliary field (a genearalized eigenfunction) 
$w$ (right) for the unsteady flow~\eqref{unsteady}. We observe that the energy decreases 
towards high frequency, with a linear decay (scaling range) in the intermediate region of wave numbers.
%
\begin{figure}
\includegraphics[width=0.5\textwidth, height = 0.4\textwidth]{w.eps}
\includegraphics[width=0.5\textwidth, height = 0.4\textwidth]{w_hat.eps}
\caption{A quasi-stationary solution (left) at time $t=2$, $\epsilon=10^{-3}$, $|B|=500$, and $\lambda=1$. 
The number of space nodes $N=2^9$ in each direction and time step $dt=1/10^8$. The power spectrum of the 
auxiliary field $w$ (right) for the unsteady flow~\eqref{unsteady} shows that 
the energy decreases towards high frequency, with a scaling range (slope = -1.1629) 
in the intermediate region of wave numbers.}\label{w_hat}
\end{figure}

After calculating $\mu$ by the algorithm discussed in (\ref{lyapunov}), the minimization problem (\ref{eq:inf}) is solved 
by the golden-section method\cite{B73}. Golden-section method is a method to find the minimal point of a particular function called unimodal function. The idea of Golden-section method is by successively narrowing down the range of searching interval. A common definition of unimodal function is for some value $m$, it is monotonically decreasing when $x<m$ and monotonically increasing when $x\geq m$. 
%The left plot of Figure (\ref{fig:minimize}) 
In~\cite{NX05}, the authors show that the function $\frac{\mu(\lambda)}{\lambda}$ is indeed a unimodal function which is  
strictly convex with respect to $\lambda$, 
decreasing over interval $[0, \lambda^*]$ and strictly increasing 
over interval $[\lambda^*, \infty]$, with $\lambda^*$ the unique
minimal point of $\frac{\mu(\lambda)}{\lambda}$. Notice that Newton's method is 
not applicable since $\mu'(\lambda)$ is not known.

%According to figure (\ref{fig:unsteady}), we can see the 
%damped oscillation of $\text{ln}\int_0^L u(x,y,T)dy / T$ ( $T=?$). In order to 
%compute the result faster, the limit can be approximated 
%by the average between $[2\pi, 4\pi]$ and it's a reasonable approximation 
%based on figure (\ref{fig:unsteady}).\\
%
%\begin{figure}[ht!]
%\includegraphics[width=1\textwidth]{figure/eps/unsteady.eps}
%\caption{Relation between time t and $\int_0^T u(x,y,T)dy /T$. $\epsilon = 1e-5$ and $\lambda$ is selected as 21.980675, corresponding to the minimal point of $\mu(\lambda)/\lambda$. }
%\label{fig:unsteady}	
%\end{figure}

For time periodic cellular flow (\ref{tdcell}), EAFE method and spectral method are combined together to obtain the result in Fig.~\ref{residual}. The initial searching interval is $[0, 1000]$. EAFE is used first to narrow down the searching interval to an neighborhood of $\lambda^*$ with length 2. The space size and time step is $h = 1/2^{10}$ and $dt=1/10^7$, respectively. Since the minimal diffusion parameter $\epsilon=\mathcal{O}(10^{-3})$ and amplitude of advection term is $\mathcal{O}(10^{3})$, the choice of $dt = 1/10^7$ leads to numerical method stable. The spectral method with space node $N = 2^9$ in each direction and time step $dt=1/10^8$ is used to obtain the minimal point when the length of searching interval is less than $2$. The stopping criterion of searching algorithm is that the difference of two consecutive value of $\frac{\mu(\lambda)}{\lambda}$ is less than $5\times10^{-3}$. The benefit of using two methods together is that we can find the minimal point faster without the sacrifice  of the accuracy. 

In Fig.~\ref{residual}, we plot $c^*$ as a function of $\epsilon$ for 
three values of parameter $\theta =10^{-2},10^{-1},1$. Recall that $\theta$ is the perturbation parameter which determines how close the flow is to the steady cellular flow; see \eqref{tdcell}. 
From Fig.~\ref{residual}, we observe that for all three positive values of $\theta$, $c^*(\epsilon) = O(1)$ 
as $\epsilon \downarrow 0$. In contrast, it is well known that $c^{*}(\epsilon) = O(\epsilon^{1/4})$ in 
steady cellular flows\cite{NR05}. The presence of chaotic 
trajectories in time periodic cellular flows contributed to this phenomenon. 
They are much more mobile and far reaching than their counterparts of steady cellular flows.
Fig.~\ref{residualGrid} shows the persistence of residual front speeds 
under grid refinement at $\theta = 1$ and the close proximity of the data points 
between those at grid sizes of  $N = 2^8$ and $N = 2^{9}$ validates the convergence of the results numerically. 

\begin{figure}[ht!]
\centering
\includegraphics[width=1\textwidth, height = 0.8\textwidth]{residual.eps}
\caption{Residual KPP front speeds in the small diffusion limit $\epsilon \downarrow 0$.}
\label{residual}	
\end{figure}

\begin{figure}[ht!]
\centering
\includegraphics[width=1\textwidth, height = 0.8\textwidth]{h_convergence.eps}
\caption{Residual KPP front speeds at $\theta =1$ under grid refinement.}
\label{residualGrid}	
\end{figure}

%Then we can see from the table that the minimal speed $c^*$ almost stays the same.\\

%\begin{tabular}{|c|c|c|c|c|c|c|}%\label{tab:1}
%  \hline
%  $\epsilon$ & 1e-5& 5e-5& 1e-4& 5e-4& 1e-3& 5e-3 \\
%  \hline
%  $c^*$ & 0.912845 & 0.913268 & 0.914533 & 0.922101 & 0.930152 &  0.975405\\
%  \hline
%\end{tabular}
%\label{tab:1}\\


%\section{time step of different scheme}
%The following table shows numerical restriction of time while $\epsilon= 1e-2$ and $Amp = 500$.\\

%\begin{tabular}{|p{3cm}|c|p{3cm}|c|c|}%\label{tab:1}
%  \hline
%  Scheme & FEM  & Semi-implicit & Shrinkage & High order splitting \\
%  \hline
%  order h & h  & $Ch^n$ (n depends on smoothness) & $Ch^n$ & $Ch^n$ \\
%  \hline
%  order $\triangle t$ & $\triangle t$ &$\triangle t$ &$\triangle t$& $\triangle t^2$ \\
%  \hline
%  Theoretic dt & unconditional& $\frac{2\epsilon}{||B||_{\infty}}$  &  $\frac{2\epsilon}{||B||_{\infty}}$ & $C\frac{2\epsilon}{||B||_{\infty}}$\\
%  \hline
%  Numerical dt & 1e-4 & 1e-7 & 1e-7 & 1e-7\\
%  \hline
%  shrinkage parameter & & & $(\triangle t)^{1+0.3}$ & $(\triangle t)^{2+0.01}$\\
%  \hline
%  1/h   & $2^{10}$ & $2^9$ &$2^9$&$2^9$\\
%  \hline
%  CPU time(s) & 42.375980 & 2251.478948 & 1933.723103 & 4958.21018\\
%  \hline
%\end{tabular}
%\label{tab:1}\\

\section{Concluding Remarks}
We have studied the KPP front speed asymptotics computationally 
in time periodic cellular flows with chaotic streamlines in the 
small molecular diffusivity limit. The chaotic streamlines statistically resemble a sub-diffusion process and 
enhance the KPP front speed $c^*$ significantly in the sense that $c^*$ has a positive limit 
as molecular diffusivity tends to zero. Such residual transport phenomenon is absent in steady 
cellular flows with ordered streamlines. To facilitate effective computation in the advection dominated regime, 
we combined an upwinding finite element methods with the spectral method for computing the principal eigenvalue 
of a time periodic parabolic operator with small diffusion and for a subsequent minimization. 
In future work, we plan to study residual KPP front speeds in three space dimensional chaotic flows.  


\section{Acknowledgements}
The work was partially supported by NSF grant DMS-1211179. We thank Tyler McMillen for 
helpful conversations and visualization of chaotic dynamical systems. L.C 
was also supported by NSF grant DMS-1115961 and DOE prime award \# DE-SC0006903.

\newpage 
\begin{thebibliography}{9999999}
\bibliographystyle{plain}
\bibitem{A} E. Abraham, et al., {\em Importance of stirring 
in the development of an ironfertilized phytoplanton bloom}, Nature 407 (2000), pp. 727--730.

\bibitem{ACVV02} M. Abel, M. Cencini, D. Vergni, A. Vulpiani,
{\em Front speed enhancement in cellular flows}, Chaos 12 (2002), pp. 481--488.

\bibitem{Aud}B. Audoly, H. Berestycki, Y. Pomeau,
{\em R\'eaction diffusion en \'ecoulement
stationnaire rapide}, C.R. Acad. Sci. Paris, S\'erie IIb, 328 (2000), pp. 255--262.

\bibitem{B73} R. Brent, ``Algorithms for Minimization Without Derivatives'', Prentice Hall,
1973, DFMIN module http://gams.nist.gov/serve.cgi/Module/NMS/DFMIN/5671/.

\bibitem{BBFS} R. Bank, J. Burger, W. Fichtner, R. Smith, 
{\em Some upwinding techniques for finite element approximations of convection diffusion equations}, 
Numer. Math., 58 (1990), pp. 185--202.

\bibitem{BCVV95} L. Biferale, A. Cristini, M. Vergassola, A. Vulpiani, 
{\em Eddy diffusivities in scalar transport}, Physics Fluids 7(11), 1995, pp. 2725--2734.

\bibitem{BH82} A. Brooks, T. Hughes, {\em Streamline upwind/Petrov-Galerkin formulations for 
convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations.} 
Computer methods in applied mechanics and engineering 32(1), 1982, pp. 199-259.

\bibitem{CW91}R. Camassa, S. Wiggins, {\em Chaotic advection in a Rayleigh-B\'enard flow,} 
Physical Review A, 43(2), 1990, pp. 774--797.

\bibitem{ChSo}S. Childress, A. Soward,
{\em Scalar transport and alpha-effect for a family of cat's eye
flows}, J. Fluid Mech. 205 (1989), pp. 99--133.

\bibitem{CB83} J. Cooper, J. Butcher, {\em An Iterative Scheme for Implicit Runge-Kutta Methods,}
IMA J. Numer. Anal. 3, pp. 127--140, 1983.

\bibitem{ESW} H. Elman, D. Silvester, A. Wathen, 
``Finite Elements and Fast Iterative Solvers with Applications in Incompressible Fluid Dynamics'',
Oxford University Press, 2005.

\bibitem{FP} A. Fannjiang, G. Papanicolaou, 
{\em Convection enhanced diffusion for periodic flows,} SIAM J. Applied Math, 54(2), 1994, pp. 333-408.

\bibitem{KXZ03}  H. Kim, J. Xu, and L. Zikatanov, 
{\em A multigrid method based on graph matching for convection-diffusion equations,} Numer. Linear Algebra 
Appl., (2003), pp. 181 -- 195. More details are at: http://epubs.siam.org/doi/ref/10.1137/060659545

\bibitem{Nov11}Y. Gorb, D. Nam, and A. Novikov,  
{\em Numerical simulations of diffusion in cellular flows
at high Peclet numbers}, DCDS-B, vol. 51(1), 2011, pp. 75--92.

\bibitem{H} S. Heinze, {\em Diffusion-advection in cellular flows with large Peclet numbers,}
 Archive Rational Mech. Analysis, 168(4), 2003, pp. 329--342.

\bibitem{Liu_13} Y. Liu, J. Xin, Y. Yu, 
{\em A Numerical Study of Turbulent Flame Speeds of Curvature and Strain G-equations in Cellular Flows}, 
Physica D, 243(1), pp. 20-31, 2013.  

\bibitem{NX09} J. Nolen, J. Xin, {\em Asymptotic Spreading of KPP Reactive Fronts 
in Incompressible Space-Time Random Flows}, Ann Inst. H. Poincar\'e, Analyse Non Lineaire, 26(3), 2009, pp. 815-839.

\bibitem{NX05} J. Nolen, M. Rudd, J. Xin, {\em Existence of KPP fronts 
in spatially-temporally periodic advection and variational principle for propagation speeds}, 
 Dynamics of PDE, 2(1), pp. 1-24, 2005.

\bibitem{N} Y. Notay, AGMG software and documentation; see http://homepages.ulb.ac.be/~ynotay/AGMG

\bibitem{NR05}A. Novikov, L. Ryzhik,
{\em  Boundary layers and KPP fronts in a cellular flow}, Arch. Ration. Mech. Anal., 
184(1), 2007, pp. 23 -- 48.

\bibitem{PS05} M. Paoletti, T. Solomon, {\em Experimental studies of front propagation and mode-locking 
in an advection-reaction-diffusion system}, Europhys. Lett. 69 (2005), pp. 819--825.

\bibitem{R}P. Ronney, {\em Some Open Issues in Premixed Turbulent Combustion},
Modeling in Combustion Science (J. D. Buckmaster and T. Takeno, Eds.),
Lecture Notes In Physics, Vol. 449, Springer-Verlag, Berlin, 1995, pp. 3-22.

\bibitem{RZ07} L. Ryzhik, A. Zlatos, {\em KPP pulsating front speed-up by flows}, Commun. Math. Sci. 5 (2007), 
pp. 575--593.

\bibitem{S03}W. van Saarloos, 
{\em Front propagation into unstable states}, Physics
Reports 386(2), 2003, pp. 29-222.

%\bibitem{SCHO} H. Schaeffer, R. Caflisch, C. Hauck, S. Osher, 
%{\em Sparse dynamics for partial differential equations}, Proc. National Acad. Sci, 
%March, 2013, doi:10.1073/pnas.1302752110

\bibitem{SZX13} L. Shen, A. Zhou, J. Xin, {\em Finite Element Computation of KPP Front Speeds in Cellular and 
Cat's Eye Flows}, Journal of Scientific Computing, 55(2), 2013, pp. 455--470.

\bibitem{TPB71} P. Tong, T.H.H. Pian, L.L. Bucciarelli, {\em Mode shapes and frequencies by finite element method using consistent and lumped masses}, Comput. Struct.
1 (1971), pp. 623-638.

\bibitem{W85}F. Williams, ``Combustion Theory'', Benjamin-Cummings, Menlo Park, 1985.

\bibitem{X09}J. Xin, ``An Introduction to Fronts in Random Media'', Surveys and Tutorials in Applied Math Sciences, 
Vol. 5, Springer, 2009. 

\bibitem{Yu1}J. Xin, Y. Yu,
{\em Front Quenching in G-equation Model Induced by Straining of Cellular Flow}, 
Arch. Rational Mechanics Analysis, to appear.  

\bibitem{Yu2} J. Xin, Y. Yu,
{\em Asymptotic growth rates and strong bending of turbulent flame speeds of 
G-equation in steady two dimensional incompressible periodic flows}, SIAM J. Math Analysis, to appear.  

\bibitem{XZ99} J. Xu, L. Zikatanov, {\em A monotone finite element scheme 
for convection diffusion equations}, Math. Comp., 68 (1999), pp. 1429--1446.


%\bibitem{wiki} $http://en.wikipedia.org/wiki/Golden\_section\_search$.
%
\end{thebibliography}

\end{document}
	

1406.2268.extracted_bib
{A E. Abraham, et al., {\em Importance of stirring in the development of an ironfertilized phytoplanton bloom, Nature 407 (2000), pp. 727--730.
1406.2268.extracted_bib
{ACVV02 M. Abel, M. Cencini, D. Vergni, A. Vulpiani, {\em Front speed enhancement in cellular flows, Chaos 12 (2002), pp. 481--488.
1406.2268.extracted_bib
{AudB. Audoly, H. Berestycki, Y. Pomeau, {\em R\'eaction diffusion en \'ecoulement stationnaire rapide, C.R. Acad. Sci. Paris, S\'erie IIb, 328 (2000), pp. 255--262.
1406.2268.extracted_bib
{B73 R. Brent, ``Algorithms for Minimization Without Derivatives'', Prentice Hall, 1973, DFMIN module http://gams.nist.gov/serve.cgi/Module/NMS/DFMIN/5671/.
1406.2268.extracted_bib
{BBFS R. Bank, J. Burger, W. Fichtner, R. Smith, {\em Some upwinding techniques for finite element approximations of convection diffusion equations, Numer. Math., 58 (1990), pp. 185--202.
1406.2268.extracted_bib
{BCVV95 L. Biferale, A. Cristini, M. Vergassola, A. Vulpiani, {\em Eddy diffusivities in scalar transport, Physics Fluids 7(11), 1995, pp. 2725--2734.
1406.2268.extracted_bib
{BH82 A. Brooks, T. Hughes, {\em Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Computer methods in applied mechanics and engineering 32(1), 1982, pp. 199-259.
1406.2268.extracted_bib
{CW91R. Camassa, S. Wiggins, {\em Chaotic advection in a Rayleigh-B\'enard flow, Physical Review A, 43(2), 1990, pp. 774--797.
1406.2268.extracted_bib
{ChSoS. Childress, A. Soward, {\em Scalar transport and alpha-effect for a family of cat's eye flows, J. Fluid Mech. 205 (1989), pp. 99--133.
1406.2268.extracted_bib
{CB83 J. Cooper, J. Butcher, {\em An Iterative Scheme for Implicit Runge-Kutta Methods, IMA J. Numer. Anal. 3, pp. 127--140, 1983.
1406.2268.extracted_bib
{ESW H. Elman, D. Silvester, A. Wathen, ``Finite Elements and Fast Iterative Solvers with Applications in Incompressible Fluid Dynamics'', Oxford University Press, 2005.
1406.2268.extracted_bib
{FP A. Fannjiang, G. Papanicolaou, {\em Convection enhanced diffusion for periodic flows, SIAM J. Applied Math, 54(2), 1994, pp. 333-408.
1406.2268.extracted_bib
{KXZ03 H. Kim, J. Xu, and L. Zikatanov, {\em A multigrid method based on graph matching for convection-diffusion equations, Numer. Linear Algebra Appl., (2003), pp. 181 -- 195. More details are at: http://epubs.siam.org/doi/ref/10.1137/060659545
1406.2268.extracted_bib
{Nov11Y. Gorb, D. Nam, and A. Novikov, {\em Numerical simulations of diffusion in cellular flows at high Peclet numbers, DCDS-B, vol. 51(1), 2011, pp. 75--92.
1406.2268.extracted_bib
{H S. Heinze, {\em Diffusion-advection in cellular flows with large Peclet numbers, Archive Rational Mech. Analysis, 168(4), 2003, pp. 329--342.
1406.2268.extracted_bib
{Liu_13 Y. Liu, J. Xin, Y. Yu, {\em A Numerical Study of Turbulent Flame Speeds of Curvature and Strain G-equations in Cellular Flows, Physica D, 243(1), pp. 20-31, 2013.
1406.2268.extracted_bib
{NX09 J. Nolen, J. Xin, {\em Asymptotic Spreading of KPP Reactive Fronts in Incompressible Space-Time Random Flows, Ann Inst. H. Poincar\'e, Analyse Non Lineaire, 26(3), 2009, pp. 815-839.
1406.2268.extracted_bib
{NX05 J. Nolen, M. Rudd, J. Xin, {\em Existence of KPP fronts in spatially-temporally periodic advection and variational principle for propagation speeds, Dynamics of PDE, 2(1), pp. 1-24, 2005.
1406.2268.extracted_bib
{N Y. Notay, AGMG software and documentation; see http://homepages.ulb.ac.be/~ynotay/AGMG
1406.2268.extracted_bib
{NR05A. Novikov, L. Ryzhik, {\em Boundary layers and KPP fronts in a cellular flow, Arch. Ration. Mech. Anal., 184(1), 2007, pp. 23 -- 48.
1406.2268.extracted_bib
{PS05 M. Paoletti, T. Solomon, {\em Experimental studies of front propagation and mode-locking in an advection-reaction-diffusion system, Europhys. Lett. 69 (2005), pp. 819--825.
1406.2268.extracted_bib
{RP. Ronney, {\em Some Open Issues in Premixed Turbulent Combustion, Modeling in Combustion Science (J. D. Buckmaster and T. Takeno, Eds.), Lecture Notes In Physics, Vol. 449, Springer-Verlag, Berlin, 1995, pp. 3-22.
1406.2268.extracted_bib
{RZ07 L. Ryzhik, A. Zlatos, {\em KPP pulsating front speed-up by flows, Commun. Math. Sci. 5 (2007), pp. 575--593.
1406.2268.extracted_bib
{S03W. van Saarloos, {\em Front propagation into unstable states, Physics Reports 386(2), 2003, pp. 29-222. %
1406.2268.extracted_bib
{SCHO H. Schaeffer, R. Caflisch, C. Hauck, S. Osher, %{\em Sparse dynamics for partial differential equations, Proc. National Acad. Sci, %March, 2013, doi:10.1073/pnas.1302752110
1406.2268.extracted_bib
{SZX13 L. Shen, A. Zhou, J. Xin, {\em Finite Element Computation of KPP Front Speeds in Cellular and Cat's Eye Flows, Journal of Scientific Computing, 55(2), 2013, pp. 455--470.
1406.2268.extracted_bib
{TPB71 P. Tong, T.H.H. Pian, L.L. Bucciarelli, {\em Mode shapes and frequencies by finite element method using consistent and lumped masses, Comput. Struct. 1 (1971), pp. 623-638.
1406.2268.extracted_bib
{W85F. Williams, ``Combustion Theory'', Benjamin-Cummings, Menlo Park, 1985.
1406.2268.extracted_bib
{X09J. Xin, ``An Introduction to Fronts in Random Media'', Surveys and Tutorials in Applied Math Sciences, Vol. 5, Springer, 2009.
1406.2268.extracted_bib
{Yu1J. Xin, Y. Yu, {\em Front Quenching in G-equation Model Induced by Straining of Cellular Flow, Arch. Rational Mechanics Analysis, to appear.
1406.2268.extracted_bib
{Yu2 J. Xin, Y. Yu, {\em Asymptotic growth rates and strong bending of turbulent flame speeds of G-equation in steady two dimensional incompressible periodic flows, SIAM J. Math Analysis, to appear.
1406.2268.extracted_bib
{XZ99 J. Xu, L. Zikatanov, {\em A monotone finite element scheme for convection diffusion equations, Math. Comp., 68 (1999), pp. 1429--1446. %
1406.2268.extracted_bib
{wiki $http://en.wikipedia.org/wiki/Golden\_section\_search$. %