{'fraction_non_alphanumeric': 0.027747867141555537, 'fraction_numerical': 0.00430713161600265, 'mean_word_length': 4.6925978312116925, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computer Algebra for Microhydrodynamics

I describe a method for computer algebra that helps with laborious calculations typically encountered in theoretical microhydrodynamics. The program mimics how humans calculate by matching patterns and making replacements according to the rules of algebra and calculus. This note gives an overview and walks through an example, while the accompanying code repository contains the implementation details, a tutorial, and more examples. The code repository is attached as supplementary material to this note, and maintained at https://github.com/jeinarsson/matte


I. INTRODUCTION
Microhydrodynamics concerns the motions of fluids and particles at low Reynolds numbers, where viscous forces dominate over inertial. In this regime the flow velocity and pressure fields u and p are governed by Stokes equations Although Eq. (1) is linear, solving it is often complicated by of one or more of • boundary conditions imposed on geometry in which Eq. (1) is not separable, • non-linear boundary conditions on fluid-fluid (or more complicated) interfaces, • non-linear forcing due to non-Newtonian, inertial, or other non-linear effects.
One option is to produce a numerical solution to the governing equations. When the nonlinearities dominate the solutions that may be the only option. Theoretical progress, however, may often be made by perturbation theory when there is a small parameter controlling the magnitude of the nonlinearity. Perturbation problems describe the first deviations from simple Stokes flow due to some nonlinearity. For example, the first effects of viscoelasticity is described by a 'second-order fluid' constitutive equation which renders the fluid stress nonlinear in the flow velocity gradients. Weak non-linear effects that directly affect the momentum balance in the fluid generally result in an equation of the form where f is any functional and is a small parameter. Another example is when considering fluid-fluid interfaces with surface tension, for example a drop suspended in fluid. In this case there are boundary conditions on the stress jump across the interface, applied at the unknown position of the interface: where σ andσ are the stress tensors on the two sides of the interface, and n is the unknown normal vector of the drop. This perturbation is controlled by a Capillary number Ca, which in the limit of Ca = 0 decouples the interface shape from the flow equations.
In perturbation theory we seek solutions to these problems on the form of a series, for example and similarly for all relevant variables 1 . By inserting the series into the governing equations and comparing order by order in we typically find a sequence of equations that begins with a simple Stokes flow problem, yielding u (0) , p (0) , and any other lowest order variables. Then follows a sequence of linear, but inhomogenous equations and/or boundary conditions for the higher orders, for example Even for low values of n, evaluating f (u (n−1) ) can be tedious. Solving for u (n) is even worse, and even if one avoids that by use of a reciprocal theorem the resulting integrals require substantial mathematical labor. Although this labor is straightforward in principle, the sheer amount of work may be discouraging, and prone to error if carried through. Nevertheless, in order to test our models -the perturbation equations and their associated assumptions about physical reality -we must produce a prediction.
In this note and the accompanying code I describe a set of software routines, implemented in Mathematica, that helps to manipulate algebraic expressions typically encountered in perturbation theory near Stokes flow. I describe the main idea in Section II, but leave the implementation details and usage to the tutorial notebooks. I walk through a simple example calculation in Section III.

II. ALGEBRA BY PATTERN MATCHING
We consider now solutions to Stokes equations in spherical geometry. As will be clear the method is not limited to spherical geometry, but it is certainly easier to deal with than, say, spheroidal geometry.
In terms of Cartesian tensors, the building block for solutions to Laplace's, and therefore Stokes', equations are For example, the fundamental solution to Stokes' equation, the Stokeslet, is Suppose we want to compute the derivative of the Stokeslet, ∂ i 3 G i 1 i 2 : A trained human immediately recognizes the pattern of a derivative of a sum, and knows that the linear operation applies term-wise: Focusing on the second term, a human recognizes the pattern of a derivative of a product, and applies the product rule: , and once more, A human performs these steps intuitively by training and experience, recognizing patterns and replacing them with the appropriate expression according to the rules of algebra and calculus. The next pattern may be that the derivative ∂ j r i = δ ij , and then perhaps δ ij r j = r i . The idea is that we humans perform calculations by matching patterns in the mathematical expressions to a catalogue of known, legal, rules of algebra and calculus. The goal is to implement a computer program that contains a number of patterns and their associated replacement rules, and let the computer iterate the pattern matching on a mathematical expression until it ceases to change.
In practice, however, we must carefully choose which patterns we implement. First of all, legal algebraic manipulations such as the examples above, are equalities and therefore valid in both directions. If we were to implement both the pattern δ ij r j → r i , and r i → δ ij r j the program will get stuck in an infinite loop switching back and forth between the two forms. A human practitioner has an opinion on what form is appropriate, and that opinion must be implemented in the choice of patterns. Second, a trained human typically knows a vast catalogue of different patterns and rules applicable in different situations, and many times disambiguate meaning based on context (for example, the implied summation in index notation may be over 2 or 3 or 7 dimensions, known only by context). It is far out of scope to implement any generally valid algebra software, but within a limited context with a reasonable set of mathematical objects, we can make substantial progress.
In the associated Mathematica notebooks I show that by careful choice of patterns and conventions we can implement a program that is useful for calculations in theoretical microhydrodynamics. We already discussed the role of patterns and rules for tensor contraction, differentiation, and so forth. Additionally the program must handle basic symmetry operations. For example, for a symmetric strain tensor E ij we must have E ij − E ji = 0. Moreover, we must recognize that E ij r j − E ik r k = 0, because dummy indices may be renamed. Finally, the code contains routines to collect terms and prettyprint them in a human-friendly format.
The resulting program is only around 350 lines of code, and the implementation is documented through a sequence of follow-along tutorial notebooks. I refer the interested reader to the tutorial in order to learn the details about syntax and usage.

III. EXAMPLE
Peery [1] perturbatively calculated the flow and pressure fields around a torquefree rigid spherical particle suspended in a linear flow of a second-order fluid. Here we extend that calculation to include the effects of particle translation, as induced by for example an external force [2]. Here we focus on the mathematical labor, and I will not dwell on details about the physics and interpretation. For convenience I include some screenshots of the program input and output here, but must refer the reader to the accompanying Mathematica notebook for the full program and output.
Consider a spherical particle moving with velocity v and rotating with angular velocity ω in a linear flow u ∞ = Er + Or. The second-order fluid constitutive equation for the stress is where The momentum balance ∇ · σ = 0 is thus The torque-free condition is The boundary conditions read u = v + ω × r on the particle surface, At O(1), with Wi = 0, this is a textbook Stokes flow problem. We solve it by an ansatz of the Stokeslet and its derivatives: This is an implementation of the ansatz in code: The functions partialr[g, i] and laplacer[g] represent the differential operators ∂g/∂r i and ∇ 2 g, respectively. The function getStokeslet is simply defined as The code is obviously not as dense as mathematical notation, but it is reasonably expressive. We form the boundary condition where we used that the particle rotates with the flow vorticity to lowest order 2 , that is ω (0) × r = Or. In code this tensorial boundary condition is: The program can collect all scalar prefactors of the unique tensorial groups. All these prefactors must vanish separately, which determines a 1 ..a 5 , like so: Now we may form u (0) by substituting the coefficients: Given u ij , the second-order fluid stress, and ultimately the right hand side of the next order equation. The momentum balance at this order is where u (1) = 0 , on the particle surface, where we used that the O(Wi) angular velocity is zero (this can be shown using a reciprocal theorem [2], details omitted here for brevity).
Computing the divergence of the second-order fluid stress is certainly a tedious task, as evident by the pages upon pages of algebra presented by Peery [1]. In the current program, we have That's it. With the inhomogenous Stokes equation (18) fully specified, it remains to find a solution. We could proceed like Peery [1]: take the divergence of Eq.
is given by where F denotes Fourier transform and k is the reciprocal variable. The patterns and rules for the Fourier transform of polyadic expressions like r i r j ../r m are in fact rather involved, because it turns out that the angular eigenfunctions of the 3D Fourier transform are the spherical harmonics of degree l, but the polyadic expressions are typically linear combinations of different values of l. Therefore we must separate the parts of the polyadics of different degree l, transform, and then recombine. See Einarsson and Mehlig [2] for a detailed derivation. However, once that is implemented, Eq. (22) is a convenient route to a particular solution, namely: (note that output is cut off in screenshot) Finally, we add a homogenous solution to u (1)p such that the sum satisfies the boundary conditions (20). In this case it is again convenient to make an ansatz of G ij and its derivatives: See the notebook for the implementation of this ansatz. The program readily simplifies and evaluates this expression at r = 1: The runtime of this entire calculation is about ten seconds on a 2016 laptop.

IV. DISCUSSION
In this note I have presented a method of computer algebra that helps with some of the most laborious tasks in theoretical microhydrodynamics. Perhaps this demonstration can inspire attempts at new results by computer algebra, and push the limit of how complicated models we tackle analytically.
The examples shown in this note all pertain to microhydrodynamics in spherical geometry, but it should be clear that the general method applies as long as one can invent patterns and rules for the mathematical objects in question. For example, I employed an early version of this method for my thesis work on the rotation rate of spheroidal particles [3]. In that case the solutions to Stokes' equation are polyadics in the position vector r, and the spheroid axis of symmetry n, multiplied with the scalar functions where c is a constant related to the spheroidal geometry. While these scalar functions are much more complicated than the spherical counterpart 1/r m , they still obey certain rules, for example Eq. (B3) in [3]: The point is, if you can perform the algebra by pattern matching and replacement, this method can help with the mechanical labor.

ACKNOWLEDGMENTS
I would like to thank Kristian Gustafsson who once upon a time got me started on the pattern matching path, and Joe Barakat for inspiring me to finally clean up