\section{Introduction}

Efficient markets hypothesis (EMH) has been a cornerstone of financial economics for decades and it has been brought to the centre by the influential paper of \cite{Fama1970} summarizing empirical findings following the idea of the efficient markets hypothesis by \cite{Fama1965} and \cite{Samuelson1965}. Even though the actual definitions differ, the former study builds on a random walk definition and the latter one on a martingale definition, the qualitative consequences are the same -- the efficiency of a market originates in impossibility of systematic beating of the market, usually in a form of above-average risk-adjusted returns. \cite{Fama1991} later separated the efficiency hypothesis into three forms -- weak, medium and strong -- which differ by different information sets taken into consideration and all are based on inclusion of the information sets in market prices. The weak-form EMH says that all past price movements (and associated statistics) are already reflected in the market prices. Prediction of market movements based on historical time series (technical analysis) is thus not possible for this form. The medium-form EMH states that all publicly available information are already contained in the prices, the strong-form EMH adds all (even privately) available information. The medium-form thus discards fundamental analysis as well and the strong-form eliminates even insiders from making profit. Evidently, a weaker form of EMH is always a subset of a stronger form. Even though EMH has been repeatedly disparaged both empirically and theoretically \citep{Cont2001,Malkiel2003}, and even more so after an outbreak of the Global Financial Crisis in 2007/2008, its validity remains an open issue, yet still it persists in standard textbooks of financial economics \citep{Elton2003}.

Comparison of efficiency across various assets has been discussed in different studies. In a series of papers, \cite{DiMatteo2003,DiMatteo2005} and \cite{DiMatteo2007} study long-term memory and multi-scaling of a wide portfolio of stock indices, foreign exchange rates, Treasury rates and Eurodollar interbank interest rates using various estimators of long-term memory. They show that stock indices of more developed countries are also more efficient yet showing a weak signs of anti-persistence (properties of long-term memory are described in detail in the Methodology section), finding no deviations from EMH for all analyzed maturities of Eurodollar and Treasury rates. For US dollar exchange rates, the authors find diverse results with no evident pattern connecting exchange rate efficiency level with geographical or geopolitical properties. In another series of papers, \cite{Cajueiro2004b,Cajueiro2004,Cajueiro2004a,Cajueiro2005} compare stock market indices from different continents finding that the US and Japanese markets are the most efficient ones whereas the Asian and Latin American ones are detected as the least efficient ones. \cite{Lim2007} studies non-linear dependencies, their evolution in time and connection to market efficiency for a set of stock markets. The author finds the US market to be the most efficient one followed by Korea, Taiwan and Japan. On the other end of the ranking, Malaysia, Chile and Argentina are placed. \cite{Zunino2010} utilize the complexity-entropy causality plane to rank stock market indices to show that the emergent markets are less efficient than the developed ones as one would expect. The difference is attributed to a lower entropy and a higher complexity of the emergent markets. \cite{Kristoufek2013} introduce the Efficiency Index and come up with a ranking of stock market indices finding that the most efficient markets are located in Western Europe, USA and Japan whereas the least efficient markets are situated in Latin America and Asia.

However, up to our best knowledge, proper attention has not been given to a comparison of the efficiency of commodity markets. In this paper, we analyze futures markets for a wide range of commodities -- energy, metals and various agricultural commodities -- and compare their efficiency using the Efficiency Index proposed by \cite{Kristoufek2013}. The paper is structured as follows. Section 2 covers literature dealing with the efficiency of commodities. Section 3 describes the methodology in detail. Section 4 describes the analyzed dataset and brings the results. Section 5 concludes. We show that efficiency is related to a type of commodity (energy commodities being the most efficient ones and other agricultural commodities being the least efficient ones). In addition, we find a nonstandard relationship between the local and global properties of the series as most of the series show local persistence while in global, they are mean-reverting. The series thus follow quite strong local trends but in a long term, they return to their fundamental value.

\section{Literature review}

Testing the market efficiency in commodities markets has a long history. \cite{Roll1972} examines the commodity price index and argues that the market is inefficient due to significant serial correlations of its returns. \cite{Danthine1977} disputes such claim and shows that the violation of the standard martingale condition does not imply inefficiency in the commodity spot markets with support of risk aversion and no arbitrage opportunities. \cite{Gjolberg1985} analyzes oil spot prices at the Rotterdam market, rejects the efficiency hypothesis and constructs a profitable trading rule for daily, weekly and monthly price changes. \cite{Panas1991} studies the Rotterdam oil market as well together with the Italian market and based on leptokurtic monthly price changes, he rejects the markets' efficiency. \cite{Herbert1996} examine the US spot (cash) and futures markets for natural gas and find these to be inefficient. They argue that such inefficiency is caused by the specific structure of the US gas markets.

More recently, \cite{Tabak2007} analyze the efficiency of Brent and WTI crude oil using the rescaled range analysis and show that the markets are becoming more efficient in time. \cite{Alvarez-Ramirez2008} study the auto-correlation structure of the crude oil process using the detrended fluctuation analysis. They show that in long-term, the market is efficient but in short-term, the auto-correlation structure leads to rejection of the efficiency. \cite{Alvarez-Ramirez2010} further inspect the crude oil markets using lagged detrended fluctuation analysis and argue that multi-scaling and deviations from the random walk behavior cause the spot prices to be inefficient. The research on evolution of efficiency in time is further extended by \cite{Wang2010} where the authors study short-, medium- and long-term efficiency for various scales of the detrended fluctuation analysis approach. They show that the WTI crude oil becomes more efficient in time for all three analyzed scales. Using also the detrended fluctuation analysis, \cite{Wang2011} argue that WTI crude oil spot and futures are not efficient for short time scales below one month. Crude oil markets (Brent and WTI) are also analyzed by \cite{Charles2009} who use the variance ratio tests to show that the Brent market is weak-form efficient but the WTI market is not while providing some discussion about effects of deregulation on the markets.

\cite{Lee2009} study four energy commodities -- coal, oil, gas and electricity -- using panel data stationarity tests to uncover that none of the studied markets is efficient in the strict stationarity sense. \cite{Lean2010} study WTI crude oil spot and futures prices using mean-variance and stochastic dominance approaches finding no arbitrage opportunities between spot and futures prices while the findings are robust for various sub-periods and critical events. \cite{Narayan2010} study long-term relationship between spot and futures prices of gold and oil. They find that investors use the gold market to hedge against inflation and for our purposes also more importantly that crude oil market predicts the gold market and vice versa implying inefficiency.

\cite{Wang2010a} study high-frequency futures data of crude oil, heating oil, gasoline and natural gas using various nonlinear models. For heating oil and natural gas, the authors find market inefficiencies which are profound mainly during the bull market conditions. \cite{Gebre-Mariam2011} focuses on the US natural gas market (spot and futures) finding no arbitrage opportunities for daily prices but in general, the author claims that the markets can be seen as efficient only for contracts with approximately a month to maturity. \cite{Martina2011} utilize entropy approaches to WTI crude oil spot prices and find various cycles in its prices. Entropy is also applied by \cite{Ortiz-Cruz2012} who again study daily WTI prices finding the market to be efficient with two episodes of inefficiency connected to the early 1990s and late 2000s US recessions. The authors stress that deregulation of the market has helped improving its efficiency.

\cite{Zunino2011} apply information theory methods (specifically the permutation entropy and permutation statistical complexity) to the commodity markets allowing them for efficiency ranking finding silver, copper and cotton to be the most efficient commodities. \cite{Wang2011} study the gold market using the multifractal detrended fluctuation analysis to show that the market becomes more efficient in time especially after 2001. \cite{Kim2011} use the random matrix theory and network analysis to show that stock and commodity markets are well decoupled except for oil and gold showing signs of inefficiency. \cite{Kim2011a} then focus on the Korean agricultural market using the detrended fluctuation analysis finding anti-correlated series with strong volatility clustering aiming at inefficiency.

Out of these selected papers, it is evident that analysis of efficiency of commodity markets is fruitful with many approaches to the topic. However, the studies usually focus on a single (or a pair) of efficiency measures to test whether the specific markets are or are not efficient. Moreover, the analysis is usually strongly focused on a single commodity or a small group of commodities. Here, we contribute to the literature by applying various efficiency measures on a wide portfolio of commodities ranging from energy and agricultural (with several subgroups) commodities to metals. Moreover, we utilize the efficiency measure introduced by \cite{Kristoufek2013} to rank the commodities according to their efficiency.

\section{Methodology}

Efficient market can be defined in several ways. The main distinction roots back to 1965 when \cite{Fama1965} and \cite{Samuelson1965} used different definitions -- a random walk and a martingale, respectively. We stick to the martingale definition efficiency because it is less restrictive. Based on this definition, we assume that the returns of a financial asset are serially uncorrelated and with finite variance for the efficient market situation. Such a simple definition allows to use various measures of market efficiency, which are described in this section. Eventually, we refer to the Efficient Index which takes these statistics into consideration and it helps to rank different assets according to their efficiency while using various dynamic properties of the time series under study.

\subsection{Long-term memory}

Long-term memory (long-range dependence) series are characteristic with values in (even distant, in theory infinitely distant) past influencing the present and future values. These processes are standardly described with the long-term memory parameter $H$ (Hurst exponent) which ranges between $0 \le H <1$ for stationary invertible processes. The midpoint, $H=0.5$, holds for uncorrelated (or in general short-term memory) processes, i.e. processes of the efficient market. For $H>0.5$, the processes are positively correlated with long-term memory and are usually referred to as persistent. These processes systematically follow local trends while still remaining stationary. For $H<0.5$, we have long-term memory processes with negative correlations -- anti-persistent processes. Such processes switch the direction more often than a random process does.

More formally, the long-term memory processes are defined in both time and frequency domains. In the time domain, it is connected to a power-law decaying auto-correlation function. For the auto-correlation function $\rho(k)$ with time lag $k$, the decay is characterized as $\rho(k) \propto k^{2H-2}$ for $k\rightarrow +\infty$. In the frequency domain, the spectrum $f(\lambda)$ with frequency $\lambda$ of the long-range dependent process diverges at the origin so that $f(\lambda)\propto \lambda^{1-2H}$ for $\lambda \rightarrow 0+$. These definitions further lead to non-summable auto-correlations and diverging covariance of partial sums of the process for the persistent series. These properties are used in various estimators of parameter $H$. For comparison of both time and frequency domain estimators, see \cite{Beran1994,Taqqu1995,Taqqu1996,Robinson1995a,Geweke1983,DiMatteo2003,DiMatteo2007,Barunik2010,Teverovsky1999}. Out of theses estimators, we opt for the local Whittle and GPH estimators which are suitable for short time series with a possible weak short-term memory \citep{Taqqu1995,Taqqu1996}, which can easily bias the time domain estimators \citep{Teverovsky1999,Kristoufek2012}. Moreover, these estimators have well-defined asymptotic properties -- they are consistent and asymptotically normal estimators.

\subsubsection*{Local Whittle estimator}
The local Whittle estimator \citep{Robinson1995a} is a semi-parametric maximum likelihood estimator utilizing a likelihood function of \cite{Kunsch1987} and focusing only on a part of the spectrum $f(\lambda)$ near the origin. The full parametric specification is thus not needed and one does not need to assume any specific underlying long-term memory model but only a model with divergent at origin spectrum. This way, the estimator does not take into consideration high frequencies and it is in turn resistant to the short-term memory bias. As an estimator of the spectrum of series $\{x_t\}$, we use the periodogram defined as $I(\lambda_j)=\frac{1}{T}\sum_{t=1}^{T}{\exp(-2\pi i t\lambda_j)x_t}$ with $j=1,2,\ldots,m$ where $m\le T/2$ and $\lambda_j=2\pi j/T$. The local Whittle estimator is defined as
\begin{equation}
\label{eq:LWX}
\widehat{H}=\arg \min_{0\le H <1} R(H),
\end{equation}
where
\begin{equation}
\label{eq:LWX_R}
R(H)=\log\left(\frac{1}{m}\sum_{j=1}^m{\lambda_j^{2H-1}I(\lambda_j)}\right)-\frac{2H-1}{m}\sum_{j=1}^m{\log \lambda_j}.
\end{equation}
The local Whittle estimator is consistent and asymptotically normal, specifically
\begin{equation}
\sqrt{m}(\widehat{H}-H^0) \rightarrow_d N(0,1/4).
\end{equation}

