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

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.


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 25 th 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 ε), 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 ε 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 ε 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.

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 augmented with the no-slip, no-penetration boundary condition u(X) = Ẋ for boundary points X, where overdot denotes time-derivative.Regularized stokeslet methods involve representing the flow field around a body B by an integral of the form, where S ε jk 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 ε > 0. The most widely-studied example [5] is for 3D flow and takes the form, The limiting form of this kernel is the classical stokeslet or Oseen tensor, Regularized stokeslet methods are implemented numerically by imposing equation eq. ( 2) for x = X ∈ B together with the condition u(X) = Ẋ (collocation), followed by discretization of the unknown traction f (y) and the numerical quadrature.The original (Nyström) discretization of eq. ( 2), by Cortez et al. [5], takes the form, where {X [1], . . ., X[Q] ∈ B} is a set of quadrature points, and the discretized force at X[q] is written as ).This method has the major advantage of implementational simplicity, a property which has resulted in widespread adoption; caveats are that the regularization parameter ε and discretization size h must be chosen in proportion, and that the typical size of the linear system (3Q × 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, where {B 1 , . . ., B N } is a partitioning of the surface B (mesh) with centroids x[n] ∈ 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 + ε 2 ) −1/2 as r → 0), this traction discretization need not be as refined as the quadrature discretization in equation eq. ( 5), i.e. one can take N 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 F = {x [1], . . ., x[N ]} and a 'fine quadrature' set Q = {X[1], . . ., X[Q]}, with N < Q.The force at the quadrature points f (y)dS(X[q]) is approximated by its value at the nearest force point via a nearest-neighbor projection f (X ), where ν[q, n] is a binary matrix with precisely one 1 in each row.The resulting linear system is then, where Smith [15] conducted an initial analysis of the error associated with the nearest-neighbor method eq. ( 7).In addition to a regularization error O(ε) 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(ε −2 h 2 f h q ) + O(h −1 f 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 ε, 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 (F) to the quadrature discretization (Q), denoted by δ.The error of the full problem is then discussed in terms of three distinct cases (detailed in fig.1): i) when δ > 0, i.e. the force and quadrature sets are disjoint, ii) when δ = 0, i.e. every force point is also a quadrature point, and iii) when the force and quadrature sets are non-disjoint, but δ = 0 for some points.

Analysis of the quadrature error for a single kernel evaluation
The principal challenge regarding numerical quadrature concerns evaluation where the kernel is rapidlyvarying, i.e.where |x − 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 ε → 0 because S ε jk (x, x) → ε −1 as ε → 0. This divergence also underlies the O(ε −2 h 2 f h q ) term in the nearestneighbor 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, Figure 1: Schematic detailing the characterisations for potential force (F) and quadrature (Q) discretization sets with δ, 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.
If any overlap x[n] = X[q] occurs, then clearly δ = 0. We also recall from [15] the definitions characterising the fineness of the force and quadrature discretizations, We note then, that we can characterise the choices of force (F) and quadrature (Q) discretizations for the nearest-neighbor method as one of three possibilities: These cases are illustrated in the schematic provided in fig. 1, with the classic (Nyström) discretization (which has δ = 0) included for comparison.
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 2π 0 1 0 K(r)rdrdθ where K(r) = (r 2 + ε 2 ) −1/2 , which captures the near-singular behaviour of the kernel S ε jk (x, y) for small |x − 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.

Previous analysis
As discussed by Smith [15], an estimate of quadrature error can be made using the mean value inequality we have for all constants a > ε, 1.In the region 0 r a, the bound 2. In the region a r, the bound Smith [15] used the above to split the quadrature into three regions, (i 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)

