# Sharp quadrature error bounds for the nearest-neighbor discretization of the
regularized stokeslet boundary integral equation

Meurig T. Gallagher1<EMAIL_ADDRESS><EMAIL_ADDRESS>School
of Mathematics, University of Birmingham,Birmingham. B15 2TT. UK.
<EMAIL_ADDRESS>Department of Mathematics, National Institute of
Technology, Rourkela, Odisha, India. 769 008.  Debajyoti Choudhuri2
<EMAIL_ADDRESS><EMAIL_ADDRESS>School of Mathematics,
University of Birmingham,Birmingham. B15 2TT. UK<EMAIL_ADDRESS>Department of Mathematics, National Institute of Technology, Rourkela, Odisha,
India. 769 008.  David J. Smith3<EMAIL_ADDRESS><EMAIL_ADDRESS>School of Mathematics, University of
Birmingham,Birmingham. B15 2TT. UK<EMAIL_ADDRESS>Department of
Mathematics, National Institute of Technology, Rourkela, Odisha, India. 769
008.

###### Abstract

The method of regularized stokeslets is a powerful numerical method to solve
the Stokes flow equations for problems in biological fluid mechanics. A recent
variation of this method incorporates a nearest-neighbor discretization to
improve accuracy and efficiency while maintaining the ease-of-implementation
of the original meshless method. This method contains three sources of
numerical error, the regularization error associated from using the
regularized form of the boundary integral equations (with parameter
$\varepsilon$), and two sources of discretization error associated with the
force and quadrature discretizations (with lengthscales $h_{f}$ and $h_{q}$).
A key issue to address is the quadrature error: initial work has not fully
explained observed numerical convergence phenomena. In the present manuscript
we construct sharp quadrature error bounds for the nearest-neighbor
discretisation, noting that the error for a single evaluation of the kernel
depends on the smallest distance ($\delta$) between these discretization sets.
The quadrature error bounds are described for two cases: with disjoint sets
($\delta>0$) being close to linear in $h_{q}$ and insensitive to
$\varepsilon$, and contained sets ($\delta=0$) being quadratic in $h_{q}$ with
inverse dependence on $\varepsilon$. The practical implications of these error
bounds are discussed with reference to the condition number of the matrix
system for the nearest-neighbor method, with the analysis revealing that the
condition number is insensitive to $\varepsilon$ for disjoint sets, and grows
linearly with $\varepsilon$ for contained sets. Error bounds for the general
case ($\delta\geq 0$) are revealed to be proportional to the sum of the errors
for each case.

##  1 Introduction

The development of numerical methods for the solution of Stokes flow has had
significant impact on the study of problems in biological fluid dynamics [1,
11, 12, 16, 17, 18] and vice-versa. While there have been many powerful
methods developed over the past few decades, one of the most effective and
accessible tools for solving such problems is the method of regularized
stokeslets, conceived of and developed by Cortez and colleagues [2, 4, 5, 6,
7, 8], and recently extended to incorporate the use of the fast multipole
method [13]. A key advantage of this method over previous offerings is the
meshless nature of the implementation, saving the significant investment of
time and effort it takes to generate a mesh (particularly when dealing with
complex biomolecular or cellular structures), potentially assisting with
automation for applications in image analysis. These methods have had
significant impact on a wide-ranging set of applications; a Google Scholar
search on 25th May 2018 with the term “regularized stokeslets” yielded $147$
results over the past year alone.

A new variation on the method of regularized stokeslets was recently proposed
by Smith [15], who uses a nearest-neighbor discretization of the regularized
stokeslet boundary integral equation to improve the accuracy and efficiency of
the classic Nyström discretization [5] while retaining the advantages of a
meshless method. The computational efficiency of this new method, together
with its extension to problems of locomotion in Stokes flow [9], enables the
study of previously computationally intractable problems such as improving the
detailed modelling of the embryonic node of mice and zebrafish to incorporate
Brownian and other effects, thus improving physiological accuracy.

Mathematical details of the nearest-neighbor discretization will be provided
in section 1.1, for now we note that the sources of error for the method are
threefold: the regularization error associated with using a regularized form
of the boundary integral equations (with regularization parameter
$\varepsilon$), and two sources of discretization error associated with the
approximating the the integral of the kernel at quadrature points with fine
discretization lengthscale $h_{q}$, and with approximating the forces with
coarser discretization lengthscale $h_{f}$. While the original paper of Smith
[15] provided an initial estimate of the error for the nearest-neighbor method
in terms of these two discretization lengthscales and regularization
parameter, these error bounds were noted not to be sharp as they did not fully
replicate the sensitivity to $\varepsilon$ suggested by the analysis (and seen
in the classic Nyström discretization). In the present work we not only
provide the detailed analysis for calculating sharp error bounds for the
nearest-neighbor discretization (see section 2), in doing this we uncover why
the original work did not see the dependence on the regularization parameter
$\varepsilon$ and detail the situations where this dependence exists. We then
also consider how these errors scale up in solving a practical problems
(section 3). Each of these analyses are then confirmed with numerical
experiments (in section 4). These results will provide clear guidance for the
best choices of discretization lengthscales and regularization parameter for
given computational scenarios.

###  1.1 Mathematical background

