In the Brownian dynamics simulation with hydrodynamic interactions, one needs to generate the total displacement vectors of Brownian particles consisting of two parts: a deterministic part which is proportional to the product of the Rotne-Prager-Yamakawa (RPY) tensor D [1, 2] and the given external forces F; and a hydrodynamically correlated random part whose covariance is proportional to the RPY tensor. To be more precise, one needs to calculate Du for a given vector u and compute √Dv for a normally distributed random vector v. For an arbitrary N-particle configuration, D is a 3N x 3N matrix and u, v are vectors of length 3N. Thus, classical algorithms require O(N2) operations for computing Du and O(N3) operations for computing √Dv, which are prohibitively expensive and render large scale simulations impossible since one needs to carry out these calculations many times in a Brownian dynamics simulation.
In this dissertation, we first present two fast multipole methods (FMM) for computing Du. The first FMM is a simple application of the kernel independent FMM (KIFMM) developed by Ying, Biros, and Zorin , which requires 9 scalar FMM calls. The second FMM, similar to the FMM for Stokeslet developed by Tornberg and Greengard , decomposes the RPY tensor into harmonic potentials and its derivatives, and thus requires only four harmonic FMM calls. Both FMMs reduce the computational cost of Du from O(N2) to O(N) for an arbitrary N-particle configuration.
We then discuss several methods of computing √Dv, which are all based on the Krylov subspace approximations, that is, replacing √Dv by p(D)v with p(D) a low degree polynomial in D. We first show rigorously that the popular Chebyshev spectral approximation method (see, for example, [5, 6]) requires √κ log 1/ε terms for a desired precision E, where K is the condition number of the RPY tensor D. In the Chebyshev spectral approximation method, one also needs to estimate the extreme eigenvalues of D. We have considered several methods: the classical Lanczos method, the Chebyshev-Davidson method, and the safeguarded Lanczos method proposed by Zhou and Li . Our numerical experiments indicate that K is usually very small when the particles are distributed uniformly with low density, and that the safeguarded Lanczos method is most effective for our cases with very little additional computational cost. Thus, when combined with the FMMs we described earlier, the Chebyshev approximation method with safeguarded Lanczos method as eigenvalue estimators essentially reduces the cost of computing √Dv from O(N3) to O(N) for most practical particle configurations. Finally, we propose to combine the so-called spectral Lanczos decomposition method (SLDM) (see, for example, ) and the FMMs to compute √Dv. Our numerical experiments show that the SLDM is generally more efficient than the popular Chebyshev spectral approximation method.
The fast algorithms developed in this dissertation will be useful for the study of diffusion limited reactions, polymer dynamics, protein folding, and particle coagulation as it enables large scale Brownian dynamics simulations. Moreover, the algorithms can be extended to speed up the computation involving the matrix square root for many other matrices, which has potential applications in areas such as statistical analysis with certain spatial correlations and model reduction in dynamic control theory.