The disjoint case
To improve on this analysis we first address the case for which δ > 0. This entails that there is an inner region 0 r < δ 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, This error estimate is valid provided 0 < δ, ε 1 regardless of the relative sizes of ε and δ.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 δ r < h φ1 q , h φ1 q r < h φ2 q , . .., where φ 1 = 1 and φ 1 > φ 2 > . . .(this analysis deals with the case δ < h q ; if δ h q the error is no worse).The quadrature error for the first annulus is δ −2 h 2φ1+1 q and for the nth annulus is O(h 2(φn−φn−1)+1 q ).It is clear therefore that taking φ n − φ 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 φ n = 1 − (n − 1)/P for n = 1, . . ., P + 1.The error for the first annulus is O(δ −2 h 3 q ) = O((h q /δ) 2 h q ), and for the remaining annuli is O(h ).The total error over P annuli is therefore O((h q /δ) 2 h q ) + O(P h ).By taking increasingly large values of P the latter term approaches linear convergence.Therefore provided h q /δ = O (1), quadrature convergence is linear in h q and insensitive to ε.
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, for any integer P > 3, where δ is the minimum distance from the force discretization to the quadrature discretization, as defined in equation eq. ( 8).

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 λ > ε/ √ 2 so that the local maximum of |K (r)| appears inside the inner circle 0 r < λ.By similar arguments to the above (based on taking annuli of radius λ 1−(n−1)/P ) we then have the total error estimate, q , yields an error, As ε → 0, the dominant term in the above is O(ε −1 h 2 q ), which has the very advantageous property of being quadratic in h q , but an unwanted inverse dependence on ε.It is therefore clear for the contained case that we cannot expect to be able to reduce ε independently of h q .
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χ = b the relative error in the calculation of χ in terms of small deviations in the construction of matrix A, ∆A, is bounded by where cond (A) represents the condition number of A, and ∆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 ∆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.Each row of the matrix A consists of a diagonal entry which comes from evaluation of of S ε ij 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 In the case that Q becomes large then the change in the diagonal entries of A due to the evaluation of S ε ij 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 (ε −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 for a contained quadrature set, and, for a disjoint quadrature set, as ε → 0, where δ min and δ max are the smallest and largest of the distances δ (as calculated in eq. ( 8)) between each force and quadrature discretization sets, with δ max 1 (and assuming δ max /δ min = O (1)).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.

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 ε (r) = (r 2 + ε 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.

Quadrature convergence
The function K ε (r) is illustrated in fig.2a.Two types of quadrature method are illustrated: disjoint quadrature (fig.2b) for which the quadrature set does not include the origin, and contained quadrature (fig.2c) 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.3a) exhibits approximately linear convergence with h q and is insensitive to ε, as expected from equation eq. ( 12).The contained quadrature set performs very well for ε = 0.001, exhibiting approximately quadratic convergence, as expected from equation eq. ( 14).However, as also predicted, the absolute error shows an approximate ε −1 dependence, becoming highly inaccurate for ε 10 −5 .
Figure 5: Analysis of the condition number and diagonal entries of the matrix A. Panel (a) plots the condition number cond (A) against ε for a disjoint quadrature set, and shows the robustness of the condition number to changes as ε → 0. Panel (b) plots cond (A) − 1 against ε 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/ε 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 ε for large numbers of quadrature points Q with a contained quadrature set.

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.5a we plot the condition number of A against decreasing values of ε for a disjoint quadrature set and, as predicted, we see that the condition number plateaus rapidly as ε → 0. For a contained quadrature set we plot in fig.5b cond (A) − 1, against the same values of ε, where we see the approximate linear dependence on this quantity with ε (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.5c we plot the minimum increase in diagonal elements of A from 1/ against increasing numbers of quadrature points Q; here we clearly see that the diagonal terms are indeed bounded with A ii ≥ 1/ε.In fig.5d we plot the maximum row sum of off-diagonal entries C Q multiplied by ε.Here we see that εC Q grows slowly with Q, however for computationally practical values of Q the size of ε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.

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 δ −2 h 3 q and, if ε ≤ 10 −3 , is very robust to the choice of ε.This robustness in illustrated more clearly in table 1 which contains a selection of the values used to plot fig.6a.The contained case exhibits near linear convergence in ε −1 h 2 q for moderate values of , 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 ε with the relative error approaching 100% for ε ≤ 10 −5 .While we may naively expect the error to blow up for large ε −1 h 2 q , the limit of small ε 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 ε → 0.

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 (δ) from the force discretization (F) to the quadrature discretization (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.

The
When the force and quadrature discretizations are disjoint, the total quadrature error estimate for the nearest-neighbor regularized stokeslet method is for any integer constant P > 3, where ε 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 ε, with errors that are approximately linear in h q when δ ∼ h q .
When solving practical problems (of the form Aχ = b) with a disjoint quadrature set, the condition number of the matrix A is also robust to the choice of ε.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.6a.When the force and quadrature discretizations are contained, the total quadrature error estimate for the nearest-neighbor regularized stokeslet method is for any integer constant P > 3, where ε 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 ε.It is clear that in this case we are not able to reduce ε independently of h q .
When solving practical problems (of the form Aχ = b) with a contained quadrature set, the condition number of the matrix A grows linearly with ε.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.6b.

5.3
The general case, F ⊂ Q (δ ≥ 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 .

Discussion
The analyses were confirmed via numerical experiments: we have tested the convergence of numerical quadrature of the kernel K ε (r) (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 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.

Figure 2 :Figure 3 :
Figure 2: Set up for numerical quadrature experiments.(a) The near-singular kernel K ε (r) plotted for r < 1, with ε = 0.01.(b) of a disjoint quadrature set with h q = 1.33.(c) Depiction of an contained quadrature set with h q = 0.125.

Figure 4 :
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.

Figure 6 : 3 q
Figure 6: Absolute error in calculating the resistance tensor of a prolate spheroid undergoing body motion for (a) disjoint and (b) contained quadrature sets.In (a) we see convergence with decreasing δ −2 h 3 q with an insert showing a zoomed view emphasising the lack of ε dependence for ε < 10 −3 .Panel (b) shows convergence with decreasing ε −1 h 2 q .Each panel has been calculated for five values of the regularization parameter ε with fixed h f ≈ 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.6a).The error is calculated for five choices of the regularization parameter ε and increasing δ −2 h 3 q , with fixed h f ≈ 0.249.