The dimensionless form of the Stokes flow equations, which describe the very
low Reynolds number fluid dynamics associated with sperm and cilia, is given
by

$-\bm{\nabla}p+\nabla^{2}\bm{u}=0,\quad\nabla\cdot\bm{u}=0,$ (1)

augmented with the no-slip, no-penetration boundary condition
$\bm{u}(\bm{X})=\dot{\bm{X}}$ for boundary points $\bm{X}$, where overdot
denotes time-derivative.

Regularized stokeslet methods involve representing the flow field around a
body $B$ by an integral of the form,

$u_{j}(\bm{x})=-\frac{1}{8\pi}\iint_{B}S_{jk}^{\varepsilon}(\bm{x},\bm{y})f_{k}(\bm{y})dS_{\bm{y}},$
(2)

where $S_{jk}^{\varepsilon}$ is the velocity part of the solution to the
Stokes flow equations eq. 1 driven by a smoothed point force in the
$k-$direction, with regularization parameter $\varepsilon>0$. The most widely-
studied example [5] is for 3D flow and takes the form,

$S_{ij}^{\varepsilon}(\bm{x},\bm{y})=\delta_{ij}\frac{|\bm{x}-\bm{y}|^{2}+2\varepsilon^{2}}{(|\bm{x}-\bm{y}|^{2}+\varepsilon^{2})^{3/2}}+\frac{(x_{i}-y_{i})(x_{j}-y_{j})}{(|\bm{x}-\bm{y}|^{2}+\varepsilon^{2})^{3/2}}.$
(3)

The limiting form of this kernel is the classical stokeslet or Oseen tensor,

$S_{jk}(\bm{x},\bm{y})=\frac{\delta_{jk}}{|\bm{x}-\bm{y}|}+\frac{(x_{j}-y_{j})(x_{k}-y_{k})}{|\bm{x}-\bm{y}|^{3}}.$
(4)

Regularized stokeslet methods are implemented numerically by imposing equation
eq. 2 for $\bm{x}=\bm{X}\in B$ together with the condition
$\bm{u}(\bm{X})=\dot{\bm{X}}$ (collocation), followed by discretization of the
unknown traction $\bm{f}(\bm{y})$ and the numerical quadrature. The original
(Nyström) discretization of eq. 2, by Cortez et al. [5], takes the form,

$u_{j}(\bm{X}[m])=\frac{1}{8\pi}\sum_{q=1}^{Q}S_{jk}^{\varepsilon}(\bm{X}[m],\bm{X}[q])F_{k}[q],$
(5)

where $\\{\bm{X}[1],\ldots,\bm{X}[Q]\in B\\}$ is a set of quadrature points,
and the discretized force at $\bm{X}[q]$ is written as
$F_{k}[q]=-f_{k}(\bm{X}[q])dS(\bm{X}[q])$. This method has the major advantage
of implementational simplicity, a property which has resulted in widespread
adoption; caveats are that the regularization parameter $\varepsilon$ and
discretization size $h$ must be chosen in proportion, and that the typical
size of the linear system ($3Q\times 3Q$) may be rather larger than would be
required by a classical boundary integral method to achieve converged results.
This scaling then limits the applications of the method when considering large
or complex problems.

The constant-panel boundary element discretization suggested by Smith et al.
[14] takes the form,

$u_{j}(\bm{x}[m])=\frac{1}{8\pi}\sum_{n=1}^{N}f_{k}[n]\iint_{B_{n}}S_{jk}^{\varepsilon}(\bm{x}[m],\bm{y})dS_{\bm{y}},$
(6)

where $\\{B_{1},\ldots,B_{N}\\}$ is a partitioning of the surface $B$ (mesh)
with centroids ${\bm{x}[n]\in B_{n}}$, and the discretization of the traction
on $B_{n}$ is denoted $f_{k}[n]$. Because the near-field of the regularized
stokeslet is rapidly-varying (resembling the function
$(r^{2}+\varepsilon^{2})^{-1/2}$ as $r\rightarrow 0$), this traction
discretization need not be as refined as the quadrature discretization in
equation eq. 5, i.e. one can take $N\ll Q$. The stokeslet integral in equation
eq. 6 is still evaluated numerically via quadrature, however in contrast to
equation eq. 5, the quadrature discretization does not affect the number of
degrees of freedom of the resulting linear system. The boundary element
discretization is therefore more efficient and accurate than the Nyström
method, however it has the disadvantage of requiring true mesh generation
rather than a simple list of surface points.

To attempt to combine the implementational simplicity of the Nyström method
with the efficiency and accuracy of the boundary element discretization, Smith
[15] proposed the use of a meshless nearest-neighbor method. Two
discretizations are generated, a ‘coarse force’ set
$\mathcal{F}=\\{\bm{x}[1],\ldots,\bm{x}[N]\\}$ and a ‘fine quadrature’ set
$\mathcal{Q}=\\{\bm{X}[1],\ldots,\bm{X}[Q]\\}$, with $N<Q$. The force at the
quadrature points $\bm{f}(\bm{y})dS(\bm{X}[q])$ is approximated by its value
at the nearest force point via a nearest-neighbor projection
$\bm{f}(\bm{X}[q])dS(\bm{X}[q])\approx\sum_{q=1}^{Q}\nu[q,n]\bm{f}(\bm{x}[n])dS(\bm{x}[n])$,
where $\nu[q,n]$ is a binary matrix with precisely one $1$ in each row. The
resulting linear system is then,