\subsubsection*{GPH estimator}

Contrarily to the local Whittle estimator, the GPH estimator, named after the authors of \cite{Geweke1983}, is based on a full functional specification of the underlying process as the fractional Gaussian noise implying a specific spectral form:
\begin{equation}
\log f(\lambda) \propto -(H-0.5)\log (4\sin^2(\lambda/2))
\end{equation}
Again, the spectrum is estimated using the periodogram and the Hurst exponent is estimated using the ordinary least squares on
\begin{equation}
\label{GPH}
\log I(\lambda_j) \propto -(H-0.5)\log (4\sin^2(\lambda_j/2)).
\end{equation}
The GPH estimator is consistent and asymptotically normal \citep{Beran1994}, specifically
\begin{equation}
\sqrt{T}(\widehat{H}-H^0) \rightarrow_d N(0,\pi^2/6).
\end{equation}

The GPH estimator is thus asymptotically infinitely more efficient than the local Whittle estimator. However, this is true only if the true underlying process is in fact the fractional Gaussian noise. In financial and economic time series, this is frequently not the case as the processes are mostly a combination of short-term (such us autoregressive moving average -- ARMA -- processes of various specifications) and long-term memory (such as the aforementioned fractional Gaussian noise of fractionally integrated ARMA) processes. In this case, the GPH estimator becomes biased. To avoid the bias, the GPH estimator is based only on a part of the spectrum (periodogram) close to the origin as for the local Whittle estimator. The regression in Eq. \ref{GPH} is then applied only for a part of the periodogram based on the same parameter $m$ as for the local Whittle estimator\footnote{In our analysis, we apply $m=T^{0.6}$ \citep{Phillips2004}.}.

\subsection{Fractal dimension}

Contrary to the long-term memory, which can be seen as a characteristic of global dependence and correlation structure, the fractal dimension $D$ can be taken as a measure of local memory of the series as it is a measure of roughness of the series \citep{Kristoufek2013}. As the series can be differently rough or smooth for its specific parts, it can be locally serially correlated even though on a global level, the correlations might vanish and are not necessarily observable or detectable.

For a univariate series, the fractal dimension ranges between $1<D\le 2$. For self-similar processes, the fractal dimension is tightly connected to the Hurst exponent (long-term memory) of the series so that $D=2-H$. In economic terms, this can be understood as a perfect transmission of a local behavior (fractal dimension) to a global behavior (long-term memory). However, the relation usually does not hold perfectly for the financial series so that both $D$ and $H$ give different insights into the dynamics of the series making it worth studying them separately.

In general, $D=1.5$ holds for an uncorrelated series with no local trending or no local anti-correlations and thus it is also a value of $D$ for the efficient market. For a low fractal dimension $D<1.5$, the roughness of the series is lower than for an uncorrelated process so that we observe local trending and the series is said to be locally persistent. Reversely, a high fractal dimension $D>1.5$ characterizes a series rougher than the uncorrelated one, which is connected to local anti-persistence, i.e. the series are negatively auto-correlated locally. For purposes of the Efficiency Index introduced later in this section, we utilize Hall-Wood and Genton estimators \citep{Gneiting2004,Gneiting2010}.

\subsubsection*{Hall-Wood estimator}

Hall-Wood estimator \citep{Hall1993} is a box-counting procedure which utilizes scaling of absolute deviations between steps. Formally, we have
\begin{equation}
\widehat{A(l/n)}=\left\lfloor\frac{l}{n}\right\rfloor\sum_{i=1}^{\lfloor n/l \rfloor}{|x_{i \lfloor l/n \rfloor}-x_{(i-1)\lfloor l/n \rfloor}|}
\end{equation}
representing the absolute deviations. Using the definition of the fractal dimension \citep{Gneiting2004,Gneiting2010}, the Hall-Wood estimator is given by
\begin{equation}
\widehat{D_{HW}}=2-\frac{\sum_{l=1}^{L}{(s_l-\bar{s})\log(\widehat{A(l/n)})}}{\sum_{l=1}^{L}{(s_l-\bar{s})^2}}
\end{equation}
where $L \ge 2$, $s_l=\log(l/n)$ and $\bar{s}=\frac{1}{L}\sum_{l=1}^{L}{s_l}$. To minimize potential bias, \cite{Hall1993} propose using $L=2$ so that we obtain the estimate of the fractal dimension $\widehat{D_{HW}}$ as
\begin{equation}
\widehat{D_{HW}}=2-\frac{\log\widehat{A(2/n)}-\log\widehat{A(1/n)}}{\log2}.
\end{equation}

\subsubsection*{Genton estimator}

\cite{Gneiting2004} and \cite{Gneiting2010} propose a method of moments estimator based on the robust variogram of \cite{Genton1998}. The variogram is defined as
\begin{equation}
\widehat{V_2(l/n)}=\frac{1}{2(n-l)}\sum_{i=l}^{n}{(x_{i/n}-x_{(i-l)l/n})^2},
\end{equation}
and the Genton estimator is obtained as
\begin{equation}
\widehat{D_{G}}=2-\frac{\sum_{l=1}^{L}{(s_l-\bar{s})\log(\widehat{V_2(l/n)})}}{2\sum_{l=1}^{L}{(s_l-\bar{s})^2}}
\end{equation}
where again $L \ge 2$, $s_l=\log(l/n)$ and $\bar{s}=\frac{1}{L}\sum_{l=1}^{L}{s_l}$. \cite{Davies1999} and \cite{Zhu2002} again suggest to use $L=2$ to reduce the potential bias so that the estimate $\widehat{D_{G}}$ reads
\begin{equation}
\widehat{D_{G}}=2-\frac{\log\widehat{V_2(2/n)}-\log\widehat{V_2(1/n)}}{2\log2}.
\end{equation}

\subsection{Approximate entropy}

Entropy can be considered as a measure of complexity of the considered system. The systems with high entropy can be characterized by no information flows and are thus random up to uncertainty and reversely, the series with low entropy can be seen as deterministic \citep{Pincus2004}. The efficient market can be then seen as the one with maximum entropy and the lower the entropy, the less efficient the market is. For purposes of the Efficiency Index, we need an entropy measure which is bounded. Therefore, we utilize the approximate entropy introduced by \cite{Pincus1991}.

Let $m$ be a positive integer and let $r$ be a positive real number. For a time series $\left\{ u_{1},u_{2},...,u_{T}\right\}$, with a time series length $T$, let us form a sequence of vectors $\mathbf{X}_{1}, \mathbf{X}_{2},...,\mathbf{X}_{T-m+1}$ in $\mathbf{R}^{m}$ where $\mathbf{X}_{i}=\left( u_{i},u_{i+1},...,u_{i+m-1}\right)$. Using the Takens metrics of distance
\begin{equation}
d\left[ \mathbf{X}_{i},\mathbf{X}_{j}\right] =\underset{k=1,...,m}{\max }\left( \left\vert u\left( i+k-1\right) -u\left(j+k-1\right) \right\vert \right),
\end{equation}
and defining a characteristic function $\chi_{i}^{m}(r)$ as a number of times $d\left[ \mathbf{X}_{i},\mathbf{X}_{j}\right] \leq r/\left( T-m+1\right)$ for each $1\leq i\leq N-m+1$, we define
\begin{gather}
\Phi ^{m}(r) = \frac{1}{T-m+1} \sum_{i=1}^{T-m+1}\log \left[ \chi _{i}^{m}(r) \right]
\end{gather}
which is further used in
\begin{gather}
ER_{m}=\underset{r\rightarrow 0}{\lim }\underset{T\rightarrow \infty }{\lim }\left[ \Phi ^{m}(r) -\Phi ^{m+1}(r) \right].
\end{gather}
The approximate entropy (ApEn) is then defined as
\begin{gather}
ApEn=\underset{m\rightarrow \infty }{\lim }ER_{m}.
\end{gather}
Since $r$ can be seen as an discriminating factor for the distance measured by the Takens metrics and $m$ is the number of elements whose closeness is measured, the approximate entropy measures whether different segments of the series follow similar patterns. For an identically uniformly independently distributed random process, the approximate entropy converges to $-\log \left( r/\sqrt{3}\right) $ for all $m$ \citep{Pincus1991}. For a completely deterministic process, the entropy goes to 0. Therefore, we can rescale the approximate entropy so that $0 \le ApEn \le 1 $, where 0 characterizes a completely deterministic process and 1 a completely uncertain process characteristic for the efficient market. In turn, it can be utilized in the Efficiency Index, definition of which follows.

\subsection{Capital market efficiency measure}

\cite{Kristoufek2013} introduce the Efficiency Index (EI) is defined as

\begin{equation}
EI=\sqrt{{\sum_{i=1}^n{\left(\frac{\widehat{M_i}-M_i^{\ast}}{R_i}\right)^2}}},
\end{equation}
where $M_i$ is the $i$th measure of efficiency, $\widehat{M_i}$ is an estimate of the $i$th measure, $M_i^{\ast}$ is an expected value of the $i$th measure for the efficient market and $R_i$ is a range of the $i$th measure. In words, $EI$ is simply a distance from the efficient market situation. Here, we base the index on three measures of market efficiency -- Hurst exponent $H$ with an expected value of 0.5 for the efficient market ($M_H^{\ast}=0.5$), fractal dimension $D$ with an expected value of 1.5 ($M_D^{\ast}=1.5$) and the approximate entropy with an expected value of 1 ($M_{AE}^{\ast}=1$). Hurst exponent is taken as an average of the GPH and the local Whittle estimates. In the same way, the fractal dimension is set as an average of the Hall-Wood and Genton estimates. For the approximate entropy, we utilize the estimate described in the corresponding section. The approximate entropy need to be rescaled as it ranges between 0 and 1 with the efficient market of $ApEn=1$. We thus have $R_{AE}=2$ and $R_D=R_H=1$.

\section{Data description and results}

