{'fraction_non_alphanumeric': 0.04987123750449977, 'fraction_numerical': 0.020878907872511284, 'mean_word_length': 4.080039386692924, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Regularized Single and Double Layer Integrals in 3D Stokes Flow

We present a numerical method for computing the single layer (Stokeslet) and double layer (stresslet) integrals in Stokes flow. The method applies to smooth, closed surfaces in three dimensions, and achieves high accuracy both on and near the surface. The singular Stokeslet and stresslet kernels are regularized and, for the nearly singular case, corrections are added to reduce the regularization error. These corrections are derived analytically for both the Stokeslet and the stresslet using local asymptotic analysis. For the case of evaluating the integrals on the surface, as needed when solving integral equations, we design high order regularizations for both kernels that do not require corrections. This approach is direct in that it does not require grid refinement or special quadrature near the singularity, and therefore does not increase the computational complexity of the overall algorithm. Numerical tests demonstrate the uniform convergence rates for several surfaces in both the singular and near singular cases, as well as the importance of corrections when two surfaces are close to each other.


Introduction
Stokes flows are of relevance in many practical problems where the length scales are very small, the fluid is viscous, or the velocity is very small (i.e., creeping flows), all resulting in a small Reynolds number. In dimensionless form, the incompressible Stokes equations are where p is the pressure and u is the flow velocity. The Stokeslet and stresslet are the primary fundamental solutions for the velocity u: where δ ij is the Kronecker delta and i, j, k = 1, 2, 3 are Cartesian coordinates. When used in boundary integral methods, these lead to the single and double layer representations of Stokes flow, respectively: where n k are the components of the unit outward normal vector to the surface. The integral in (3a) is continuous across ∂Ω, and the integral in (3b) is discontinuous and has a jump of ∓4πq i (x 0 ) in the limit from either the interior or exterior of the domain. The importance of boundary integral equations in Stokes flow models is well recognized. In particular, they have been used extensively in interfacial dynamic simulations [21,29], such as vesicle flows [28], drop dynamics [34,15], and particle motion [22]. For such applications, the jump conditions across the interface are incorporated into the integral formulation naturally, the dimensionality of the problem is reduced, and high accuracy can be achieved on the boundary and for points away from the boundary. Evaluating the integrals accurately for points near the boundary, e.g., when two interfaces are close together, is the most difficult case and is an active area of research [25,16,8,19,33,9,14,4,10]. If the viscosities inside and outside are different, the interface velocity is found from an integral equation with single and double layers [20,21,23,25,31].
In dealing with the evaluation of nearly singular integrals, Ying, Biros and Zorin [33] proposed an interpolation procedure. The authors considered 3D elliptic problems, with domain boundaries given by overlapping patches parameterized using coordinate charts. In particular they computed double layer integrals for Stokes velocity as in (3b). Since the errors decay rapidly away from the surface, the values near the boundary were obtained by interpolating between the values at points on the surface and points that are sufficiently separated from the surface along the surface normals. This algorithm was adapted and optimized by Sorgentone and Tornberg [25] for close interactions of viscous drops with surface tension, where a spherical harmonics expansion was used to parameterize the surface. A quadrature by expansion method was developed by Klöckner et al. [16] and Barnett [3] for evaluation of Laplace and Helmholtz potentials, through local expansions. This method achieves exponential accuracy but requires upsampling the density on a finer grid. Af Klinteberg and Tornberg [1] applied the quadrature by expansion method to simulate spheroidal particles in periodic Stokes flow, using precomputations and the fast Ewald summation method for faster computations. A target-specific QBX method was developed by Siegel and Tornberg [24], where the same accuracy is achieved using fewer terms. Another approach was followed by Bruno and Kunyansky in [9] for surface scattering problems, where partitions of unity were used along with an analytical resolution of the singularity by a change to polar coordinates.
The method of Beale and coworkers [5,6,7] follows a different approach. The integrals are first regularized to remove the singularity. This introduces a regularization error, which is reduced by adding correction terms derived analytically using local expansions. Discretizing the integrals will introduce an additional error component, and discretization corrections can be derived as well to reduce this error. This work has led to further explorations for the Laplace's equation [6,8], the Stokes equations [26,19], and the Helmholtz equation [18]. In [26], the Stokes single layer integral, which has the Stokeslet as the kernel, was written using the gradient of the Laplacian Green's function as the kernel, and then regularization and discretization corrections were derived for the normal and tangential parts of this integral. The method of regularized Stokeslets [12,13] is also based on desingularizing the velocity, but the approach is slightly different. It can be viewed as replacing the entire problem by a regularized one by smoothing out the forces acting on the fluid particles, whereas we choose the regularization specifically to compute the surface integrals accurately. In [19], Nguyen and Cortez derived regularization corrections for the regularized Stokeslets with an orthogonal system of coordinates, while paying special attention to preserving the incompressibility condition.
In this paper we extend the framework described above to evaluate both the single and double layer integrals of Stokes flow in three dimensions on or near smooth, closed surfaces. We only treat the regularization component of the error, and our numerical results show that, with proper choice of the regularization parameter, high order uniform convergence of the integrals (3a), (3b) is achieved, uniformly with respect to y near the surface. For evaluation at points on the surface, as needed in solving integral equations, we design regularizations for the Stokeslet and stresslet kernels with high order accuracy, approaching O(h 5 ). Here h is the grid spacing chosen in three-space and used in coordinate planes for the discretization of the integrals, as explained in Section 2. For evaluation near the surface we use simpler regularization and derive corrections to achieve accuracy uniformly about O(h 3 ) in practice. With this approach, there are no parameters to fine-tune except the regularization parameter δ. Through experimentation with this parameter we find that δ/h = 3 for the case on the surface and δ/h = 2 near the surface are reliable choices. Because the correction formulas are precomputed analytically using local asymptotic analysis, the overall computational complexity of the algorithm does not increase. Also, no special gridding or quadrature is needed near the singularity, so that the spacing does not change with proximity to the boundary. Our method here is more direct than in [26] and treats the double layer or stresslet integral as well as the Stokeslet. As in [8], we use a quadrature rule for closed surfaces, introduced in [32], which works well for general surfaces without requiring coordinate charts.
The numerical method is summarized in Section 2, including the formulas for the regularization corrections used for points near the surface. These formulas are derived analytically in Sections 3.1 and 3.2. The higher order regularizations for evaluation on the surface are given in Section 3. 3

Numerical Method
At the heart of the numerical method presented here is the regularization of the singularities that develop in the kernels of (3a)-(3b) as r = |y − x| approaches 0. As a first step, we use subtraction to reduce the singularity in the double layer (3b), where x 0 is the boundary point closest to y, and we have applied the well known identity (see, for example, [20] sec. 2.1-2.3, or [21]) where χ(y) = 8π, 4π, 0 if y is inside, on, and outside the boundary, respectively. We then regularize the Stokeslet as follows where the smoothing factors are chosen so that lim ρ→∞ s 1 (ρ) = 1, lim ρ→∞ s 2 (ρ) = 1, s 1 (ρ) = O(ρ) and s 2 (ρ) = O(ρ 3 ) for ρ small, and s 1 (r/δ)/r, s 2 (r/δ)/r 3 are smooth as functions of y − x with r = |y − x| and fixed δ > 0. The regularized version of the stresslet is similar, whereq(x) = q(x) − q(x 0 ), and s 3 is chosen with lim ρ→∞ s 3 (ρ) = 1, s 3 (ρ) = O(ρ 5 ) for small ρ, and s 3 (r/δ)/r 5 smooth for δ > 0. There are a number of possible choices for the smoothing factors; several were explored in [19]. In this work, we will use where erf is the error function. This choice of regularization is simple; (s 1 − 1)/r decays rapidly in the far field, and s 2 , s 3 are derived by modifying the error function to get s 2 (ρ) = O(ρ 3 ) and s 3 (ρ) = O(ρ 5 ) for small ρ as mentioned above.
The regularization error, defined as u δ − u for the Stokeslet and w δ − w for the stresslet, is at best O(δ). Our approach to reducing this error is twofold. When the integrals are evaluated on the surface (e.g., when solving integral equations), we modify the smoothing functions (8a)-(8c) to readily achieve a higher accuracy of O(δ 5 ). This is explained in Section 3.3. When the evaluation point is near the surface, however, such a modification is not helpful, and we increase the accuracy by adding corrections. These corrections represent the dominant terms in the regularization error, and improve the accuracy to O(δ 3 ). Specifically, we compute the single layer and the double layer as with corrections, to be explained, Here, assuming that y is near the surface, x 0 is the surface point closest to y, so that y = x 0 + bn for some b, where n is the unit outward normal at x 0 . H is the mean curvature and λ = b/δ. The superscript nl denotes the normal part of a vector, a nl ≡ (a · n)n, and superscript tan denotes the tangential part, a tan ≡ a − a nl . For any local coordinate system x = x(α 1 , α 2 ) on the surface, we have tangent vectors T i = ∂x/∂α i , dual vectors T * i , i = 1, 2, such that T * i · T j = δ ij , and the metric tensor g ij = T i · T j . Note that a tan = (a · T 1 )T * 1 + (a · T 2 )T * 2 . The surface gradient of a scalar function f , the surface divergence of a vector function v, and the surface Laplacian of a scalar function f are defined as They are independent of coordinates. It can be checked that (11b) agrees with formula (9.41.1) in [2] for the surface divergence of a tangential vector field. Section 4.1 has some details on how these quantities are computed for a Monge parameterization. In (10a), (10b) I 1 , . . . , I 3b are certain integrals that occur in the derivations of Section 3. With our choice of s 1 , s 2 , s 3 , where erfc(x) = 1 − erf(x). If we chose (8a,8b,8c) differently, the expressions for the corrections would be the same, but I 1 , . . . , I 3b would be different.
Once the integrands are smoothed out, we discretize the integrals using a surprisingly simple quadrature method for closed surfaces introduced in [32] and explained in [8]. We choose an angle θ and define a partition of unity on the unit sphere, consisting of functions Fig. 1. Sets R 1 and R 2 are defined similarly. For a function f on the surface the integral is computed as The partition of unity functions ψ i are constructed from the C ∞ bump function b(r) = exp ((r 2 /(r 2 − 1)) for |r| < 1 and zero otherwise. The quadrature is effectively reduced to the trapezoidal rule without boundary. Thus for regular integrands the quadrature has arbitrarily high order accuracy, limited only by the degree of smoothness of the integrand and surface. The points in R i can be found by a line search since they are well separated; see [32] and [8].
Finally we discuss the accuracy of this method, relying on the analogy with the case of Laplace's equation treated in [6], [8]. Error estimates for the harmonic double layer potential were derived in [6] and extended to the single layer in [8]. The present work is similar except that discretization corrections were included in the earlier case. Based on the theory for the harmonic potentials, we expect that the error in the present method, for evaluation of (3a, 3b) at points near the surface, can be estimated uniformly as The first term represents the regularization error remaining after the correction (10a, 10b), and the second term is the discretization error. For evaluation on the surface, using the formulas of Sec. 3.3, the error estimates for the single and double layer integrals should be The improvement in the double layer results from the subtraction in equation (4). These estimates are discussed further in Sec. 3.4 of [8]. The accuracy depends critically on the relationship between δ and h. In this work we have taken δ/h = constant, such as 2 or 3, for simplicity. In the numerical examples reported in Section 4, we see errors about O(h 3 ) near the surface and significantly higher order on the surface, indicating that the regularization error is dominant. However, this pattern could not continue as h → 0 with δ/h fixed. In principle the discretization error can be controlled as h → 0 by increasing δ/h, e.g. δ = ch p for any chosen p < 1. Then the exponential in the error estimates decreases rapidly, and the second term, representing the discretization error, is dominated by the first term as h → 0. In this way convergence can be achieved which is O(h 3p ) near the surface and higher on the surface.

Regularization Corrections
In this section, we derive the regularization corrections for the Stokeslet and stresslet integrals, and discuss the special case of evaluating them on the surface with very high accuracy.

Stokeslet
Computing the regularization correction in the first part of the Stokeslet (6) is similar to the derivation for the Laplace's single layer potential in [8], but the second part needs extra care. First write the error as where we define φ(r) = s(r) − 1. For the part of the surface near y we will use a special parameterization x(α), with x(0) = 0 and y along the normal line from x(0), so that y = bn 0 for some b, where n 0 is the unit normal at x(0). With tangent vectors to the surface T i = ∂x/∂α i , i = 1, 2 and metric tensor g ij = T i · T j , we assume that, at α = 0, g ij = δ ij and ∂g ij /∂α k = 0, i, j, k = 1, 2. Also, rotating if necessary, T 1 , T 2 have the directions of principal curvature, and n 0 = T 1 × T 2 . We use expansions near α = 0, where In the expansion of F := [f (α)·(y−x(α))](y−x(α)), only the terms even in α will contribute: We choose a new parameter ξ to replace α, defined by r 2 = b 2 + |ξ| 2 and ξ i /|ξ| = α i /|α|, F even Note that F even 1 is F even where terms with ξ 1 ξ 2 were omitted, as they lead to 2π 0 cos(θ) sin(θ)dθ or 2π 0 cos 3 (θ) sin(θ)dθ type coefficients once we change ξ to polar coordinates, and thus contribute zero. The regularization error is now written as where m 1 includes the nonradial terms, where |∂α/∂ξ| = 1 + bµ + O(|ξ| 2 + b 2 ), and |T 1 × T 2 | = 1 + O(|ξ| 2 ). Combining the terms in m 1 and neglecting the terms odd in ξ (since they contribute 0 to the integral in ξ), we get m even Similarly for m 2 , Substituting (23) and (24) into (21), and changing to polar coordinates such that ξ = δζ, b = δλ, |ζ| = η, q = κ 1 cos 2 θ + κ 2 sin 2 θ, we get which simplifies to where H = (κ 1 + κ 2 )/2 is the mean curvature, and This gives the correction expressed in a special coordinate system. To extend this to an arbitrary system, first we find from (15a)-(15b) that in our special coordinates at α = 0, we have ∂ i T i = κ i n 0 and ∂ i n = −κ i T i , i = 1, 2, where ∂ i is the partial derivative in α i . Using these we get at α = 0 Multiplying by −1, we see that the two terms on the right side of (28a) occur in (26b) and (26c). Similarly the two terms in (28b) occur in (26a) and (26c). We can combine terms to get, with the sum over i = 1, 2 implied, The dual basis in the tangent space is defined as the surface divergence at α = 0, according to (11b); we used the fact that ∂ i T i is normal. We can now write the regularization error as where f nl 0 and f tan 0 are the normal and tangential parts of f 0 , leading to the correction (10a).

Stresslet
The regularization error in the stresslet can be written as where as before,q(x) = q(x) − q(x 0 ). In what follows, we will write q forq and assume q(0) = 0 because of the subtraction. We compute the Taylor expansions of the nonradial parts: Next we compute the product of (32a) and (32b): Now, let In the expansion of this quantity, we keep only terms even in α, and we also neglect terms with α 1 α 2 , α 3 1 α 2 , and α 1 α 3 2 , as all of these will contribute 0 to the error. We then get the following: Next, we rewrite Q even in the new parameter ξ, and then write the regularization error (31) as where Substituting (37) into (36), and changing to polar coordinates as before, we get which simplifies to with H the mean curvature, as before, and The quantities in the square brackets in (39a) are the surface divergence times the normal, (∇ S ·q)n 0 , and the surface gradient ∇ S (q·n). These are similar to the Stokeslet terms except now we have q 0 = 0 because of the subtraction. We will show that the sum of terms with κ i in (39a) and (39b) is δλI 3a /2 times ∆ S [(q · n)n]. We compute the latter using the facts that q 0 = 0, and in our special coordinates ∆ S = Σ i ∂ 2 i and ∂ i n = −κ i T i , both at 0. We find and The latter matches (39a)-(39b), verifying our assertion. We note for later use that with sum over i, where tan means the tangential part. Since ∆ S q = Σ i q ii , the first term in brackets in (39c) is {∆ S q} tan . We will relate the second term to the surface gradient of the surface divergence of q. To identify the latter, we first note that where q nl = (q · n)n is the normal part of q. Also, for scalar F we have ∇ S F = Σ j (∂ j F )T * j . To find ∇ S (∇ S · q) at 0, we begin with with q i = ∂ i q etc. and sums over i, j implied. Then for the normal part we have where in the first step we have used the fact that n · T * i ≡ 0. Now combining (44)-(46) and using (43) we get In summary we have shown that the terms in brackets in (39c) equal We can now write a complete formula for , leading to the correction (10b),

Evaluation on the boundary
When evaluating the integrals on the boundary, we can remove the lower order terms in the regularization error and make it O(δ 5 ) by an appropriate choice of the smoothing function.
To find this function, we write s # (r) = s(r) + ars (r) + br 2 s (r), where constants a and b are chosen to make two moments equal to 0. In (25)-(27c) we now have λ = 0. For the first part of the Stokeslet, with φ 1 , the moment conditions are I 1 = 0 and a similar integral with η 3 in place of η. This gives us the function used in [8] for the Laplacian single layer, For the second part of the Stokeslet (see (25)-(27c) with φ 2 ), the I 2a term vanishes on the surface, and the moment conditions are I 2b = 0 and a similar integral with η 5 . The fifth order smoothing is A similar approach for the stresslet leads us to s # 3 (r) = erf(r) − 2 9 r(8r 6 − 36r 4 + 6r 2 + 9)e −r 2 / √ π.
Using these regularization functions for the case on the surface gives high order convergence without the need to compute corrections.

Numerical experiments 4.1 Monge parameterization
In the table below we describe some details on how the geometric quantities were computed. Much of this can be found in Appendix B of [8], but we summarize it here for completeness. Suppose the surface is given by φ(x 1 , x 2 , x 3 ) = 0, with φ > 0 outside, and the normal vector ∇φ is predominantly in the x 3 direction. Then the parameterization is given by , where z is the vertical coordinate on the surface.

Evaluation formula
Notes

Flow around a translating spheroid
We first define and compute the flow around a sphere (a = b) and a prolate spheroid (a > b), translating with velocity U . This allows us to test the single layer integral alone. The exact solutions for the flow are well known for both the sphere and the spheroid [11,17]. Taking U = (1, 0, 0), the velocity outside the spheroid is where e = c/a, c 2 = a 2 − b 2 , e r = (x 2 e 2 + x 3 e 3 )/r is the unit radial vector in the x 2 x 3 -plane, and e i , i = 1, 2, 3 are the unit basis vectors. The surface traction is given by (f 1 (x), 0, 0) with The expressions for the unit sphere can be found in the references mentioned above, or for example, in [26].
In this test, we compute the velocity (3a) given the surface traction in (58). First, we choose a grid size h and find the quadrature points as explained in Section 2. Second, to choose points near the surface we cover R 3 with a 3D grid of size h and select the points that are ≤ h distance outside the surface. For these points, we compute the velocity using (9a), (6), (10a), with s 1 , s 2 defined in (8a), (8b), and the regularization parameter δ chosen so that δ/h is a constant. The integral in (6) is discretized by the quadrature rule (13), evaluating the traction f at the quadrature points using (58). To compute the velocity on the surface, the same quadrature is used but without corrections and with s 1 , s 2 in (6) replaced by s # 1 , s # 2 from (50), (51). Then we compare the computed velocity with the exact value (54). We define the error at a single point as e(x) = |u computed (x) − u exact (x)|, where |·| is the vector's Euclidean norm. We then compute either the max or the L 2 norm of this error over the evaluation points. The L 2 norm is defined as e 2 = x e 2 (x)/n 1/2 , where n is the number of evaluation points. All errors reported are absolute errors; the largest velocity in magnitude is the translational velocity U , which was taken to be 1 for both the sphere and the spheroid. To demonstrate how δ/h should be chosen in practice, in our first test we keep h fixed and vary δ. For larger values of δ/h, the discretization error is small and the regularization error is dominant. For small δ/h, the discretization error will dominate, and is expected to be O(h). Figure 2 shows the errors for the unit sphere, where the grid size h = 1/64 was chosen and δ was varied. The left graph is the error at the quadrature points, where the higher regularization (50)-(51) was used and no corrections are necessary. For a spheroid, rather than a sphere, we find that the error increases with increasing δ/h. Since the integral on the surface is computed with very high accuracy in δ (Section 3.3), the regularization parameter should be chosen larger to preserve the high order of convergence. The value δ = 3h seems to suffice in all tests we performed. The graph on the right in Figure 2 is the error for points close to the surface. The graph demonstrates a typical behavior where δ/h > 1 should be chosen so that the regularization error dominates leading to a higher order of convergence in grid size h. The regularization δ = 2h seems sufficiently large to achieve O(h 3 ) convergence, as shown in our subsequent testing. We further demonstrate the predicted convergence rates by refining the grid while keeping δ/h fixed. Figure 3 shows the errors for the unit sphere for points on the surface (taken to be the quadrature points), and points that are ≤ h distance outside the surface. We repeat this test for the spheroid a = 1, b = 0.5 in Figure 4. The errors are larger for the spheroid due to the larger curvature and varied spacing: for h = 1/32 for example, there are 17070 quadrature points on the unit sphere, but only 6958 quadrature points on the spheroid. Both cases display the predicted order of convergence. As described earlier, choosing δ/h < 1 will result in slower convergence of O(h), while δ/h = 1 is in the intermediate regime where convergence can be observed as O(h) (as we see for the sphere) or faster (as it happens to be for the spheroid). A larger δ/h = 2 will typically result in the predicted O(h 3 ) convergence off the surface, as seen for both the sphere and the spheroid. On the surface, an even larger regularization is recommended, and we choose δ/h = 3. The errors decrease rapidly with observed convergence rates close to O(h 5 ).  Figure 3: Single layer integral, errors for the unit sphere: over quadrature points (left) and over points distance ≤ h outside the surface (right).  Figure 4: Single layer integral, errors for the spheroid a = 1, b = 0.5: over quadrature points (left) and over points distance ≤ h outside the surface (right).

Double layer identity
In order to now test the double layer integral alone, we use the identity (2.3.19) from [20],

8π
ilm ∂Ω where χ = 1, 1/2, 0 when x 0 is inside, on, and outside the boundary. We let l = 1 and define q i (x) = i1m x m = (0, −x 3 , x 2 ). Given this density function, we compute the left hand side of the identity (59) as the double layer potential using (4), (9b), (7), (10b). We then compare the computed values to the exact values, given by the right hand side of (59). It is worth noting that since we use subtraction in the double layer, for this example we get (q − q 0 ) · (x − x 0 ) ≡ 0, and therefore the integral is identically zero on the boundary. For a general ellipsoid, we use and test the identity at points inside and outside the surface at distance ≤ h away. Figure 5 shows the errors for the unit sphere and the ellipsoid a = 1, b = 0.6, c = 0.4. Again, we test three values of δ/h: 0.5, 1, and 2. The behavior is similar to the single layer integral near the surface. The errors reported are again absolute errors, with the largest magnitude of the solution about 1 for the sphere and about 0.6 for the ellipsoid. The comparison of the two cases is also affected by the differing number of points for given h. For reference, the number of quadrature points for the sphere and the ellipsoid for each h is given in Table 1.

Sum of single and double layer
One of the advantages of using boundary integral formulations is that jumps in the physical quantities across interfaces get incorporated into the integrals naturally. Specifically, the general integral formulation, expressed as the sum of the single and double layer integrals,  Here n is the outward unit normal, and the plus/minus signs denote the outside/inside of the boundary. To demonstrate the accuracy of this formulation, we use the following solution. On the inside, we assume the velocity is given by a point force singularity of strength b = (1, 0, 0), placed at y 0 = (2, 0, 0). The solution is given by the Stokeslet velocity and the stress tensor is whereŷ = y − y 0 , r = |ŷ|. We assume this data for the inside of the boundary, and take the solution to be u + = 0, σ + = 0 for the outside. The jumps [u] and [f ] are evaluated at the quadrature points using these inside/outside values. The single layer and double layer integrals in (61) are then computed as described in Sections 4.2 and 4.3. To compare to the exact solution on the boundary, we take the solution as the average of outside and inside, or half of the formula for u i in (62). This test allows us to check convergence when the formulation involves both the single and double layer potentials, using the high order regularization (50)-(52) for points on the surface, and the effect of corrections (10a)-(10b) for points off the surface. Figure 6 shows errors for the unit sphere, evaluated at the quadrature points only (left graph) and points inside and outside the surface that are distance ≤ h away. Figure 7 shows similar errors for the ellipsoid a = 1, b = 0.6, c = 0.4. As another test, we use the four-atom molecular surface as in [8], given by 4 k=1 exp(−|x−x k | 2 /r 2 ) = c, with centers ( √ 3/3, 0, − √ 6/12), (− √ 3/6, ±.5, − √ 6/12), (0, 0, √ 6/4) and r = .5, c = .6. The results for this surface are shown in Figure 8 and exhibit a similar behavior. The largest magnitude of the solution is about 0.079 for both the sphere and the ellipsoid, and about 0.073 for the molecular surface. See Table 1 for the number of quadrature points for these surfaces for different grid sizes. On the surface, we see high order convergence, O(h 5 ) for the sphere and O(h 4 ) for the ellipsoid, when regularization is chosen large enough, such as δ/h = 3. For smaller regularization parameter δ/h = 1 or less, observed convergence is O(h), so this is not recommended in practice. For points off the surface, the accuracy in the corrected solution is the predicted O(h 3 ) for regularization δ/h = 2.  Figure 6: Sum of single and double layer, errors for the unit sphere: over quadrature points (left) and over points distance ≤ h from the surface (right).

Flow due to an interface with different viscosities
Here we consider an example of an interface between two fluids with different viscosities, and an integral equation must be solved to find the interface velocity. The interface undergoes a discontinuity in the surface force [f ], while the velocity across the interface is continuous  Figure 8: Sum of single and double layer, errors for the molecular surface: over quadrature points (left) and over points distance ≤ h from the surface (right). [20]. The integral equation for the interface velocity is given by for x 0 ∈ ∂Ω, where µ 0 , µ 1 are the external and internal fluid viscosities and λ = µ 1 /µ 0 . The discontinuity in the surface force is given by [f ] = 2γHn − ∇ S γ, where γ is the surface tension, H is the mean curvature, and n is the outward unit normal [20,25]. In our numerical tests, we set µ 0 = 1, µ 1 = 2, and γ = 1 + x 2 1 . We solve the integral equation using successive evaluations, i.e., for N = 1, 2, ..., and u 0 = 0. We stop these iterations when the iteration error, defined as is below a prescribed tolerance, and |·| is the vector's Euclidean norm. We use the higher order regularization derived in Sec. 3.3. Since the exact solution is not known, we check the convergence rates by defining and taking either the max or the L 2 norm of this error over the surface points given by h, the larger of the two grid sizes used. These errors are shown in Table 2 for the unit sphere and the ellipsoid a = 1, b = 0.6, c = 0.4, with δ = 3h in both cases. In these tests, it took N = 8 iterations for the iteration error (66) to reach below 10 −10 for the sphere and around 12 iterations for the ellipsoid. The error in the solution (67) is larger especially for the ellipsoid, likely coming mostly from evaluating the single layer integral with the surface tension density. This can be remedied somewhat by computing the single layer integral with increased resolution before solving the integral equation, since the single layer is a nonhomogeneous term. For example, we solved the integral equation for each of the values of h, but in each case computed the Stokeslet integral at the needed points using the finer grid h = 1/256. As shown in the last section of Table 2, doing this reduces the error by over an order of magnitude.  Table 2: Flow due to an interface, for the unit sphere and the ellipsoid a = 1, b = 0.6, c = 0.4. Ellipsoid*: the single layer integral was computed using h = 1/256. Grid size h, max and L 2 norms of the error defined in (67). Regularization parameter δ = 3h.

Two interfaces close to each other
In the numerical simulation above, the integral equation is solved for a single interface, so the high order regularization of Sec. 3.3 is used and corrections are not present (unless one wishes to compute the flow off the surface). However, in the case of two or more surfaces that get close to each other, the error will deteriorate due to the near singularity issue and the high order regularization alone will not help. To demonstrate the importance of corrections, we modify the previous test in the following way. We use two unit spheres, one centered at the origin and the other at (2, 0, ), where = 1/16 3 . The integral equation is similar to (64) but with a sum of two single layer and two double layer integrals, one for each interface: for x 0 ∈ ∂Ω p , p = 1, 2, λ p = µ p /µ 0 . We take µ 0 = 1 and λ 1 = λ 2 = 2 as the viscosity ratio for each interface, and define the surface force as before using the surface tension γ = 1 + (x 1 − x c ) 2 , where x c is the x-coordinate of the center of the sphere. We again perform successive evaluations to compute the interface velocities, similar to (65). Regularization δ/h = 2 was used off the surface, while on the surface δ/h = 3 was used along with the high order regularization, as in the previous test. We solve the equation for h = 1/16, 1/32, 1/64, and estimate the convergence rates using the error definition of (67). Table 3 compares the errors for three solutions: direct, uncorrected, and corrected. In the direct solution, regularization (8a)-(8c) was used without corrections. The uncorrected solution was computed using the higher regularization (50)-(52) for same-surface integrals, i.e., when p = m in (68), but without corrections for the other-surface case, i.e., when p = m in (68). Finally, the corrected solution uses higher regularization for same-surface integrals, and corrections (10a)-(10b) for the other-surface integrals. The error reported in Table 3 is over the points on one sphere -it is essentially the same for the other sphere. Not only the error decreases in magnitude with corrections, but the convergence rate improves from first to third order. Figure 9 shows the error distribution (on a log scale) for the three cases. It is clear that the largest error is where the surfaces are near, and it does not improve overall when higher regularization is used for the same-surface integrals, unless the corrections are added for the other-surface integrals.

Conclusions
The numerical results we have performed in this paper are in agreement with the analytical prediction of uniformly third order spatial convergence for the computation of single and double layer integrals of the form (3a) and (3b). This is true when the evaluation point is near the surface, as is the case when two interfaces are close to each other. The accurate solution is obtained by regularizing the kernels and adding analytically derived correction terms to eliminate the first and second order regularization error terms. When the evaluation point is on the surface, a much higher accuracy is achieved without corrections, by improving the way the integrands are regularized. For this case, it might help to derive correction terms for the error due to the discretization of the integrals, although with an appropriate choice of the regularization parameter, this seems rarely necessary.
The error decays rapidly away from the surface, but our numerical results suggest that it might still be somewhat larger when the evaluation point is very close to the surface. One might experiment then with an interpolation technique such as [33], where the value very near the surface is interpolated from corrected values further away, rather than computed directly with corrections.
One of the advantages of the corrections method is that it does not increase the computational complexity of the overall method. Specifically, with N quadrature points and M evaluation points, computing the integrals will require O(N M ) CPU time, while the corrections add O(M ) to this. The computational efficiency of the algorithm can be improved then irrespective of the corrections, by using a fast summation algorithm such as a treecode [30] or a fast multipole method [27]. Such methods decrease the CPU time generally to O(M log N ) or O(M ), respectively.