$u_{j}(\bm{x}[m])=\frac{1}{8\pi}\sum_{n=1}^{N}\left(\sum_{q=1}^{Q}S_{jk}^{\varepsilon}(\bm{x}[m],\bm{X}[q])\nu[q,n]\right)F_{k}[n],$
(7)

where $\bm{F}[n]=\bm{f}(\bm{x}[n])dS(\bm{x}[n])$.

Smith [15] conducted an initial analysis of the error associated with the
nearest-neighbor method eq. 7. In addition to a regularization error
$O(\varepsilon)$ found by Cortez et al. [5], the error associated with
discretization of the traction is $O(h_{f})$ (where $h_{f}$ characterises the
fineness of the force points), and the error associated with numerical
quadrature was estimated as
$O(\varepsilon^{-2}h_{f}^{2}h_{q})+O(h_{f}^{-1}h_{q})$ (where $h_{q}$
characterises the fineness of the quadrature points; formal definitions in eq.
9 and eq. 10 below). It was noted that this error bound was not sharp because
numerical experiments suggested that the error does not diverge for very small
values of $\varepsilon$, indeed the choice of $h_{f}$ as the lengthscale for
quadrature discretisation error was somewhat arbitrary. In the present
manuscript we will address this issue further. It will be shown that the
quadrature error for a single evaluation of the kernel depends on the shortest
distance from the force discretization ($\mathcal{F}$) to the quadrature
discretization ($\mathcal{Q}$), denoted by $\delta$. The error of the full
problem is then discussed in terms of three distinct cases (detailed in fig.
1): i) when $\delta>0$, i.e. the force and quadrature sets are disjoint, ii)
when $\delta=0$, i.e. every force point is also a quadrature point, and iii)
when the force and quadrature sets are non-disjoint, but $\delta=0$ for some
points.

##  2 Analysis of the quadrature error for a single kernel evaluation

The principal challenge regarding numerical quadrature concerns evaluation
where the kernel is rapidly-varying, i.e. where $|\bm{x}-\bm{y}|$ is ‘small’.
For the Nyström discretization the near-field part of the integral is
primarily evaluated in the sum eq. 5 when $q=m$. It is clear that this
evaluation is problematic as $\varepsilon\rightarrow 0$ because
$S_{jk}^{\varepsilon}(\bm{x},\bm{x})\rightarrow\varepsilon^{-1}$ as
$\varepsilon\rightarrow 0$. This divergence also underlies the
$O(\varepsilon^{-2}h_{f}^{2}h_{q})$ term in the nearest-neighbor error
estimate. A key advantage of the nearest-neighbor method however is that this
situation can be avoided by ensuring that the force and quadrature
discretizations are disjoint. We will denote the minimum distance between the
discretizations by,

$\delta=\min_{q=1,\ldots,Q}\min_{n=1,\ldots,N}|\bm{x}[n]-\bm{X}[q]|.$ (8)

If any overlap $\bm{x}[n]=\bm{X}[q]$ occurs, then clearly $\delta=0$. We also
recall from [15] the definitions characterising the fineness of the force and
quadrature discretizations,

$\displaystyle h_{f}$
$\displaystyle=\max_{m=1,\ldots,N}\min_{n=1,\ldots,N;n\not=m}|\bm{x}[m]-\bm{x}[n]|,$
(9) $\displaystyle h_{q}$
$\displaystyle=\max_{p=1,\ldots,Q}\min_{q=1,\ldots,Q;q\not=p}|\bm{X}[p]-\bm{X}[q]|.$
(10)

We note then, that we can characterise the choices of force ($\mathcal{F}$)
and quadrature ($\mathcal{Q}$) discretizations for the nearest-neighbor method
as one of three possibilities:

  1. 1.

Disjoint: $\mathcal{F}\cap\mathcal{Q}=\varnothing$ $\left(\delta>0\right)$;

  2. 2.

Contained: $\mathcal{F}\subset\mathcal{Q}$ $\left(\delta=0\right)$;

  3. 3.

General: $\mathcal{F}\cap\mathcal{Q}\neq\varnothing$ with
$\mathcal{F}\not\subset\mathcal{Q}$ $\left(\delta\geq 0\right)$.

These cases are illustrated in the schematic provided in fig. 1, with the
classic (Nyström) discretization (which has $\delta=0$) included for
comparison.

Figure 1: Schematic detailing the characterisations for potential force
($\mathcal{F}$) and quadrature ($\mathcal{Q}$) discretization sets with
$\delta$, denoting the minimum distance between the discretizations, defined
in equation eq. 8. Each choice is further shown as belonging to either the set
of nearest-neighbor or classic (Nyström) discretizations.

We will develop detailed analysis of the quadrature error of, in order, the
disjoint case (section 2.2), and the contained case (section 2.3). Each
analysis will be based on analysis of the error of approximation of
$\int_{0}^{2\pi}\int_{0}^{1}K(r)rdrd\theta$ where
$K(r)=(r^{2}+\varepsilon^{2})^{-1/2}$, which captures the near-singular
behaviour of the kernel $S_{jk}^{\varepsilon}(\bm{x},\bm{y})$ for small
$|\bm{x}-\bm{y}|$. In section 3 we will discuss how the quadrature errors
scale for practical problems and discuss the general case of mixed disjoint
and contained quadrature sets.