We analyze daily front futures prices, i.e. futures with the earliest delivery, of 25 commodities in period between 1.1.2000 and 22.7.2013\footnote{The time series were obtained from http://www.quandl.com server on 23.7.2013.}. The dataset contains 4 energy (Brent crude oil, WTI crude oil, heating oil and natural gas), 5 metals (copper, gold, silver, palladium and platinum), 7 grains (corn, oats, rough rice, soybean meal, soybean oil, soybeans and wheat), 5 softs (cocoa, coffee, cotton, orange juice and sugar) and 4 other agricultural commodities (feeder cattle, lean hogs, live cattle and lumber) futures from Chicago Board of Trade (CBOT), Chicago Mercantile Exchange (CME), IntercontinentalExchange (ICE), New York Mercantile Exchange (NYMEX) and its division Commodity Exchange (COMEX), which are summarized in Tab. \ref{tab1}. We analyze logarithmic prices $S_{i,t}=\log P_{i,t}$, where $P_{i,t}$ is the price of futures $i$ at time $t$, for the fractal dimension and logarithmic returns $r_{i,t}=S_{i,t}-S_{i,t-1}$ for the long-term memory and approximate entropy. The returns of all the analyzed futures are stationary according to ADF \citep{Dickey1979} and KPSS \citep{Kwiatkowski1992} tests (we do not report the $p$-values here).

Estimated Hurst exponents, fractal dimensions and approximate entropies are summarized in Tab. \ref{tab2}. We observe that majority of commodities is characteristic with the fractal dimension below 1.5 which indicates local persistence. These series are thus locally trending. This is most evident for feeder cattle, lean hogs and live cattle, i.e. majorly livestock futures. On the contrary, the energy commodities -- namely both the crude oils and natural gas -- are close to fractal dimension of 1.5 and as such, they do not show any signs of local inefficiencies. For the long-term memory part, most of the futures are below 0.5 indicating anti-persistence which translates into a mean-reversion of prices, something that is not standardly observed for stocks, stock indices and exchange rates which are characteristic by a unit-root behavior. The strongest anti-persistence is seen for cocoa, oats and orange juice. Nonetheless, there is a portion of commodities which show signs of persistence. These are copper, palladium, platinum and sugar. Cotton and natural gas get the closest to the value of the efficient market. For the approximate entropy, several values are close to 1 for the efficient market\footnote{Several values even reach value above 1 due to the finite sample.} -- lumber, sugar and heating oil. The most complex, and thus the least efficient, series include feeder cattle and live cattle.

Putting these results together, we arrive at the Efficiency Indices and efficiency ranking which are graphically represented in Fig. \ref{EI}. The most efficient of the commodities turns out to be heating oil closely followed by WTI crude oil. Cotton, wheat and coffee come after these with a similar level of efficiency. The ranking is then supplemented by other commodities, efficiency of which increases quite steadily across the ranking. Feeder cattle is the least efficient commodity in this set quite closely followed by live cattle. The livestock futures thus seem to be rather inefficient compared to the others. Connected to this finding, we also show an average Efficiency Index for commodities according to their type. In Fig. \ref{Means}, we can see that the energy futures are the most efficient ones followed by softs, grains and metals. By far the least efficient group consists of the other agricultural commodities, i.e. feeder cattle, lean hogs, live cattle and lumber. This is well in hand with the observations about very inefficient livestock futures.

In Figs. \ref{Parts1} and \ref{Parts2}, we decompose the efficiency index into its parts. In Fig. \ref{Parts1}, the actual futures ranked according to the Efficiency Index are represented, and in Fig. \ref{Parts2}, these are sorted according to their type to better see possible patterns and regularities. We observe that for about half of the futures, the approximate entropy is the dominant inefficiency source. Interestingly, it is the most important part for both the most and the least efficient commodities. For the others, the long-term memory part is dominant. Fractal dimension forms usually only a smaller part of the inefficiency and only for wheat, it contributes the most. When we look at the whole groups of commodities, we observe that for energy, grains and other agricultural commodities, the approximate entropy forms an important or even a dominant part for most of them. For grains, fractal dimension creates a significant part for three of the group. For softs, the long-term memory is the most important of the inefficiency contributors. And for metals, the evidence is mixed.

Fig. \ref{DH} then illustrates a relationship between fractal dimension and Hurst exponent. For self-similar processes, it holds that $D=2-H$. In economic terms, self-similar processes are characteristic by translating the local properties into the global ones. Therefore, for a locally persistent process with $D<1.5$, this translates into the global persistence with $H>0.5$, and vice versa. However, we do not observe such relationship for the analyzed commodities. Actually, the dependence is reversed so that with the increasing Hurst exponent, the fractal dimension increases. This is in contrast with the results for stock indices \citep{Kristoufek2013}. Nonetheless, such result can be well explained by characteristics of commodities futures -- locally (or in the short term), the changes in futures prices are partially predictable, but globally, the prices return to their fundamental values.

\section{Conclusion}

We have analyzed the market efficiency of 25 commodities futures across various groups -- metals, energies, softs, grains and other agricultural commodities. To do so, we have utilized the recently proposed Efficiency Index to find that the most efficient of all the analyzed commodities is heating oil, closely followed by WTI crude oil, cotton, wheat and coffee. On the other end of the ranking, we have detected live cattle and feeder cattle. The efficiency also seems to be characteristic for specific groups of commodities -- energy commodities have been found the most efficient, followed by softs, grains and metals whereas the other agricultural commodities (formed mainly of livestock) form the least efficient group. Apart from that, we have also discussed the contributions of the long-term memory, fractal dimension and approximate entropy to the total inefficiency. We have uncovered that the contribution is type-dependent as well even though the regularities are not strongly pronounced. Last but not least, we have come across the nonstandard relationship between the fractal dimension and Hurst exponent. For the analyzed dataset, the relationship between these two is positive meaning that local persistence (trending) is connected to global anti-persistence. We attribute this to specifics of commodity futures which might be predictable in a short term and locally but in a long term, they return to their fundamental price, which differs from the results found for stock indices \citep{Kristoufek2013}.

\section*{Acknowledgements}

The support from the Czech Science Foundation under Grants 402/09/0965 and P402/11/0948, and project SVV 267 504 are gratefully acknowledged.

\section*{References}


\section{Introduction\label{s:intro}}
\input introduction.tex

\section{Formulation\label{s:formulation}}
\input formulation.tex

\section{Method\label{s:method}}
\input method.tex

\section{Computing Local Averages of Pressure and Stress\label{s:aver}}
\input aver.tex

\section{Results\label{s:results}}
\input results.tex

\section{Conclusions\label{s:conclusions}}
\input conclusions.tex

\begin{appendices}
\section{Error estimates for near-singular integration \label{A:AppendixA}}
\input appen1.tex

\section{Jumps in pressure and stress \label{A:AppendixB}}
\input appen2.tex

\section{Variable curvature formulation \label{A:AppendixC}}
\input appen3.tex
\end{appendices}

\bibliographystyle{plainnat}

\subsection{Integral equation formulation}
We give the formulation for the general case in which we have viscosity
contrast between the interior and exterior fluid ($\mu_0$ is the
viscosity of the exterior fluid, $\mu_p$ is the viscosity of the
interior fluid for vesicle $p$, and $\nu_p = \mu_p/\mu_0$), and solid
boundaries with a prescribed velocity.  First we introduce $\SS_{pq}$
and ${\mathcal{D}}_{pq}$, the single- and double-layer potentials for Stokes flow,
where the constant factors are chosen so that our formulation is in
agreement with Pozrikidis~\cite[equation 2.2]{pozrikidis2001b}. The
subscripts denote the potentials induced by a hydrodynamic density on
the membrane of vesicle $q$ and evaluated on the membrane of vesicle
$p$:
\begin{align*}
  \SS_{pq}[{\mathbf{f}}]({\mathbf{x}}) &:= \frac{1}{4\pi\mu_{0}}\int_{\gamma_{q}}\left(
    -\boldsymbol{I} \log \rho
  + \frac{{\mathbf{r}} \otimes {\mathbf{r}}}{\rho^{2}} \right){\mathbf{f}}\, ds_{{\mathbf{y}}},
  && {\mathbf{x}} \in \gamma_{p}, \\
  {\mathcal{D}}_{pq}[{\mathbf{u}}]({\mathbf{x}}) &:= \frac{1-\nu_{q}}{\pi}\int_{\gamma_{q}}
    \frac{{\mathbf{r}} \cdot {\mathbf{n}}}{\rho^{2}}\frac{{\mathbf{r}} \otimes {\mathbf{r}}}{\rho^{2}}{\mathbf{u}}\,
  ds_{{\mathbf{y}}}, && {\mathbf{x}} \in \gamma_{p},
\end{align*}
where ${\mathbf{r}} = {\mathbf{x}} - {\mathbf{y}},$ and $\rho = \|{\mathbf{r}}\|_{2}$.
Also, we define
\[\SS_{p} := \SS_{pp}\quad\mbox{and}\quad {\mathcal{D}}_{p} := {\mathcal{D}}_{pp},\]
to indicate vesicle self-interactions.  Next, we define
\begin{align*}
  {\mathcal{E}}_{pq}[{\mathbf{f}},{\mathbf{u}}]({\mathbf{x}}) &= \SS_{pq}[{\mathbf{f}}]({\mathbf{x}}) + {\mathcal{D}}_{pq}[{\mathbf{u}}]({\mathbf{x}}), && {\mathbf{x}} \in \gamma_{p}, \\
  {\mathcal{E}}_{p}[{\mathbf{f}},{\mathbf{u}}]({\mathbf{x}}) &= \sum_{q=1}^{M} {\mathcal{E}}_{pq}[{\mathbf{f}},{\mathbf{u}}]({\mathbf{x}}), && {\mathbf{x}} \in \gamma_{p}.
\end{align*}
${\mathcal{B}}$ is the completed double-layer operator for confined
Stokes flow with density ${\boldsymbol\eta}$
\begin{align*}
  {\mathcal{B}}[{\boldsymbol\eta}]({\mathbf{x}}) = {\mathcal{D}}_{\Gamma}[{\boldsymbol\eta}]({\mathbf{x}}) +
    \sum_{q=1}^{M}R[\xi_{q}({\boldsymbol\eta}),{\mathbf{c}}_{q}]({\mathbf{x}}) +
    \sum_{q=1}^{M}S[{\boldsymbol\lambda}_{q}({\boldsymbol\eta}),{\mathbf{c}}_{q}]({\mathbf{x}}),
    \quad {\mathbf{x}} \in \gamma \cup \Gamma.
\end{align*}
If ${\mathbf{x}} \in \Gamma_{0}$, the rank one modification
$\mathcal{N}_{0}[{\boldsymbol\eta}]({\mathbf{x}}) = \int_{\Gamma_{0}}({\mathbf{n}}({\mathbf{x}}) \otimes
{\mathbf{n}}({\mathbf{y}})) {\boldsymbol\eta}({\mathbf{y}})ds_{{\mathbf{y}}}$ is added to ${\mathcal{B}}$ which removes the
one-dimensional null space of the corresponding integral
equation~\cite{pozrikidis1992}.  The Stokeslets and rotlets are defined
as
\begin{align*}
  S[{\boldsymbol\lambda}_{q}({\boldsymbol\eta}),{\mathbf{c}}_{q}]({\mathbf{x}}) = \frac{1}{4\pi\mu_{0}}
    \left(\log \rho + \frac{{\mathbf{r}} \otimes {\mathbf{r}}}{\rho^{2}}\right)
    {\boldsymbol\lambda}_{q}({\boldsymbol\eta})
  \qquad \text{and} \qquad
  R[\xi_{q}({\boldsymbol\eta}),{\mathbf{c}}_{q}]({\mathbf{x}}) = \frac{\xi_{q}({\boldsymbol\eta})}{\mu_{0}}
    \frac{{\mathbf{r}}^{\perp}}{\rho^{2}},
\end{align*}
where ${\mathbf{c}}_{q}$ is a point inside $\omega_{q}$, $\omega_{q}$ is the
interior of vesicle $q$, ${\mathbf{r}}={\mathbf{x}} - {\mathbf{c}}_{q}$, and ${\mathbf{r}}^{\perp} =
(r_{2},-r_{1})$.  The operator ${\mathcal{B}}$ satisfies the jump condition
\begin{align*}
  \lim_{\substack{{\mathbf{x}} \rightarrow {\mathbf{x}}_{0} \\ {\mathbf{x}} \in \Omega}}
    {\mathcal{B}}[{\boldsymbol\eta}]({\mathbf{x}}) = -\frac{1}{2}{\boldsymbol\eta}({\mathbf{x}}_{0}) +
    {\mathcal{B}}[{\boldsymbol\eta}]({\mathbf{x}}_{0}), \qquad {\mathbf{x}}_{0} \in \Gamma,
\end{align*}
and the size of the Stokeslets and rotlets are
\begin{align*}
  {\boldsymbol\lambda}_{q,i} = \frac{1}{2\pi} \int_{\gamma_{q}}
    {\boldsymbol\eta}_{i}({\mathbf{y}}) ds_{{\mathbf{y}}}, \quad i=1,2
  \qquad \text{and} \qquad
  \xi_{q} = \frac{1}{2\pi}\int_{\gamma_{q}} {\mathbf{y}}^{\perp}
    {\boldsymbol\eta}({\mathbf{y}})ds_{{\mathbf{y}}}.
\end{align*}
The inextensibility constraint is written in operator form as
\begin{align*}
  \mathcal{P}[{\mathbf{u}}]({\mathbf{x}}) = {\mathbf{x}}_{s} \cdot {\mathbf{u}}_{s}.
\end{align*}

Putting everything together, the integral formulation equation of~\eqref{e:vesicles-pde} is
given by~\cite{rah:vee:bir}
\begin{subequations}
\label{e:ves:dyn}
\begin{align}
  &(1+\nu_{p}) {\mathbf{u}}({\mathbf{x}}) = {\mathcal{E}}_{p}[{\mathbf{f}},{\mathbf{u}}]({\mathbf{x}}) +
  {\mathcal{B}}_{p}[{\boldsymbol\eta}]({\mathbf{x}}), \quad &&{\mathbf{x}} \in \gamma_{p},
  \label{e:vesicle:dynamic} \\
  &\mathcal{P}[{\mathbf{u}}]({\mathbf{x}}) = 0, &&{\mathbf{x}} \in \gamma_{p},
  \label{e:inextensibility} \\
  &\mathbf{U}({\mathbf{x}}) = -\frac{1}{2}{\boldsymbol\eta}({\mathbf{x}}) + {\mathcal{E}}_{\Gamma}[{\mathbf{f}},{\mathbf{u}}]({\mathbf{x}}) +
    {\mathcal{B}}[{\boldsymbol\eta}]({\mathbf{x}}), &&{\mathbf{x}} \in \Gamma,
  \label{e:BIE}
\end{align}
\end{subequations}
where ${\mathbf{f}}$ is the hydrodynamic density, which in our case is the stress
jump on the vesicle interface, which also is referred to as the traction
jump, and is given by
\begin{align*}
  {\mathbf{f}} = -\kappa_{b}{\mathbf{x}}_{ssss} + (\sigma {\mathbf{x}}_{s})_{s}.
\end{align*}
Since ${\mathbf{u}} = d{\mathbf{x}}/dt$ and ${\mathbf{f}}$ depends on $\sigma$ and
${\mathbf{x}}$, \eqref{e:ves:dyn} is a system of integro-differential-algebraic equations
for ${\mathbf{x}},\sigma$, and ${\boldsymbol\eta}$.

\subsection{Spatial Discretization\label{s:spaceDisc}}
Let ${\mathbf{x}}(\alpha),$ with $\alpha \in (0,2\pi]$, be a parametrization of
the interface $\gamma_p$, $\{\alpha_k = 2k\pi/n\}_{k=1}^n$ be $n$
material points, and
\begin{align*}
  {\mathbf{x}}(\alpha) = \sum_{k = -n/2+1}^{n/2} \hat{{\mathbf{x}}}(k) e^{ik\alpha}.
\end{align*}
Using the FFT to compute $\hat{{\mathbf{x}}}$, we compute derivatives of ${\mathbf{x}}$
with spectral accuracy, assuming that $\gamma$ belongs to $C^{\infty}$.
In particular, we can compute the arclength derivative with spectral
accuracy since
\begin{align*}
  \pderiv{}{s} = \pderiv{\alpha}{s}\pderiv{}{\alpha} =
    \frac{1}{\|\partial{\mathbf{x}}/\partial\alpha\|}\pderiv{}{\alpha},
\end{align*}
where $s$ is the arclength parametrization.

Since the single-layer potential $\SS$ has a logarithmic singularity, we
use the hybrid Gauss-trapezoid quadrature rule given in Table 8 of
\cite{alpert1999}.  The error of this quadrature rule is $\mathcal{O}(h^{8}
\log h)$ for integrands with logarithmic singularities. Letting ${\mathbf{x}}_k =
{\mathbf{x}}(\alpha_k)$, the quadrature rule is
\begin{align*}
  \mathcal S[{\mathbf{f}}]({\mathbf{x}}) \approx \sum_{k=1}^{n+m} w_{k}
      S({\mathbf{x}},{\mathbf{x}}_{k}) {\mathbf{f}}({\mathbf{x}}_{k})\| {\mathbf{x}}_{\alpha,k}\|,
\end{align*}
where $n$ is the number of nodes, $m$ is a number of additional
quadrature nodes, $w_{k}$ are the quadrature weights, ${\mathbf{x}}_{k}$ are the
quadrature abscissae, and ${\mathbf{x}}_{\alpha}$ is $\partial {\mathbf{x}}/\partial
\alpha$.  The first $n$ abscissae are the usual equispaced quadrature
points.  Then, the other $m$, which is determined by the desired order
of convergence for the integral, are additional quadrature abscissas
clustered around the singularity.

The double-layer potential has no singularity in two dimensions since
\begin{align*}
  \lim_{\substack{{\mathbf{x}}' \to {\mathbf{x}} \\ {\mathbf{x}}' \in \gamma}} D({\mathbf{x}}',{\mathbf{x}}) =
    \frac{\kappa({\mathbf{x}})}{2\pi}({\bf t}({\mathbf{x}}) \otimes {\bf t}({\mathbf{x}})),
  \quad {\mathbf{x}} \in \gamma,
\end{align*}
where $\kappa$ is the curvature and ${\bf t}$ is the unit tangent.
Thus, a composite trapezoid rule will give spectral accuracy since the
integrand is periodic and smooth.

\subsection{Time discretization}
Time derivatives are discretized as
\begin{align*}
  \frac{d{\mathbf{x}}}{dt} \approx \frac{\beta{\mathbf{x}}^{N+1} - {\mathbf{x}}^{0}}{\Delta t},
\end{align*}
where ${\mathbf{x}}^{0}$ is a linear combination of previous time steps.
Operators and terms that are treated explicitly are discretized at
${\mathbf{x}}^{e}$ which is also a linear combination of previous time steps.
The simplest IMEX method is IMEX Euler which is given by $\beta=1$,
${\mathbf{x}}^{0} = {\mathbf{x}}^{N}$, and ${\mathbf{x}}^{e} = {\mathbf{x}}^{N}.$ The second-order time
integrator we use is given by $\beta=3/2,$
\begin{align*}
{\mathbf{x}}^{0} = 2{\mathbf{x}}^{N} - \frac{1}{2}{\mathbf{x}}^{N-1}, \hspace{20pt}
{\mathbf{x}}^{e} = 2{\mathbf{x}}^{N} - {\mathbf{x}}^{N-1}.
\end{align*}
To avoid significant time stepping constraints, we treat the bending
term ${\mathbf{f}}_{b}$ and tension term ${\mathbf{f}}_{\sigma}$ semi-implicitly.  An
approximation for the position and tension of vesicle $p$ at time $N+1$
satisfies
\begin{subequations}
  \begin{align}
    \frac{\alpha_{p}}{\Delta t}(\beta {\mathbf{x}}_{p}^{N+1} - {\mathbf{x}}_{p}^{0}) &=
      \SS_{p}^{e} {\mathbf{f}}_{p}^{N+1} +
      {\mathcal{D}}_{p}^{e}{\mathbf{u}}_{p}^{N+1} + {\mathcal{B}}_{p}[{\boldsymbol\eta}^{N}]
      + \sum_{\substack{q=1 \\ q \neq p}}^{M} {\mathcal{E}}_{pq}^{e}[{\mathbf{f}}_{q}^{N},{\mathbf{u}}_{q}^{N}],
      &&{\mathbf{x}} \in \gamma_{p},
      \label{e:vesicle:dynamic:exp} \\
      \beta \mathcal{P}^{e}{\mathbf{x}}_{p}^{N+1} &= \mathcal{P}^{e}{\mathbf{x}}_{p}^{0},
      &&{\mathbf{x}} \in \gamma_{p},
      \label{e:inextensibility:exp} \\
      U({\mathbf{x}}) &= -\frac{1}{2}{\boldsymbol\eta}^{N}({\mathbf{x}}) + {\mathcal{E}}_{\Gamma}^{e}[{\mathbf{f}}^{N},{\mathbf{u}}^{N}]({\mathbf{x}}) +
      {\mathcal{B}}[{\boldsymbol\eta}^{N}]({\mathbf{x}}) + \mathcal{N}_{0}[{\boldsymbol\eta}^{N}]({\mathbf{x}}), &&{\mathbf{x}} \in \Gamma,
      \label{e:BIE:exp} \\
      {\mathbf{u}}_{p}^{N+1} &= \frac{\beta {\mathbf{x}}_{p}^{N+1}-{\mathbf{x}}_{p}^{0}}{\Delta t},
      &&{\mathbf{x}} \in \gamma_{p},
      \label{e:velocity:exp}
  \end{align}
  \label{e:explicit:method}
\end{subequations}
where $\alpha_{p} = (1+\nu_{p})/2$, and operators with a superscript $e$
are discretized at ${\mathbf{x}}^{e}$.  Note that~\eqref{e:explicit:method} has
explicit vesicle-vesicle and vesicle-boundary interactions.  That is, we
are approximating the solution of~\eqref{e:ves:dyn} with the solution of
a problem where the vesicles are fully decoupled from one another and
from the boundary.  To solve~\eqref{e:explicit:method}, we first find
${\boldsymbol\eta}^{N}$ by solving~\eqref{e:BIE:exp} and then
solve~\eqref{e:vesicle:dynamic:exp} and~\eqref{e:inextensibility:exp}
for the new position and tension of each vesicle independently of all
the others.

In high-concentration flows, a new source of stiffness arises.  As two
vesicles approach one another, ${\mathcal{E}}_{pq}$ becomes increasingly large.
A similar result holds for vesicle-boundary interactions where
${\mathcal{B}}_{p}$ introduces stiffness.  We propose the new method with
semi-implicit vesicle-vesicle and vesicle-boundary interactions (the
differences with \eqref{e:explicit:method} being highlighted in red)
\begin{subequations}
  \begin{align}
    \frac{\alpha_{p}}{\Delta t}
    (\beta {\mathbf{x}}_{p}^{N+1} - {\mathbf{x}}_{p}^{0}) &=
      \SS_{p}^{e} {\mathbf{f}}_{p}^{N+1} + {\mathcal{D}}_{p}^{e}{\mathbf{u}}_{p}^{N+1} +
      {\mathcal{B}}_{p}[{\boldsymbol\eta}^{N\textcolor{red}{+1}}] + \sum_{\substack{q=1 \\ q \neq p}}^{M}
      {\mathcal{E}}_{pq}^{e}[{\mathbf{f}}_{q}^{N\textcolor{red}{+1}},
      {\mathbf{u}}_{q}^{N\textcolor{red}{+1}}],
      &&{\mathbf{x}} \in \gamma_{p},
      \label{e:vesicle:dynamic:imp} \\
      \beta \mathcal{P}^{e}{\mathbf{x}}_{p}^{N+1} &= \mathcal{P}^{e}{\mathbf{x}}_{p}^{0},
      &&{\mathbf{x}} \in \gamma_{p},
      \label{e:inextensibility:imp} \\
      U({\mathbf{x}}) &= -\frac{1}{2}{\boldsymbol\eta}^{N\textcolor{red}{+1}}({\mathbf{x}}) +
      {\mathcal{E}}_{\Gamma}^{e}[{\mathbf{f}}^{N\textcolor{red}{+1}},
      {\mathbf{u}}^{N\textcolor{red}{+1}}]({\mathbf{x}}) +
      {\mathcal{B}}[{\boldsymbol\eta}^{N\textcolor{red}{+1}}]({\mathbf{x}})+\mathcal{N}_{0}[{\boldsymbol\eta}^{N\textcolor{red}{+1}}]({\mathbf{x}}), &&{\mathbf{x}} \in \Gamma,
      \label{e:BIE:imp} \\
      {\mathbf{u}}_{p}^{N+1} &= \frac{\beta {\mathbf{x}}_{p}^{N+1}-{\mathbf{x}}_{p}^{0}}{\Delta t},
      &&{\mathbf{x}} \in \gamma_{p}.
      \label{e:velocity:imp}
  \end{align}
  \label{e:implicit:method}
\end{subequations}
Unlike in~\eqref{e:explicit:method}, vesicles
in~\eqref{e:implicit:method} are implicitly coupled since the operator
${\mathcal{E}}_{pq}^{e}$ is applied to ${\mathbf{f}}_{q}^{N+1}$ and ${\mathbf{u}}_{q}^{N+1}$.
Similarly, the operator ${\mathcal{B}}_{p}$ implicitly couples the vesicles to the
solid walls.  Equation~\eqref{e:implicit:method} is still an
approximation of~\eqref{e:ves:dyn}, but its benefit is that it is
experimentally a more stable method (see Sections~\ref{s:taylorGreen}
and~\ref{s:couette}).

Since~\eqref{e:implicit:method} is fully coupled, it is too expensive to
solve without preconditioning.  Roughly speaking, the number of GMRES
iterations depends on the sizes of ${\mathcal{E}}_{pp}$ and ${\mathcal{E}}_{pq}$ and the
latter depends on the inter-vesicle proximity.  We apply preconditioned
GMRES~\cite{saad} to~\eqref{e:implicit:method} where we use the
block-diagonal preconditioner.  The number of GMRES iterations depends
on the minimum inter-vesicle distance, but this preconditioned linear
system experimentally results in mesh-independence\footnote{We
experimented with multigrid preconditioners~\cite{quaife:biros2013} but
mesh-independence was only achieved for computationally expensive
preconditioners.}.

To test the preconditioner, one step of a Taylor-Green flow is done
with the number of points per vesicle ranging from $16$ to $128$
(Table~\ref{t:GMRES:preco}).  We solve~\eqref{e:implicit:method} to a
tolerance of $10^{-8}$ with and without preconditioning.  The required
GMRES iterations for the block-diagonal solver grows almost linearly
and the results are consistent with those reported
in~\cite{rah:vee:bir}.  The implicit solver without preconditioning
requires an unacceptable number of GMRES iterations, but with
preconditioning, we see that the required number of GMRES steps is
independent of $N$ \footnote{Notice, however, that the number of GMRES
iterations will depend on how close the vesicles are to one another as
this determines the overall spectrum structure of the operator. But for
a given flow, the method is independent of $N$.}.  Using the
block-diagonal preconditioner will not result in a speedup unless we
precompute and store the block diagonal preconditioner.  We perform
this precomputation for all the experiments in Section~\ref{s:results}.

\begin{table}[htps]
\centering
\begin{tabular}{>{\centering}m{0.5cm} >{\centering}m{3cm} >{\centering}m{3cm} >{\centering}m{3cm}}
$N$ & Block- \\ Diagonal & Unpreconditioned \\ Implicit & Preconditioned
\\ Implicit \tabularnewline
\hline
16 & 46 & 315 & 18 \tabularnewline
32 & 79 & 341 & 15 \tabularnewline
64 & 164 & 882 & 15 \tabularnewline
128 & 294 & 2117 & 15
\end{tabular}
\mcaption{The number of GMRES iterations required to take one time step
of a Taylor-Green flow using semi-implicit inter-vesicle
interactions~\eqref{e:implicit:method}.  $N$ is the number of
discretization points per vesicle, the second column is the number of
GMRES iterations to invert the block-diagonal preconditioner, the third
column is the number of unpreconditioned GMRES iterations, and the
final column is the number of preconditioned GMRES iterations.  Notice
that without preconditioning, the number of iterations quickly grows
while the preconditioned linear system is
mesh-independent.}{t:GMRES:preco}
\end{table}
We have also observed experimentally that when vesicles are close,
treating inter-vesicle interactions explicitly with second-order time
stepping is unconditionally unstable.  We have experimented with time
step sizes of varying orders of magnitude and each experiment
inevitably becomes unstable.  However, scheme~\eqref{e:implicit:method}
with second-order time stepping is stable once $\Delta t$ is small
enough for the simulation under consideration.

\subsection{Fast summation}
For high-concentration flows, the most expensive part of the simulation
is evaluating the integral operators $\SS$ and ${\mathcal{D}}$.  The single-layer
potential $\SS$ can be written as
\begin{align*}
  \frac{1}{4\pi\mu_{0}} \int_{\gamma}  \left(
    -\boldsymbol{I}\log \rho  +
    \frac{{\mathbf{r}} \otimes {\mathbf{r}}}{\rho^{2}} \right){\mathbf{f}} ds =
    -\frac{1}{4\pi\mu_{0}} \int_{\gamma} (\log \rho) {\mathbf{f}} ds +
    \frac{1}{4\pi\mu_{0}} \int_{\gamma}
    \left(\frac{{\mathbf{r}} \cdot {\mathbf{f}}}{\rho^{2}}\right){\mathbf{r}} ds \\
  =-\frac{1}{4\pi\mu_{0}} \int_{\gamma}(\log \rho) {\mathbf{f}} ds +
    \frac{1}{4\pi\mu_{0}} \int_{\gamma}
    \frac{{\mathbf{r}}}{\rho^{2}} ({\mathbf{x}} \cdot {\mathbf{f}}) ds_{{\mathbf{y}}} -
    \frac{1}{4\pi\mu_{0}} \int_{\gamma}
    \frac{{\mathbf{r}}}{\rho^{2}} ({\mathbf{y}} \cdot {\mathbf{f}}) ds_{{\mathbf{y}}},
\end{align*}
which can be computed with three Laplace FMM applications.  We also use
the FMM to evaluate the single-layer potential at arbitrary target
locations by assigning a charge of zero to these locations.
Unfortunately, the double-layer potential can not be decomposed into a
sequence of Laplace single- and double-layer potentials.  Future work
involves using the kernel-independent FMM~\cite{ying-biros-zorin-03} to
accelerate evaluating the double-layer potential.

\subsection{Near-singular integration}\label{s:near-singular}
Another difficulty with high-concentration simulations is the
evaluation of layer potentials at points close to a boundary.  Integral
operators ${\mathcal{D}}$ and $\SS$ can be approximated with high accuracy using
the trapezoid rule if the target point is sufficiently far from the
source points.  As the target point approaches the vesicle,
more points are required to achieve the same accuracy, and, the
required number of points for a fixed error grows without
bound (Appendix~\ref{A:AppendixA}).  A near-singular integration
strategy that guarantees a uniform accuracy for all target points is
required.

We adopt the strategy from~\cite{ying-biros-zorin06} since it is not
depend on the nature of the integrand and it is simple to implement.
Suppose that a vesicle $\gamma$ is discretized with $N$ points and that
the resulting arclength term is of size $h$.  Let $d({\mathbf{x}},\gamma) =
\inf_{{\mathbf{y}} \in \gamma} \|{\mathbf{x}} - {\mathbf{y}}\|$ be the distance from ${\mathbf{x}}$ to
$\gamma$.  We define the far zone of $\gamma$ as $\Omega_{1}=\{{\mathbf{x}}
\:|\: d({\mathbf{x}},\gamma) \geq h\}$, and the near zone of $\gamma$ as
$\Omega_{0}=\{{\mathbf{x}} \:|\: d({\mathbf{x}},\gamma) < h\}.$ The near and far zones
can be constructed efficiently using a quadtree-structure, similar to
the near-field and far-field interactions in fast multipole codes.

For ${\mathbf{x}} \in \Omega_{1}$, the trapezoid rule with $N^{3/2}$ points
guarantees that the error for $\SS$ is $\mathcal{O}(h^{M/2-2})$ and for
${\mathcal{D}}$ is $\mathcal{O}(h^{M/2-4})$, where the density function (in our
case, either ${\mathbf{f}}$ or ${\boldsymbol\eta}$) belongs to $C^{M}$.  These error bounds
are proved in Appendix~\ref{A:AppendixA}.  For a target point ${\mathbf{x}}$ in
the near zone, we first find the closest boundary point, ${\mathbf{x}}_{0}$,
using two applications of Newton's method with the initial guess being
the closest discretization point to ${\mathbf{x}}$.  We then place $m$
interpolation points at the points ${\mathbf{x}}_{j} = {\mathbf{x}}_{0} + j\beta h
({\mathbf{x}}-{\mathbf{x}}_{0})/\|{\mathbf{x}}-{\mathbf{x}}_{0}\|$ for $j=0,\ldots,m$, where $\beta$ is
slightly larger than $1$ to guarantee that all the interpolation points
are in the far zone (left plot of Figure~\ref{f:nearsing:conv}).  The
layer potential is evaluated at ${\mathbf{x}}_{0}$ using a local interpolant with
$N_{\mathrm{int}}$ points on $\gamma$ and evaluated at ${\mathbf{x}}_{j},$
$j=1,\ldots,m$ using the $N^{3/2}$-point trapezoid rule.  Then, a
one-dimensional Lagrange interpolation assigns a value to the layer
potential at ${\mathbf{x}}$.  Note that the proposed scheme can become more
accurate by using interpolation schemes that compute derivatives of the
solution and building more accurate polynomial interpolations.  This,
however, would significantly complicate the implementation.  Pseudocode
for this algorithm is given in Algorithm~\ref{a:algorithm2}.  In
Appendix~\ref{A:AppendixA}, we show that the error of our algorithm
applied to the single-layer potential $\SS$ is
\begin{align*}
  \mathcal{O}\left(h^{\min\left(N_{\mathrm{int}}-1,
      m,\frac{M}{2}-2\right)}\right),
\end{align*}
and applied to the double-layer potential ${\mathcal{D}}$ is
\begin{align*}
  \mathcal{O}\left(h^{\min\left(N_{\mathrm{int}}-1,
      m,\frac{M}{2}-4\right)}\right).
\end{align*}

\begin{algorithm}
\begin{algorithmic}
\FORALL[Loop over far zone]{${\mathbf{x}} \in \Omega_{1}$}
\STATE ${\mathbf{u}} = \mathtt{trapezoid}(\gamma,{\mathbf{x}},{\mathbf{f}})$
\COMMENT{Apply the $N^{3/2}$-point trapezoid rule}
\ENDFOR
\FORALL[Loop over near zone]{${\mathbf{x}} \in \Omega_{0}$}
\STATE ${\mathbf{x}}_{0} = \argmin d({\mathbf{x}},\gamma)$
\COMMENT{Find the point on $\gamma$ closest to ${\mathbf{x}}$}
\STATE ${\mathbf{u}}_{0} = \mathtt{interpolate}
    ({\mathbf{x}}_{0},{\mathbf{u}}_{1},\ldots,{\mathbf{u}}_{N_{\mathrm{int}}})$
\COMMENT{Interpolate ${\mathbf{u}}$ at ${\mathbf{x}}_{0}$}
\FOR{$j=1,\ldots,m$}
\STATE ${\mathbf{x}}_{j} = {\mathbf{x}}_{0} + jh\beta({\mathbf{x}}-{\mathbf{x}}_{0})/\|{\mathbf{x}}-{\mathbf{x}}_{0}\|$
\COMMENT{Define Lagrange interpolation points}
\STATE ${\mathbf{u}}_{j} = \mathtt{trapezoid}(\gamma,{\mathbf{x}}_{j},{\mathbf{f}})$
\COMMENT{Evaluate layer potential at ${\mathbf{x}}_{j}$}
\ENDFOR
\STATE ${\mathbf{u}} = \mathtt{interpolate}({\mathbf{x}},{\mathbf{u}}_{0},\ldots,{\mathbf{u}}_{m})$
\COMMENT{Compute ${\mathbf{u}}$ using Lagrange interpolation}
\ENDFOR
\RETURN ${\mathbf{u}}$ \COMMENT{Return the layer potential at all target points}
\end{algorithmic}
\mcaption{LayerPotential($\gamma$,  ${\mathbf{f}}$, ${\mathbf{u}}$, $\Omega_{0}$,
$\Omega_{1}$)} {a:algorithm2}
\end{algorithm}

\begin{figure}[htps]
  \begin{center}
    \begin{minipage}{0.45\textwidth}
    \input{near-sing.tikz}
    \end{minipage}
    \begin{minipage}{0.45\textwidth}
    \input{nearSingErrors.tikz}
    \end{minipage}
  \end{center}
  \mcaption{\emph{Left:} Lagrange interpolation points in $\Omega_{1}$
  and $\gamma$ are used to approximate layer potentials at ${\mathbf{x}} \in
  \Omega_{0}$.  The value at ${\mathbf{x}}_{0}$ uses an interpolant from
  neighboring points on the boundary.  All other interpolation points
  are in $\Omega_{1}$ so that the $N^{3/2}$-point trapezoid rule is
  sufficiently accurate.  \emph{Right:} Here we report convergence
  rates. For the blue and the black curves, the expected order of
  convergence is $5$. We obtain $4.7$ and $5.0$, respectively. As the
  green curve corresponds to a linear flow, Lagrange interpolation gives
  exact values.  Hence, machine accuracy is reached once $N$ is
  sufficiently large.}{f:nearsing:conv}
\end{figure}

We test our near-singular integration with three examples.  In all the
examples, 32 target points are placed in $\Omega_{0}$ for
$N=16,\ldots,512$.  The boundary is a $3:2$ ellipse, the boundary values
are interpolated with $N_{\mathrm{int}}=6$ points, and the Lagrange
interpolation uses $m=6$ points.  Since $m = N_{\mathrm{int}}$, we are
using one more point to interpolate in the direction coming off of
$\gamma$.  We do this because the interpolant that is coming off of the
vesicle is always evaluated near the first interpolation point.  Thus,
it may suffer from the Runge phenomenon (we have observed this
experimentally).  This is not the case for the local interpolation used
to find ${\mathbf{u}}({\mathbf{x}}_{0})$ since ${\mathbf{x}}_{0}$ is located near the middle of the
local interpolation points.  For the first two examples, we compute the
density function of $\SS$ for the two solutions of the Stokes equations:
\begin{align*}
  &{\mathbf{u}}({\mathbf{x}})=(x,-y), \\
  &{\mathbf{u}}({\mathbf{x}})=\left(-\log|{\mathbf{x}}| + \frac{{\mathbf{x}} \otimes {\mathbf{x}}}{|{\mathbf{x}}|^{2}}\right)
  \left[ \begin{array}{c}
    1 \\ 1
  \end{array}\right].
\end{align*}
For the third example, we pick an arbitrary density function and compare
the results from near-singular integration with the results from an
over-refined trapezoid rule.  The error for all three examples is
reported in~Figure~\ref{f:nearsing:conv}.  For the extensional flow ${\mathbf{u}}
= (x,-y)$, the Lagrange interpolation is exact since the flow is linear
and the error decays faster than the other examples.  For the other
examples, the error is close to fifth-order which agrees with the order
of accuracy of the near-singular integration

If we have $M$ vesicles and $N$ points per vesicle, the major costs of
the near-singular integration algorithm are:
\begin{itemize}
  \item \emph{Upsampling each vesicle from $N$ to $N^{3/2}$ points:}
  Using the FFT, this requires $\mathcal{O}(MN^{3/2}\log N)$ operations.
  \item \emph{Applying the $N^{3/2}$-point quadrature rule to
  $\Omega_{1}$:} With the FMM, this requires $\mathcal{O}(MN^{3/2} + MN)$
  operations, and without the FMM, this requires $\mathcal{O}(M^{2}N^{5/2})$
  operations.
  \item \emph{Evaluating the layer potential in $\Omega_{0}$:} Assuming
  there are $\mathcal{O}(1)$ points in the near zone, this requires
  $\mathcal{O}(MN^{3/2})$ operations.
\end{itemize}
The algorithm can be further sped up by introducing the intermediate and
far zones~\cite{ying-biros-zorin06}.  The intermediate zone is
$\Omega_{1}=\{{\mathbf{x}} \:|\: d({\mathbf{x}},\gamma) \in (\sqrt{h},h]\}$ and the far
zone is $\Omega_{2}=\{{\mathbf{x}} \:|\: d({\mathbf{x}},\gamma) \geq \sqrt{h}\}$.  Then,
the same accuracy can be achieved by using the $N^{3/2}$-point trapezoid
rule in $\Omega_{1}$ and the $N$-point trapezoid rule in $\Omega_{2}$.
While this would reduce the constant of the complexity, it does not
change the overall order and complicates the implementation.

\subsection{Collision detection}\label{s:collision}
Even with semi-implicit treatment of inter-vesicle interactions, due
to time discretization errors, the vesicles can come so close that they
collide or cross a solid wall.\footnote{The streamlines of Stokes
equations never intersect.  Therefore, it is physically impossible for
vesicles to collide with one another or with solid walls in finite
time.  Lubrication theory can be used to study the forces between
nearly touching vesicles, and therefore compute estimates of the
distance between the vesicles, but we choose not to discuss these
forces here.}  There has been a lot of work in collision
detection~\cite{jimenez-e13,redon-e04}, but they are not spectrally
accurate or require complex surface-surface intersections.  Most
methods construct piecewise approximations to the boundaries to
simplify the detection of collisions.  Since our vesicles and solid
walls are periodic, we can use some nice properties of layer potentials
to create a spectrally accurate collision detector.  We use the classic
potential theory result~\cite{kellogg}
\begin{align}
  I_{\Gamma}({\mathbf{x}}) = \frac{1}{2\pi}\int_{\Gamma}\pderiv{}{{\mathbf{n}}_{{\mathbf{y}}}}\log|{\mathbf{x}} - {\mathbf{y}}| ds_{{\mathbf{y}}} =
  \begin{cases}
    1 & {\mathbf{x}} \in \Omega, \\
    \frac{1}{2} & {\mathbf{x}} \in \Gamma, \\
    0 & {\mathbf{x}} \in \mathbf{R}^{2} \backslash \overline{\Omega},
  \end{cases}
  \label{e:layerpotentials}
\end{align}
where $\Omega \subset \mathbf{R}^{2}$ is bounded with a smooth boundary
$\Gamma$, and ${\mathbf{n}}_{{\mathbf{y}}}$ is the outward pointing normal.  We
apply~\eqref{e:layerpotentials} to the configuration of the vesicles
and the solid walls by computing $I_{X}({\mathbf{x}}) = I_{\gamma_{1}}({\mathbf{x}}) +
\cdots + I_{\gamma_{M}} + I_{\Gamma}$ for all ${\mathbf{x}} \in \gamma$.  The
maximum value of $I_{X}({\mathbf{x}})$ is $3/2$ for confined flows and $1/2$ for
unconfined flows if vesicles have not crossed and have not left the
domain.  The maximum value of $I_{X}({\mathbf{x}})$ will be at least $3/2$ for
unconfined flows, and at least $5/2$ for confined flows if vesicles
have crossed or left the domain.  We note that $I_{X}({\mathbf{x}})$ can be
evaluated accurately using our near-singular integration scheme, and in
linear time using the FMM.

We can use $I_{X}({\mathbf{x}})$ to create a simple adaptive time-stepping
method.  At each time step, we check for collisions.  If a collision has
occurred, we initialize the simulation at a previous time step with a
time step half the size.  We are currently implementing a more rigorous
adaptive time stepping method that uses ideas from spectral deferred
correction methods~\cite{dut:gre:rok,hua:jia:min,min}.  At each time
step, we estimate the local truncation error and then try to commit a
constant amount of error per time step by adjusting the time step size.
We leave the details and the implementation as future work.  Future work
also includes adaptivity in space.  If vesicles are sufficiently far
from one another, their hydrodynamic interaction can be represented with
a moderate number of points, but as they approach one another, the
hydrodynamic interaction strengthens and this requires a finer spatial
resolution.

\subsection{Shear Flow}
Here we consider two vesicles in the shear flow $(y,0)$.  We initially
place the left vesicle's center slightly above the x-axis and place the
right vesicle's center exactly on the x-axis.  Initially, we expect the
right vesicle to undergo tank-treading~\cite{kantsler,misbah2006} while
the other vesicle will travel from left to right.  Here we study the
effect of the reduced area and the viscosity contrast on the dynamics
and errors.  We expect vesicles with larger viscosity contrast and
larger reduced areas to behave more like rigid bodies.  In particular,
we expect them to come closer to one another.
Figure~\ref{f:shear:reducedAreas} indicate that a higher viscosity
contrast results in the vesicles coming closer together.  Moreover, for
larger reduced areas, the effect of the viscosity contrast is less
pronounced.  We also see that vesicles with a larger reduced area come
closer together.

\begin{figure}[htps]
\centering
\begin{tabular}{m{1cm}CCCCC}
$\nu=1$ &
\input{shearOrder1RA1viscCont1Time1.tikz} &
\input{shearOrder1RA1viscCont1Time2.tikz} &
\input{shearOrder1RA1viscCont1Time3.tikz} &
\input{shearOrder1RA1viscCont1Time4.tikz} &
\input{shearOrder1RA1viscCont1Time5.tikz} \\
$\nu=4$ &
\input{shearOrder1RA1viscCont4Time1.tikz} &
\input{shearOrder1RA1viscCont4Time2.tikz} &
\input{shearOrder1RA1viscCont4Time3.tikz} &
\input{shearOrder1RA1viscCont4Time4.tikz} &
\input{shearOrder1RA1viscCont4Time5.tikz} \\
$\nu=1$ &
\input{shearOrder1RA2viscCont1Time1.tikz} &
\input{shearOrder1RA2viscCont1Time2.tikz} &
\input{shearOrder1RA2viscCont1Time3.tikz} &
\input{shearOrder1RA2viscCont1Time4.tikz} &
\input{shearOrder1RA2viscCont1Time5.tikz} \\
$\nu=4$ &
\input{shearOrder1RA2viscCont4Time1.tikz} &
\input{shearOrder1RA2viscCont4Time2.tikz} &
\input{shearOrder1RA2viscCont4Time3.tikz} &
\input{shearOrder1RA2viscCont4Time4.tikz} &
\input{shearOrder1RA2viscCont4Time5.tikz}
\end{tabular}
\begin{tabular}{cc}
\input{shearRA1GapSizeOrder1.tikz} &
\input{shearRA2GapSizeOrder1.tikz}
\end{tabular}
\mcaption{Two vesicles are submerged in a shear flow.  In the first and
second rows, the vesicles have reduced area 0.65, and in the third and
fourth rows, the vesicles have reduced area 0.99.  Each simulation is
done with viscosity contrasts $\nu=1$ and $\nu=4$.  The bottom plots
show the distance between the vesicles for both reduced areas and both
viscosity contrasts.  Comparing semi-implicit inter-vesicle interactions
with both first- and second-order time stepping, the distance between
the vesicles never differ by more than 2.2\% for both reduced areas and
both viscosity contrasts.}{f:shear:reducedAreas}
\end{figure}

We now do a convergence study for both values of the reduced area and
viscosity contrast.  We vary the number of points per vesicle, the
number of time steps, and the time integrator.  We report the errors in
area and length at the time horizon $t=12$ in Tables~\ref{t:shear:RA65}
and~\ref{t:shear:RA99}.  First, we see that first- and second-order
convergence is achieved for all combinations of reduced area and
viscosity contrast.  Second, semi-implicit inter-vesicle interactions
do not improve the accuracy of first-order time stepping, but, we will
see later that this time integrator is more stable.

\begin{table}[htps]
\centering
\begin{tabular}{ccc|cc|cc}
Time & & & \multicolumn{2}{c|}{$\nu=1$} & \multicolumn{2}{c}{$\nu=4$} \\
Integrator   & $N$ & $\Delta t$ & $e_{A}$   & $e_{L}$   &
$e_{A}$   & $e_{L}$ \\
\hline

Explicit (1) & 32  & 0.04       & $6.73e-2$ & $5.21e-2$ &
$4.07e-2$ & $2.41e-2$ \\

Explicit (1) & 64  & 0.02       & $3.36e-2$ & $2.62e-2$ &
$2.02e-2$ & $1.20e-2$ \\

Explicit (1) & 128 & 0.01       & $1.68e-2$ & $1.31e-2$ &
$1.01e-2$ & $5.99e-2$ \\
\hline

Implicit (1) & 32  & 0.04       & $6.73e-2$ & $5.21e-2$ &
$4.06e-2$ & $2.40e-2$ \\

Implicit (1) & 64  & 0.02       & $3.36e-2$ & $2.62e-2$ &
$2.02e-2$ & $1.20e-2$ \\

Implicit (1) & 128 & 0.01       & $1.68e-2$ & $1.31e-2$ &
$1.02-2$ & $5.98e-3$ \\
\hline

Implicit (2) & 32  & 0.04       & $1.59e-4$ & $3.21e-4$ &
$2.69e-5$ & $6.33e-5$ \\

Implicit (2) & 64  & 0.02       & $3.87e-6$ & $6.52e-5$ &
$3.78e-6$ & $1.20e-5$ \\

Implicit (2) & 128 & 0.01       & $2.79e-7$ & $1.47e-5$ &
$3.85e-7$ & $2.28e-6$
\end{tabular}
\mcaption{The errors in area and length at $t=12$ for different time
integrators applied to two vesicles of reduced area $0.65$ in a shear
flow.  The parenthetic values are the time stepping order.  Results are
presented for no viscosity contrast and a viscosity contrast of
4.}{t:shear:RA65}

\begin{tabular}{ccc|cc|cc}
Time & & & \multicolumn{2}{c|}{$\nu=1$} & \multicolumn{2}{c}{$\nu=4$} \\
Integrator   & $N$ & $\Delta t$ & $e_{A}$   & $e_{L}$   &
$e_{A}$   & $e_{L}$ \\
\hline

Explicit (1) & 32  & 0.04       & $1.05e-1$ & $5.48e-2$ &
$9.67e-2$ & $4.94e-2$ \\

Explicit (1) & 64  & 0.02       & $5.15e-2$ & $2.70e-2$ &
$4.78e-2$ & $2.45e-2$ \\

Explicit (1) & 128 & 0.01       & $2.55e-2$ & $1.34e-2$ &
$2.37e-2$ & $1.22e-2$ \\
\hline

Implicit (1) & 32  & 0.04       & $1.05e-1$ & $5.48e-2$ &
$9.69e-2$ & $4.93e-2$ \\

Implicit (1) & 64  & 0.02       & $5.16e-2$ & $2.70e-2$ &
$4.77e-2$ & $2.45e-2$ \\

Implicit (1) & 128 & 0.01       & $2.55e-2$ & $1.34e-2$ &
$2.36e-2$ & $1.22e-2$ \\
\hline

Implicit (2) & 32  & 0.04       & $2.40e-5$ & $1.13e-4$ &
$5.85e-5$ & $4.68e-5$ \\

Implicit (2) & 64  & 0.02       & $1.65e-6$ & $1.50e-5$ &
$7.05e-6$ & $5.86e-6$ \\

Implicit (2) & 128 & 0.01       & $1.56e-7$ & $1.92e-6$ &
$8.62e-7$ & $7.26e-7$

\end{tabular}

\mcaption{The errors in area and length at $t=12$ for different time
integrators applied to two vesicles of reduced area $0.99$ in a shear
flow.  The parenthetic values are the time stepping order.  Results are
presented for no viscosity contrast and a viscosity contrast of
4.}{t:shear:RA99}
\end{table}

\subsection{Taylor-Green Flow}
\label{s:taylorGreen}
We consider the Taylor-Green flow ${\mathbf{u}} = (\sin x \cos y,-\cos x \sin
y)$ with 9 vesicles each with $N=64$ points of reduced area 0.65 placed
inside the periodic cell $(0,2\pi)^{2}$ (left plot of
Figure~\ref{f:taylorGreen:collision}).  We consider first-order IMEX
Euler with both semi-implicit and explicit inter-vesicle interactions.
We find the largest time step of both methods so that the vesicles do
not cross before a time horizon of $T = 1$.  With explicit
interactions, the vesicles cross if $74$ time steps are taken but not
if $75$ time steps are taken.  However, with semi-implicit
interactions, the simulation successfully completes even with one time
step.  Of course the error is too large and a time step size this large
is unreasonable.

\begin{figure}[htps]
\begin{center}
  \begin{tabular}{ccc}
    \input{taylorGreenSnapsExplicitTime1.tikz} &
    \input{taylorGreenSnapsExplicitTime2.tikz} &
    \input{taylorGreenSnapsImplicitTime2.tikz} \\
    &
    \input{taylorGreenSnapsExplicitTime2Zoomed.tikz} &
    \input{taylorGreenSnapsImplicitTime2Zoomed.tikz}
  \end{tabular}
\end{center}
\mcaption{We compare first-order time stepping with explicit and
semi-implicit inter-vesicle interactions for a time step size of 0.02
and 64 points per vesicle with the initial configuration in the
left-most plot.  With explicit interactions (middle plot), the vesicles
cross after only 13 time steps.  However, with semi-implicit
interactions (right plot), the vesicles do not cross.  The bottom plots
show magnifications of the region where the explicit solver first
crosses.}{f:taylorGreen:collision}
\end{figure}

For this example, we found that using explicit interactions with 75
time steps requires approximately the same amount of CPU time as using
semi-implicit interactions with 30 time steps.  Therefore, if one
permits the error to be approximately 2.5 times larger, using
semi-implicit inter-vesicle interactions is justified.

We do a convergence study for first- and second-order time stepping
with semi-implicit inter-vesicle interactions.  If these interactions
are treated explicitly, the coarsest simulation results in the vesicles
crossing.  We achieve the desired first-order convergence, but
second-order convergence is not achieved at these resolutions.
However, we ran the simulation with $\Delta t = 3.13e-4$ and $\Delta t
= 1.56e-4$ and this is where second-order convergence is first
observed; the error in area reduced by $4.67$ and the error in length
reduced by $3.96$.

We have also observed a qualitative difference in the
behaviour of the errors of first- and second-order time stepping.  The
errors of the first-order methods continue to grow throughout the
simulation, while for second-order time stepping, the errors
plateau.  Thus, given a fixed time step size, second-order time stepping
can simulate longer time horizons than first-order time stepping.  We
demonstrate this behaviour in the bottom plot of
Figure~\ref{f:taylorGreen:summary} where we plot the errors in area and
length for the first- and second-order time stepping.  Snapshots of the
simulation are in the top plots of Figure~\ref{f:taylorGreen:summary}.

\begin{table}[htps]
\begin{center}
\begin{tabular}{cc|cc|cc}
& & \multicolumn{2}{c|}{First-order} & \multicolumn{2}{c}{Second-order} \\
$N$ & $\Delta t$ & $e_{A}$ & $e_{L}$ & $e_{A}$ & $e_{L}$ \\
\hline
32  & $2e-2$ & $5.24e-2$ & $3.03e-2$ & $5.25e-4$ & $3.98e-3$ \\
64  & $1e-2$ & $2.77e-2$ & $1.64e-2$ & $2.08e-4$ & $1.96e-3$ \\
128 & $5e-3$ & $1.43e-2$ & $8.68e-3$ & $7.76e-5$ & $9.42e-4$ \\
256 &$2.5e-3$& $7.24e-3$ & $4.53e-3$ & $2.73e-5$ & $4.20e-4$ \\
\end{tabular}
\mcaption{The errors in area and length at $t=5$ for a Taylor-Green
flow.  All runs were done with semi-implicit inter-vesicle interactions
since explicit interactions with $N=32$ and $\Delta t = 2e-2$ results
in crossing vesicles.}{t:taylorGreen:convergence}
\end{center}
\end{table}

\begin{figure}[htps]
\begin{center}
  \begin{tabular}{cccccc}
    \input{taylorGreenSnapsTime1.tikz} &
    \input{taylorGreenSnapsTime2.tikz} &
    \input{taylorGreenSnapsTime3.tikz} &
    \input{taylorGreenSnapsTime4.tikz} &
    \input{taylorGreenSnapsTime5.tikz} &
    \input{taylorGreenSnapsTime6.tikz}
  \end{tabular}
  \begin{tabular}{cc}
    \input{taylorGreenErrorsOrder1.tikz} &
    \input{taylorGreenErrorsOrder2.tikz}
  \end{tabular}
\end{center}
\mcaption{The top plots show several snapshots of a Taylor-Green
simulation with a long time horizon using second-order time stepping.
In the bottom plots are the errors in area and length for the first- and
the second-order time integrators.  Notice that the errors for the
first-order method continue to grow and the simulation stops near $t=20$
(the simulation is stopped when a 10\% error has been incurred), whereas
the errors for the second-order method plateau and the simulation is
able to run much longer.}{f:taylorGreen:summary}
\end{figure}

\subsection{Extensional Flow}
We consider two vesicles with reduced area $0.65$ symmetrically placed
around the origin (top plots of Figure~\ref{f:extensional:summary}).
The background velocity is given by ${\mathbf{u}} = (-x,y)$.  We first run the
simulation with and without near-singular integration.  We compute the
distance between the vesicles, and the errors in area and length
(bottom plots of Figure~\ref{f:extensional:summary}).  We see that when
$t \approx 5$, the simulation without near-singular integration begins
to introduce large errors in the area of the vesicles.  This indicates
that the approximation of the single-layer potential is not adequate.
Then, around $t \approx 12$, the vesicles are too close for the
$N$-point trapezoid rule and the simulations fails.  This agrees with
our error estimate from Appendix~\ref{A:AppendixA}  since for this
simulation, an arclength term is approximately 0.24.  However, with
near-singular integration, the simulation successfully runs to
completion.

\begin{figure}[htps]
\begin{center}
  \begin{tabular}{ccccc}
  \input{extensionalSnapsTime1.tikz} &
  \input{extensionalSnapsTime2.tikz} &
  \input{extensionalSnapsTime3.tikz} &
  \input{extensionalSnapsTime4.tikz} &
  \input{extensionalSnapsTime5.tikz}
  \end{tabular}
  \begin{tabular}{cc}
  \input{extensionalGapSize.tikz} &
  \input{extensionalErrors.tikz}
  \end{tabular}
\end{center}
\mcaption{Top: Two vesicles discretized with $N=32$ points at several
time steps using near-singular integration.  Bottom left: The distance
between the two vesicles both with (solid) and without (dashed)
near-singular integration.  Bottom right: The error in area (blue) and
length (red) with (solid) and without (dashed) near-singular
integration.}{f:extensional:summary}
\end{figure}

We report the errors in area and length at $t=24$ in
Table~\ref{t:extensional:convergence}.  Also reported is the error in
the distance between the vesicles at $t=24$, for which we take the
``true'' distance from an overrefined simulation.  As before, we
achieve first-order convergence and there is little difference between
the accuracy when comparing explicit and semi-implicit inter-vesicle
interactions.  We also achieve second-order convergence, but it
requires a finer time step size when compared to the first-order
methods before the asymptotic rate is achieved.

In Table~\ref{t:extensional:convergence}, we also list the required CPU
time per time step relative to the smallest simulation.  When the
number of points is doubled from 32 to 64, the amount of work per time
step goes up by less than a factor of two.  The increase in CPU time
when $N$ is doubled from 64 to 128 is attributed to extra work required
by near-singular integration since, to simplify the implementation, we
only work with powers of 2.  Therefore, we upsample 64 to 256, but 128
is upsampled to 4096.  We also see that semi-implicit inter-vesicle
interactions are computationally not much more expensive.  This is
because a lot of work is spent precomputing the block-diagonal
preconditioner which is done once per time step for all the time
integrators.

\begin{table}[htps]
\begin{centering}
\begin{tabular}{ccccccc}
Integrator & $N$ &   $\Delta t$ & $e_{A}$ & $e_{L}$ & $e_{\mathrm{gap}}$ & CPU Time \\ \hline
Explicit (1) & $32$  & $0.04$ & $1.99e-4$ & $6.53e-4$ & $1.97e-1$ & --- \\
Explicit (1) & $64$  & $0.02$ & $6.34e-5$ & $3.42e-4$ & $3.73e-3$ & $1.15$ \\
Explicit (1) & $128$ & $0.01$ & $3.13e-5$ & $1.78e-4$ & $1.16e-3$ & $2.83$ \\
\hline
Implicit (1) & $32$  & $0.04$ & $2.12e-4$ & $6.64e-4$ & $1.92e-1$ & $2.65$ \\
Implicit (1) & $64$  & $0.02$ & $6.43e-5$ & $3.46e-4$ & $2.23e-3$ & $2.72$ \\
Implicit (1) & $128$ & $0.01$ & $3.15e-5$ & $1.79e-4$ & $4.35e-4$ & $7.40$ \\
\hline
Implicit (2) & $32$  & $0.04$ & $3.26e-4$ & $1.38e-4$ & $3.96e-2$ & $2.43$ \\
Implicit (2) & $64$  & $0.02$ & $4.14e-6$ & $4.45e-5$ & $8.32e-2$ & $2.89$ \\
Implicit (2) & $128$ & $0.01$ & $5.08e-7$ & $1.27e-5$ & $2.69e-2$ & $8.13$ \\
Implicit (2) & $128$ & $0.005$& $1.02e-7$ & $3.13e-6$ & $1.27e-4$ & $7.70$ \\
\end{tabular}
\mcaption{The errors in area, length, and the error in the gap size at
$t=24$, and the CPU time per time step for the extensional flow.  The
parenthetic values are the order of the time integrator.  The error in
the gap size is found by computing a ``true'' solution with $N=256$ points per
vesicle and a time step size so that the error in area is $9.0e-9$ and
in length is $7.3e-11$.  The CPU times are relative to the simulation
with $N=32$, $\Delta t = 0.04$ which took approximately
$2.62e-1$ seconds per time step.}{t:extensional:convergence}
\end{centering}
\end{table}

\subsection{Stenosis}

We consider a single vesicle of reduced area $0.65$ in a constricted
tube.  The solid wall is parameterized as ${\mathbf{x}}(\theta) =
10r(\theta)\cos(\theta)$ and $y(\theta) =
3r(\theta)\sin(\theta)\eta(\theta)$, where $r(\theta) =
(\cos(\theta)^{8} + \sin(\theta)^{8})^{-1/8}$, and
\begin{align*}
  \eta(\theta) = \left\{
    \begin{array}{cl}
      \frac{1-0.6\cos(x)}{1.6} & |x(\theta)| \leq \pi, \\
      1 & |x(\theta)| > \pi.
    \end{array}
  \right.
\end{align*}
For this geometry, the error of the trapezoid rule is $\mathcal{O}(10^{-4})$
with $N_{\mathrm{wall}}=256$ points and is $\mathcal{O}(10^{-6})$ with
$N_{\mathrm{wall}}=512$ points.  We impose a Gaussian-shaped boundary
condition at the intake and outake of the tube which is plotted along
with the geometry in Figure~\ref{f:stenosis:snaps}.  The vesicle's
initial height is $2.3$ times larger than the size of the constriction.
We report a convergence study for second-order time stepping in
Table~\ref{t:stenosis:errors} and we exceed the expected rate of
convergence.  We note that for this example, $N_{\mathrm{wall}}$ can not
be less than 256.  If it is, the Lagrange interpolation points required
by near-singular integration for wall-to-vesicle interactions are on
both sides of the solid wall.  Since the double-layer potential has a
jump across the solid wall, this creates a discontinuity in the velocity
field and the method becomes unstable.

This example is ideally situated for adaptive time stepping.  As the
vesicle passes through the constriction, large jumps in the errors are
committed.  In the simulation with $\Delta t = 5.0e-3$, the error in
length jumps by four orders of magnitude, and the error in area jumps
by three orders of magnitude near $t=3$.  These jumps are less severe
in the simulation with $\Delta t = 2.5e-3$, which has a jump of one
order of magnitude in length and two orders of magnitude in area.

\begin{figure}[htps]
\centering
\begin{tabular}{ccc}
\input{stenosisSnapsTime1.tikz} &
\input{stenosisSnapsTime2.tikz} &
\input{stenosisSnapsTime3.tikz} \\
\input{stenosisSnapsTime4.tikz} &
\input{stenosisSnapsTime5.tikz} &
\input{stenosisSnapsTime6.tikz} \\
\input{stenosisSnapsTime7.tikz} &
\input{stenosisSnapsTime8.tikz} &
\input{stenosisSnapsTime9.tikz}
\end{tabular}
\mcaption{Snapshots of a single vesicle in a constricted tube.  The
boundary condition on the solid wall (plotted in blue) at the inlet and
outlet is Gaussian-shaped, and on the remainder of the solid wall is
zero.  The vesicle is initially $2.3$ times larger than the size of the
constriction.  Results from~\cite{rah:vee:bir} were only able to
simulate vesicles that were $1.12$ times larger than the
constriction.}{f:stenosis:snaps}
\end{figure}

\begin{table}[htps]
\begin{centering}
  \begin{tabular}{cccccc}
    $N$ & $N_{\mathrm{wall}}$ & $\Delta t$ &
      $e_{A}$ & $e_{L}$ \\ \hline
    128 & 256 & $5.0e-3$ & $6.55e-4$ & $1.17e-2$ \\
    256 & 512 & $2.5e-3$ & $9.59e-7$ & $1.05e-6$
  \end{tabular}
  \mcaption{The errors in area and length at $t=8$ for the second-order
  stenosis flow.}{t:stenosis:errors}
  \end{centering}
\end{table}

We also report the minimum distance, and the time when the minimum
distance occurs, between the vesicle and the solid wall in
Table~\ref{t:stenosis:gapSize}.  In order to focus the error on the time
step size, we take a fine spatial resolution of $N=256$ and
$N_{\mathrm{wall}}=512$.  We report the minimum over time of the
distance
\begin{align*}
  d(\gamma,\Gamma) = \inf_{{\mathbf{x}} \in \gamma, {\mathbf{y}} \in \Gamma} \|{\mathbf{x}} - {\mathbf{y}}\|
\end{align*}
rather than the error since an exact solution is not available.  The
time integrators predict different locations, $t_{\min}$, where the
minimum occurs, but they both are converging.  For the first-order
results, if we compare the minimum distances with the distance of the
most accurate second-order simulation at $t=5.26$, which is $9.6369e-2$,
then we see that first-order convergence has been achieved.  Without a
much more accurate numerical solution, we are unable to compute
convergence rates for the second-order results.  However, it appears
that we have resolved the distance up to at least three digits of
accuracy.

\begin{table}[htps]
  \begin{centering}
  \begin{tabular}{c|cc|cc}
  & \multicolumn{2}{c|}{First-order} &
    \multicolumn{2}{c}{Second-order} \\
    $\Delta t$ & $t_{\min}$ & $d(\Gamma,\gamma)$ &
      $t_{\min}$ & $d(\Gamma,\gamma)$  \\ \hline
    $1.0e-2$ & $4.37$ & $1.0318e-1$ & $4.34$ & $1.0196e-1$ \\
    $5.0e-3$ & $4.41$ & $1.0262e-1$ & $4.39$ & $9.8055e-2$ \\
    $2.5e-3$ & $4.42$ & $1.0143e-1$ & $5.77$ & $9.5809e-2$ \\
    $1.3e-3$ & $5.26$ & $9.9541e-2$ & $5.77$ & $9.5804e-2$ \\
    $6.3e-4$ & $5.26$ & $9.7938e-2$ & $5.77$ & $9.5803e-2$ \\
    $3.1e-4$ & $5.26$ & $9.7137e-2$ & $5.77$ & $9.5802e-2$
  \end{tabular}
 \mcaption{The minimum distance, $d(\Gamma,\gamma)$, between the vesicle
 and the solid wall, and the time, $t_{\min}$, when this distance is
 achieved.  Rather than forming an overrefined ''true" solution, we
 report the actual distances.}{t:stenosis:gapSize}
 \end{centering}
\end{table}

\subsection{Couette Apparatus}
\label{s:couette}
Here we consider two simulations of a Couette apparatus with different
volume fractions.  We first consider a Couette apparatus where the
inner boundary, which is slightly off-centered, is rotating with
constant angular velocity and the outer boundary is fixed.  We randomly
place 42 vesicles of reduced area 0.65 inside the apparatus which
corresponds to a volume fraction of 20\%.  Each vesicle is discretized
with 64 points, each solid wall with 64 points, the time step size is
$0.01$, and the time horizon is $T=30$.  Snapshots of the simulation
are in the top plots of Figure~\ref{f:couetteLow}.  Using first-order
time stepping with explicit inter-vesicle interactions, the error in
area is $1.70e-2$ and the error in length is $7.94e-2$ at $t=30$.
Semi-implicit inter-vesicle interactions are not required for this time
step size.  The simulation took a little less than 6 hours, or about
7.1 seconds per time step.

We repeat this experiment with second-order time stepping which
requires semi-implicit inter-vesicle interactions.  The simulation took
about 83 hours, or about 100 seconds per time step.  The errors in area
and length at $t=30$ are $2.95e-3$ and $7.13e-4$.  To achieve a similar
accuracy with first-order time stepping, a time step size 11 times
smaller would be required which would require more than 83 hours of CPU
time due to the double-layer potentials.  Moreover, the errors in area
and length for second-order time stepping have plateaued while for
first-order time stepping, they continue to grow (bottom plots of
Figure~\ref{f:couetteLow}).  Again, this means that with second-order
time stepping, we can take much longer time horizons without requiring
a smaller time step size.

\begin{figure}[htps]
\begin{center}
  \begin{tabular}{c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}}
    \input{couetteSnapsTime1.tikz} &
    \input{couetteSnapsTime2.tikz} &
    \input{couetteSnapsTime3.tikz} &
    \input{couetteSnapsTime4.tikz} &
    \input{couetteSnapsTime5.tikz} &
    \input{couetteSnapsTime6.tikz}
  \end{tabular}
  \begin{tabular}{cc}
    \input{couetteErrorsOrder1.tikz} &
    \input{couetteErrorsOrder2.tikz}
  \end{tabular}
\end{center}
\mcaption{The top plots show snapshots of a Couette flow with 42
vesicles and a volume fraction of 20\%.  One vesicle is coloured in blue
to view its time history.  The outer boundary is fixed while the inner
boundary completes one revolution every 10 time units.  The bottom plots
show the errors in area and length for the first- and second-order time
integrators.  Notice that the errors for first-order method  continue to
grow whereas the errors for the second-order method
plateau.}{f:couetteLow}
\end{figure}

We now consider a Couette apparatus with the inner boundary exactly in
the middle of the outer boundary.  We use 150 vesicles of reduced area
0.65 resulting in a volume fraction of approximately $40.5\%$.  With
this higher concentration, a finer discretization of the solid walls is
required to resolve its density function ${\boldsymbol\eta}$.  Therefore, we have
discretized both boundaries with 256 points.  We also used our
collision detection algorithm to alert us if vesicles had crossed.
With $\Delta t = 0.01$ and explicit inter-vesicle interactions, the
first-order method results in crossing vesicles.  Therefore, we report
the results using first-order semi-implicit inter-vesicle interactions
which does not result in vesicles crossing.  Snapshots are illustrated
in Figure~\ref{f:couetteHigh}.  On average, each time step required 27
GMRES iterations. The total simulation took about 10 days, but this
number will drop dramatically once the double-layer potential is
implemented with a fast summation method.  At $t=10$, the error in area
is $9.62e-3$ and the error in length is $1.92e-2$.  We also ran this
simulation with second-order time stepping and achieved an error in
area of $9.39e-4$ and an error in length of $4.19e-4$.
\begin{figure}[htps]
\begin{center}
  \begin{tabular}{c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}}
    \input{couetteDenseSnapsTime1.tikz} &
    \input{couetteDenseSnapsTime2.tikz} &
    \input{couetteDenseSnapsTime3.tikz} &
    \input{couetteDenseSnapsTime4.tikz} &
    \input{couetteDenseSnapsTime5.tikz} &
    \input{couetteDenseSnapsTime6.tikz}
  \end{tabular}
\end{center}
\mcaption{Snapshots of a Couette flow with 150 vesicles.  One vesicle
is coloured in blue to view its time history.  The inner boundary
completes one full revolution every 10 time units.}{f:couetteHigh}
\end{figure}