###  2.1 Previous analysis

As discussed by Smith [15], an estimate of quadrature error can be made using
the mean value inequality $|K(r)-K(0)|\leqslant M_{1}r$, where $M_{1}$ is a
bound on $|K^{\prime}(r)|$. Because
$K^{\prime}(r)=-r(r^{2}+\varepsilon^{2})^{-3/2}$, we have for all constants
$a>\varepsilon$,

  1. 1.

In the region $0\leqslant r\leqslant a$, the bound
$M_{1}=O(\varepsilon^{-2})$.

  2. 2.

In the region $a\leqslant r$, the bound $M_{1}=|K^{\prime}(a)|$.

Smith [15] used the above to split the quadrature into three regions, (i)
$0<r<h_{f}$, (ii) $h_{f}<r<h_{f}^{1/2}$ and (iii) $h_{f}^{1/2}<r<1$. The
quadrature errors can be estimated from the values of $M_{1}$, the area of the
region, and the quadrature spacing. The resulting error estimates, for each
region in turn, are then, (i) $O(\varepsilon^{-2}h_{f}^{2}h_{q})$, (ii)
$O(h_{f}^{-1}h_{q})$, (iii) $O(h_{f}^{-1}h_{q})$.

###  2.2 The disjoint case

To improve on this analysis we first address the case for which $\delta>0$.
This entails that there is an inner region $0\leqslant r<\delta$ which
contains no quadrature points, i.e. the region is neglected from the numerical
quadrature. The error associated with this neglect can be calculated as,

$L(\delta,\varepsilon)=\int_{0}^{2\pi}\int_{0}^{1}K(r)rdrd\theta=2\pi\left[(\delta^{2}+\varepsilon^{2})^{1/2}-\varepsilon\right]=O(\delta).$
(11)

This error estimate is valid provided $0<\delta,\varepsilon\ll 1$ regardless
of the relative sizes of $\varepsilon$ and $\delta$.

The remaining error can be calculated by a similar approach to section 2.1. To
achieve a sharp error estimate, we will consider a sequence of annuli
$\delta\leqslant r<h_{q}^{\phi_{1}}$, $h_{q}^{\phi_{1}}\leqslant
r<h_{q}^{\phi_{2}},\ldots$, where $\phi_{1}=1$ and $\phi_{1}>\phi_{2}>\ldots$
(this analysis deals with the case $\delta<h_{q}$; if $\delta\geqslant h_{q}$
the error is no worse). The quadrature error for the first annulus is
$\delta^{-2}h_{q}^{2\phi_{1}+1}$ and for the $n$th annulus is
$O(h_{q}^{2(\phi_{n}-\phi_{n-1})+1})$.

It is clear therefore that taking $\phi_{n}-\phi_{n-1}$ to be small and
negative will yield a close-to-optimal error estimate. For example, for any
fixed integer $P>3$ we may take $\phi_{n}=1-(n-1)/P$ for $n=1,\ldots,P+1$. The
error for the first annulus is
$O(\delta^{-2}h_{q}^{3})=O((h_{q}/\delta)^{2}h_{q})$, and for the remaining
annuli is $O(h_{q}^{1-2/P})$. The total error over $P$ annuli is therefore
$O((h_{q}/\delta)^{2}h_{q})+O(Ph_{q}^{1-2/P})$. By taking increasingly large
values of $P$ the latter term approaches linear convergence. Therefore
provided $h_{q}/\delta=O\left(1\right)$, quadrature convergence is linear in
$h_{q}$ and insensitive to $\varepsilon$.

In summary, the total error estimate (including regularization error, force
discretization error, and quadrature error) for the nearest-neighbor
regularized stokeslet method with disjoint discretizations is,

$E_{1}=O(\varepsilon)+O(h_{f})+O((h_{q}/\delta)^{2}h_{q})+O(Ph_{q}^{1-2/P}),$
(12)

for any integer $P>3$, where $\delta$ is the minimum distance from the force
discretization to the quadrature discretization, as defined in equation eq. 8.

###  2.3 The contained case

The analysis of section 2.1 is based on three regions parameterised by
$h_{f}$. The argument is in fact valid with $h_{f}$ replaced by any
lengthscale $\lambda>\varepsilon/\sqrt{2}$ so that the local maximum of
$|K^{\prime}(r)|$ appears inside the inner circle $0\leqslant r<\lambda$. By
similar arguments to the above (based on taking annuli of radius
$\lambda^{1-(n-1)/P}$) we then have the total error estimate,

$O(\varepsilon)+O(h_{f})+O(\varepsilon^{-2}\lambda^{2}h_{q})+O(P\lambda^{-2/P}h_{q}).$
(13)

Taking $\lambda=\varepsilon^{1/2}h_{q}^{1/2}$, yields an error,

$E_{2}=O(\varepsilon)+O(h_{f})+O(\varepsilon^{-1}h_{q}^{2})+O(P\varepsilon^{-1/P}h_{q}^{1-1/P}).$
(14)

As $\varepsilon\rightarrow 0$, the dominant term in the above is
$O(\varepsilon^{-1}h_{q}^{2})$, which has the very advantageous property of
being quadratic in $h_{q}$, but an unwanted inverse dependence on
$\varepsilon$. It is therefore clear for the contained case that we cannot
expect to be able to reduce $\varepsilon$ independently of $h_{q}$.

Figure 2: Set up for numerical quadrature experiments. (a) The near-singular
kernel $K^{\varepsilon}(r)$ plotted for $r<1$, with $\varepsilon=0.01$. (b)
Depiction of a disjoint quadrature set with $h_{q}=1.33$. (c) Depiction of an
contained quadrature set with $h_{q}=0.125$.

##  3 Practical implications of the quadrature error

Having developed the analysis to understand the quadrature error inherent in a
single evaluation of the kernel, it is of practical use to assess how this
error scales in a full application of the nearest-neighbor discretization
(when solving a resistance problem for example). When numerically constructing
and solving the matrix system $A\bm{\chi}=\bm{b}$ the relative error in the
calculation of $\bm{\chi}$ in terms of small deviations in the construction of
matrix $A$, $\Delta A$, is bounded by

${\Big{\|}}\frac{\Delta\bm{\chi}}{\bm{\chi}}{\Big{\|}}\leq{\big{\|}}A^{-1}{\big{\|}}{\big{\|}}A{\big{\|}}\frac{\big{\|}\Delta
A\big{\|}}{\big{\|}A\big{\|}}=\mathrm{cond}\left(A\right)\frac{\big{\|}\Delta
A\big{\|}}{\big{\|}A\big{\|}},$ (15)

where $\mathrm{cond}\left(A\right)$ represents the condition number of $A$,
and ${{\|}\Delta A{\|}}/{{\|}}A{\|}$ is the relative error in the numerical
construction of $A$. The analysis of section 2 provides the error estimates
for the size of $\|\Delta A\|/\|A\|$ for both disjoint (section 2.2) and
contained (section 2.3) quadrature sets; to build an understanding of the
error in a practical application of the nearest-neighbor method it thus
remains to understand how the condition number of the matrix $A$ behaves.

Figure 3: Convergence of numerical quadrature of the kernel $K^{\varepsilon}$
with number of points for four values of the regularization parameter
$\varepsilon$. (a) Disjoint quadrature set, (b) Contained quadrature set.

Each row of the matrix $A$ consists of a diagonal entry which comes from
evaluation of of $S_{ij}^{\varepsilon}$ eq. 3 at the force and associated
nearest-neighbor quadrature points. Consequently, for numerically tractable
numbers of quadrature points $Q$, the diagonal entries of $A$ have lower bound

$A_{ii}\geq\overline{A}\propto\begin{cases}\left(\delta^{2}+\varepsilon^{2}\right)^{-1/2}&\text{for
a disjoint quadrature set},\\\ \varepsilon^{-1}&\text{for a contained
quadrature set}.\end{cases}$ (16)

In the case that $Q$ becomes large then the change in the diagonal entries of
$A$ due to the evaluation of $S_{ij}^{\varepsilon}$ at many points will become
significant; we will explore this numerically in section 4, however for
practical densities of quadrature points this source of error is insignificant
compared to the dominant ($\varepsilon^{-1}$) term. Denoting the sum of the
off-diagonal elements (corresponding to a surface integral over a fixed area)
by $C_{Q}$, which will grow with increasing numbers of quadrature points $Q$,
we can apply the Gershgorin circle theorem [3] to show that all eigenvalues of
$A$ lie in a circle of radius $C_{Q}$ about the diagonal values in eq. 16. The
ratio between the largest and smallest eigenvalues (the condition number) is
therefore bounded by

$\frac{1/\varepsilon+C_{Q}}{1/\varepsilon-C_{Q}}\sim 1+2C_{Q}\varepsilon,$
(17)

for a contained quadrature set, and, for a disjoint quadrature set, as
$\varepsilon\rightarrow 0$,

$\frac{1/\delta_{\mathrm{min}}+C_{Q}}{1/\delta_{\mathrm{max}}-C_{Q}}\sim\frac{\delta_{\mathrm{max}}}{\delta_{\mathrm{min}}}\left(1+C_{Q}\left(\delta_{\mathrm{min}}+\delta_{\mathrm{max}}\right)\right),$
(18)

where $\delta_{\mathrm{min}}$ and $\delta_{\mathrm{max}}$ are the smallest and
largest of the distances $\delta$ (as calculated in eq. 8) between each force
and quadrature discretization sets, with $\delta_{\mathrm{max}}\ll 1$ (and
assuming $\delta_{\mathrm{max}}/\delta_{\mathrm{min}}=O\left(1\right)$).
Provided again that $Q$ is not too large (which will cause $C_{Q}$ to
correspondingly increase), is is clear that this method resolves a problem
that affects boundary element methods for Stokes flow: ensuring that the
condition number remains bounded as the size of the force elements approaches
zero.

The bounds that this analysis places on the condition number of the matrix $A$
are practically very useful when solving problems with the nearest-neighbor
discretization. There may be situations where it is desirable to discretize a
subject with both disjoint and contained quadrature sets; when considering a
biological swimmer, for example, it may be helpful to consider separately the
discretization of flagellum and body. We can thus consider the error of the
general case as being composed of the error from the disjoint case, $E_{1}$
eq. 12, plus the error from the contained case, $E_{2}$ eq. 14, multiplied by
the condition number of the matrix $A$.

##  4 Numerical experiments

We will now confirm the analysis through numerical experiments. Section 4.1
will consider the convergence of numerical quadrature of the function
$K^{\varepsilon}(r)=(r^{2}+\varepsilon^{2})^{-1/2}$ for a region including
$r=0$; section 4.2 will investigate the condition number of the matrix $A$,
and section 4.3 the resistance tensor, each owing to the problem of a prolate
spheroid undergoing rigid body motion.

###  4.1 Quadrature convergence

The function $K^{\varepsilon}(r)$ is illustrated in fig. 2. Two types of
quadrature method are illustrated: disjoint quadrature (fig. 2) for which the
quadrature set does not include the origin, and contained quadrature (fig. 2)
for which the quadrature set does include the origin. Numerical results with
these quadrature sets are shown in fig. 3. The disjoint quadrature set (fig.
3) exhibits approximately linear convergence with $h_{q}$ and is insensitive
to $\varepsilon$, as expected from equation eq. 12. The contained quadrature
set performs very well for $\varepsilon=0.001$, exhibiting approximately
quadratic convergence, as expected from equation eq. 14. However, as also
predicted, the absolute error shows an approximate $\varepsilon^{-1}$
dependence, becoming highly inaccurate for $\varepsilon\leqslant 10^{-5}$.

Figure 4: Sketch of the nearest-neighbor discretization of a prolate spheroid.
Here, the red dots show the force discretization, with a disjoint quadrature
set shown in light green.

###  4.2 Condition number

To assess the condition number analysis of section 3 on a relevant problem we
follow Smith [15] and construct the matrix $A$ owing to the resistance problem
of a prolate spheroid associated with rigid body motion. In fig. 5 we plot the
condition number of $A$ against decreasing values of $\varepsilon$ for a
disjoint quadrature set and, as predicted, we see that the condition number
plateaus rapidly as $\varepsilon\rightarrow 0$. For a contained quadrature set
we plot in fig. 5 $\mathrm{cond}\left(A\right)-1$, against the same values of
$\varepsilon$, where we see the approximate linear dependence on this quantity
with $\varepsilon$ (the slope in the figure is calculated as approximately
$1.02$). To assess the predictions regarding the diagonal entries of the
matrix $A$, in fig. 5 we plot the minimum increase in diagonal elements of $A$
from $1/\epsilon$ against increasing numbers of quadrature points $Q$; here we
clearly see that the diagonal terms are indeed bounded with $A_{ii}\geq
1/\varepsilon$. In fig. 5 we plot the maximum row sum of off-diagonal entries
$C_{Q}$ multiplied by $\varepsilon$. Here we see that $\varepsilon C_{Q}$
grows slowly with $Q$, however for computationally practical values of $Q$ the
size of $\varepsilon C_{Q}$ (and thus the condition number of $A$) remain
manageable. These numerical results agree with the analysis of section 3, and
give confidence to the calculated error bounds.

Figure 5: Analysis of the condition number and diagonal entries of the matrix
$A$. Panel (a) plots the condition number $\mathrm{cond}\left(A\right)$
against $\varepsilon$ for a disjoint quadrature set, and shows the robustness
of the condition number to changes as $\varepsilon\rightarrow 0$. Panel (b)
plots $\mathrm{cond}\left(A\right)-1$ against $\varepsilon$ for a contained
quadrature set, and shows approximate linear dependence with a slope of
$1.02$. Panel (c) shows the minimum deviation of the diagonal entries of $A$
away from $1/\varepsilon$ for large numbers of quadrature points $Q$ with a
contained quadrature set. Panel (d) shows the $C_{Q}$, the maximum sum of the
off-diagonal entries of $A$, multiplied by $\varepsilon$ for large numbers of
quadrature points $Q$ with a contained quadrature set.

###  4.3 Resistance problem

To assess the analysis of section 2 on a relevant problem we follow Smith [15]
and calculate the resistance tensor of a prolate spheroid, with an axis ratio
of $5$, associated with rigid body motion. This problem has the added benefit
of an analytical solution with which to compare [10]. A rendering of the force
and quadrature discretizations for this problem is shown in fig. 4. We again
test both disjoint and contained quadrature sets with results provided in fig.
6. The disjoint case again shows approximately linear convergence with
$\delta^{-2}h_{q}^{3}$ and, if $\varepsilon\leq 10^{-3}$, is very robust to
the choice of $\varepsilon$. This robustness in illustrated more clearly in
table 1 which contains a selection of the values used to plot fig. 6. The
contained case exhibits near linear convergence in $\varepsilon^{-1}h_{q}^{2}$
for moderate values of $\epsilon$, and the error collapses onto a single
curve. For this case we see, as predicted by the analysis in section 2, a
clear dependence on $\varepsilon$ with the relative error approaching $100\%$
for $\varepsilon\leq 10^{-5}$. While we may naively expect the error to blow
up for large $\varepsilon^{-1}h_{q}^{2}$, the limit of small $\varepsilon$ for
a contained quadrature set leads to a calculation of zero force, and thus a
zero resistance tensor, over the prolate spheroid. This results in a relative
error of $100\%$ in the limit $\varepsilon\rightarrow 0$.

Figure 6: Absolute error in calculating the resistance tensor of a prolate spheroid undergoing rigid body motion for (a) disjoint and (b) contained quadrature sets. In (a) we see convergence with decreasing $\delta^{-2}h_{q}^{3}$ with an insert showing a zoomed view emphasising the lack of $\varepsilon$ dependence for $\varepsilon<10^{-3}$. Panel (b) shows convergence with decreasing $\varepsilon^{-1}h_{q}^{2}$. Each panel has been calculated for five values of the regularization parameter $\varepsilon$ with fixed $h_{f}\approx 0.249$. Table 1: Error in calculating the resistance tensor of a prolate spheroid undergoing rigid body motion for a disjoint quadrature set (as depicted in fig. 6). The error is calculated for five choices of the regularization parameter $\varepsilon$ and increasing $\delta^{-2}h_{q}^{3}$, with fixed $h_{f}\approx 0.249$. | $\delta^{-2}h_{q}^{3}$
---|---
$\varepsilon$ | $2.282$ | $3.355$ | $4.805$ | $5.975$ | $9.101$ | $13.75$ | $23.40$
$10^{-2}$ | $1.603$ | $1.619$ | $1.929$ | $2.712$ | $3.666$ | $10.40$ | $20.66$
$10^{-3}$ | $0.4600$ | $0.8547$ | $1.560$ | $2.072$ | $3.699$ | $10.55$ | $27.35$
$10^{-4}$ | $0.4657$ | $0.8602$ | $1.564$ | $2.081$ | $3.702$ | $10.55$ | $27.52$
$10^{-5}$ | $0.4658$ | $0.8603$ | $1.564$ | $2.081$ | $3.702$ | $10.55$ | $27.52$
$10^{-6}$ | $0.4658$ | $0.8603$ | $1.564$ | $2.081$ | $3.702$ | $10.55$ | $27.53$

##  5 Conclusions

This paper has calculated sharp quadrature error bounds for the nearest-
neighbor regularized stokeslet discretization. We have shown that this error
depends on the shortest distance ($\delta$) from the force discretization
($\mathcal{F}$) to the quadrature discretization ($\mathcal{Q}$), and that the
behaviour of the quadrature error can be characterised by two discrete cases.
The total error in solving a Stokes flow problem using the nearest-neighbor
discretization is either described by one of these cases, or by a general,
mixed, case. We will now detail the characteristics of each of these cases in
turn.

###  5.1 The disjoint case, $\mathcal{F}\cap\mathcal{Q}=\varnothing$
$\left(\delta>0\right)$

When the force and quadrature discretizations are disjoint, the total
quadrature error estimate for the nearest-neighbor regularized stokeslet
method is

$E_{1}=O\left(\varepsilon\right)+O\left({h_{f}}\right)+O\left({\delta^{-2}h_{q}^{3}}\right)+O\left({Ph_{q}^{1-2/P}}\right),$

for any integer constant $P>3$, where $\varepsilon\ll 1$ is the stokeslet
regularization parameter, and $h_{f}$ and $h_{q}$ are given by equations eq. 9
and eq. 10 respectively. In this case the quadrature error is very robust to
the choice of $\varepsilon$, with errors that are approximately linear in
$h_{q}$ when $\delta\sim h_{q}$.

When solving practical problems (of the form $A\bm{\chi}=\bm{b}$) with a
disjoint quadrature set, the condition number of the matrix $A$ is also robust
to the choice of $\varepsilon$. This ensures that, providing $Q$ is not too
large, the error in solving a practical problem should remain $E_{1}$; this
error result has been validated through solving the resistance problem of a
prolate spheroid undergoing rigid body motion, displayed in fig. 6.

###  5.2 The contained case, $\mathcal{F}\subset\mathcal{Q}$
$\left(\delta=0\right)$

When the force and quadrature discretizations are contained, the total
quadrature error estimate for the nearest-neighbor regularized stokeslet
method is

$E_{2}=O\left(\varepsilon\right)+O\left({h_{f}}\right)+O\left({\varepsilon^{-1}h_{q}^{2}}\right)+O\left({P\varepsilon^{-1/P}h_{q}^{1-1/P}}\right),$

for any integer constant $P>3$, where $\varepsilon\ll 1$ is the stokeslet
regularization parameter, and $h_{f}$ and $h_{q}$ are given by equations eq. 9
and eq. 10 respectively. In this case the quadrature error is approximately
quadratic in $h_{q}$, but also has an unwanted inverse dependence of
$\varepsilon$. It is clear that in this case we are not able to reduce
$\varepsilon$ independently of $h_{q}$.

When solving practical problems (of the form $A\bm{\chi}=\bm{b}$) with a
contained quadrature set, the condition number of the matrix $A$ grows
linearly with $\varepsilon$. This ensures that, providing $Q$ is not too
large, the error in solving a practical problem should remain $E_{2}$; this
error result has been validated through solving the resistance problem of a
prolate spheroid undergoing rigid body motion, displayed in fig. 6.

###  5.3 The general case, $\mathcal{F}\not\subset\mathcal{Q}$ ($\delta\geq
0$)

We have provided analyses for each of the disjoint and contained cases. It is
conceivable that for practical purposes one may wish to use separate
discretizations for different elements of a problem (for example discretizing
a swimmer’s body and flagella differently) inducing the combination of both
disjoint and contained quadrature sets. The analysis of the condition number
of the matrix system indicates that the total error for the solution of
practical problems with general force and quadrature sets should be
proportional to the sum $E_{1}+E_{2}$.

###  5.4 Discussion

The analyses were confirmed via numerical experiments: we have tested the
convergence of numerical quadrature of the kernel
$K^{\varepsilon}\left(r\right)$ (section 4.1), assessed the change in
condition number of the matrix system (section 4.2), and tested the
calculation of the resistance tensor of a prolate spheroid undergoing rigid
body motion (section 4.3). Each of these numerical experiments closely
replicated the predictions of the analysis.

The analysis contained within the present work provides useful insight into
the error inherent in using the nearest-neighbor discretization; we believe
that the application of this analysis when choosing the parameters of the
method for a given problem of interest will be valuable in ensuring that the
desired convergence criteria are met.

The nearest-neighbor discretization For

The approach of the nearest-neighbor method reduces degrees of freedom while
retaining near field accuracy. For better scaling to large problems involving
many far field evaluations it may be interesting to explore whether fast
multipole implementations can be integrated into the nearest-neighbor method.
If it is possible for such adaptations to be made while keeping the ease-of-
implementation and simplicity of the present method then this will surely be
valuable.

One area where the nearest-neighbor discretisation may prove useful is in the
simulation of biological microswimmers: the meshless nature of this method
leaves open the possibility for automated swimmer generation from the analysis
of experimental imaging data, an option which is far from straightforward for
methods which require the generation of a true mesh.

## Acknowledgments

This work was supported by Engineering and Physical Sciences Research Council
award EP/N021096/1.

## References

  * [1] L. a. Klinteberg and A.-K. Tornberg, A fast integral equation method for solid particles in viscous flow using quadrature by expansion, J. Comp. Phys., 326 (2016), pp. 420–445.
  * [2] J. Ainley, S. Durkin, R. Embid, P. Boindala, and R. Cortez, The method of images for regularized Stokeslets, J. Comput. Phys., 227 (2008), pp. 4600–4616.
  * [3] H. E. Bell, Gershgorin’s theorem and the zeros of polynomials, Am. Math. Mon., 72 (1965), pp. 292–295.
  * [4] R. Cortez, The method of regularized Stokeslets, SIAM J. Sci. Comput., 23 (2001), pp. 1204–1225.
  * [5] R. Cortez, L. Fauci, and A. Medovikov, The method of regularized Stokeslets in three dimensions: Analysis, validation, and application to helical swimming, Phys. Fluids, 17 (2005), p. 031504.
  * [6] R. Cortez and F. Hoffmann, A fast numerical method for computing doubly-periodic regularized Stokes flow in 3D, J. Comput. Phys., 258 (2014), pp. 1–14.
  * [7] R. Cortez and M. Nicholas, Slender body theory for Stokes flows with regularized forces, Comm. Appl. Math. Comput. Sci., 7 (2012), pp. 33–62.
  * [8] R. Cortez and D. Varela, A general system of images for regularized Stokeslets and other elements near a plane wall, J. Comput. Phys., 285 (2015), pp. 41–54.
  * [9] M. T. Gallagher and D. J. Smith, Meshfree and efficient modelling of swimming cells, Phys. Rev. Fluids, In press. (2018).
  * [10] S. Kim and S. J. Karrila, Microhydrodynamics: principles and selected applications, Courier Corporation, 2013.
  * [11] C. Pozrikidis, Boundary integral and singularity methods for linearized viscous flow, Cambridge Univ Press, 1992.
  * [12] C. Pozrikidis, A Practical Guide to Boundary Element Methods with the Software Library BEMLIB, CRC, 2002.
  * [13] M. W. Rostami and S. D. Olson, Kernel-independent fast multipole method within the framework of regularized stokeslets, J. Fluid. Struct., 67 (2016), pp. 60–84.
  * [14] D. Smith, A boundary element regularized Stokeslet method applied to cilia-and flagella-driven flow, Proc. R. Soc. Lond. A, 465 (2009), pp. 3605–3626.
  * [15] D. Smith, A nearest-neighbour discretisation of the regularized stokeslet boundary integral equation, J. Comput. Phys., 358 (2018), pp. 88–102.
  * [16] S. Veerapaneni, D. Gueyffier, D. Zorin, and G. Biros, A boundary integral method for simulating the dynamics of inextensible vesicles suspended in a viscous fluid in 2D, J. Comp. Phys., 228 (2009), pp. 2334–2353.
  * [17] S. Veerapaneni, A. Rahimian, G. Biros, and D. Zorin, A fast algorithm for simulating vesicle flows in three dimensions, J. Comp. Phys., 230 (2011), pp. 5610–5634.
  * [18] A. Zinchenko and R. Davis, An efficient algorithm for hydrodynamical interaction of many deformable drops, J. Comp. Phys., 157 (2000), pp. 539–587.