Against Modularity

This paper introduces the elliptic package of R routines, for numerical calculation of elliptic and related functions. Elliptic functions furnish interesting and instructive examples of many ideas of complex analysis, and the package illustrates these numerically and visually. A statistical application in fluid mechanics is presented.


Introduction
The elliptic functions crop up here and there in diverse areas of applied mathematics such as cosmology (Kraniotis and Whitehouse 2002), chemical engineering (Koopman and Lee 1991), dynamical systems (Kotus and Urbański 2003), and quantum mechanics (Chow 2002); here they are applied to fluid mechanics (Johnson andMcDonald 2004, 2005).They also provide interesting and instructive non-elementary examples of many results in complex analysis such as Cauchy's integral theorem and the residue theorem.
In this paper I introduce elliptic, a new R package for numerical calculation of Weierstrass and Jacobi elliptic functions, theta functions and modular functions.The emphasis is on efficient numerical calculation, and informative visualization techniques.

Elliptic functions
This section gives a very brief introduction to elliptic functions.For more detail and rigorous derivations, the reader is referred to the classic literature: the standard reference would be Whittaker and Watson (1952).Chandrasekharan (1985) approaches the field from a more modern perspective, and Abramowitz and Stegun (1965) provide the definitive reference work for the case of real invariants.
A meromorphic function f is said to be elliptic if ∃ ω 1 , ω 2 ∈ C with ω 2 /ω 1 ∈ C\R such that (1) whenever f (z) is defined and m, n ∈ Z. Notation in this paper is consistent with that of Abramowitz and Stegun (1965); ω 1 and ω 2 are called the half periods.In 1862, Weierstrass introduced his ℘ function which is defined as (2) The ℘ function is, in a well defined sense, the simplest nontrivial elliptic function (Whittaker and Watson 1952).Given this, we have a Taylor expansion of the form with where a prime indicates summation over Z 2 excluding (m, n) = (0, 0).For reasons to be made clear in section 4.1.1,g 2 and g 3 are known as the invariants.Other equivalent forms for ℘ include its differential equation and the relation which is equivalent to w = ℘(z).
Related functions include the zeta function ζ(z), defined by and the sigma function σ(z), defined by is analytic except for points on the lattice of periods, at which it has simple poles with residue 1.One classic result is due to Legendre: if ω 1 , ω 2 is a pair of basic periods 1 , with Im(ω 2 /ω 1 ) > 0, then where 1 A pair of basic periods is one that generates the period lattice.Basic periods are not unique as many pairs of periods may generate the same lattice.However, there is one pair of basic periods, the fundamental periods that are, in a well-defined sense, optimal (Chandrasekharan 1985).

Jacobian elliptic functions
Jacobi approached the description of elliptic functions from a different perspective (Weisstein 2005).Given m = k 2 and m 1 with m + m 1 = 1, Jacobi showed that if are elliptic with periods and The Jacobian elliptic functions are encountered in a variety of contexts and bear a simple analytical relation with the Weierstrass ℘ function.

Numerical evaluation and Jacobi's theta functions
Although equation 2 is absolutely convergent, it converges too slowly to be of use in practical work, and an alternative definition is needed.
Jacobi presented his four theta functions in 1829 and, although they have interesting properties in their own right, here they are used to provide efficient numerical methods for calculation of the elliptic functions above.They are defined as follows: q n(n+1) cos(2n + 1)z ( 14) It may be seen that, provided |q| < 1, the series converges for all z ∈ C. Indeed, the convergence is very rapid: the number of correct significant figures goes as the square of the number of terms.It should be noted that there are different notations in use, both for the four function names, and for the independent variables.
All the functions described in section 2 may be expressed in terms of the theta functions.This is the default method in elliptic, although alternative algorithms are implemented where possible to provide a numerical and notational check.
For example, Weierstrass's ℘ function is given by where q = e iω2 /ω 1 ; the other functions have similar theta function definitions.

Package elliptic in use
This section shows elliptic being used in a variety of contexts.First, a number of numerical verifications of the code are presented; then, elliptic and related functions are visualized using the function view(); and finally, the package is used to calculate flows occurring in an oceanographic context.

Numerical verification
Work in the field of elliptic functions is very liable to mistakes 2 , and package elliptic includes a number of numerical checks to guard against notational inexactitude.These checks generally use the convenient (trivial) function maxdiff() that shows the maximum absolute difference between its two arguments: For example, we may compare the output of P(), which uses equation 17, against the straightforward Laurent expansion, used by P.laurent(): showing reasonable agreement; note that function P() uses the conceptually distinct theta function formula of equation 17.Package elliptic includes a large number of such numerical verification tests in the test suite provided in the package, but perhaps more germane is the inclusion of named identities appearing in Abramowitz and Stegun (1965).For example, consider function e18.10.9(),named for the equation number of the identity appearing on page 650.This function returns the difference between the (algebraically identical) left and right hand side of three grouped identities: where q = e −πK /K .From the manpage: > abs(e18.10.9(parameters(g= g))) [1] 3.614624e-15 4.440892e-15 3.552714e-15 again showing reasonably accurate numerical results, but perhaps more importantly explicitly verifying that the notational scheme used is internally consistent.
Although the examples above use the invariants g2 and g3 to define the elliptic function and its periods, sometimes the fundamental periods are known and the invariants are desired.This is done by function g.fun(), which takes the fundamental periods as arguments and returns the two invariants g 2 and g 3 .Observe that there are many pairs of basic periods that generate the same lattice-see figure 1-but it usual to specify the unique fundamental periods as this pair usually has desirable numerical convergence properties.

Unimodularity
Many functions of the package are unimodular.The invariants g 2 and g 3 are defined in terms of a pair of basic periods ω 1 and ω 2 .However, any pair of basic periods should have the same invariants, because any pair of basic periods will define the same elliptic function (hence the name).Basic period pairs are related by a unimodular transformation: if ω 1 , ω 2 and ω1 , ω2 are two pairs of basic periods then there exist integers a, b, c, d with ad − bc = 1 and where M is unimodular: that is, an integer matrix with a determinant of unity.In this context, unimodular matrices (and the transformations they define) are interesting because any two pairs of basic periods are related by a unimodular transformation.
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Re(z) The package includes functions that generate unimodular matrices.The underlying function is congruence(), which generates 2 × 2 integer matricies with a determinant of 1, given the first row.For example: (observe that the determinant of M is unity) and thus we may verify the unimodularity of, for example, g.fun() by evaluating the invariants for a pair of fundamental periods, and comparing this with the invariants calculated for a pair of basic periods that are related to the fundamental periods by a unimodular transformation (here M).In R idiom: > maxdiff(g.fun(o),g.fun(M %*% o, maxiter = 800)) [1] 1.634621e-13 showing that the invariants for period pair o = (1, i) T are almost identical to those for period pair o = Mo = (4 + 9i, 3 + 7i) T .Observe that the latter evaluation requires many more iterations for accurate numerical evaluation: this behaviour is typically encountered when considering periods whose ratio is close to the real axis.
In addition, function unimodular() generates unimodular matrices systematically, and function unimodularity() checks for a function's being unimodular.

Contour integration and the residue theorem
As noted in section 2, the zeta function ζ(z) possesses a simple pole of residue 1 at the origin.The residue theorem would imply that when the contour is taken round a path that encircles the origin but no other poles.This may be verified numerically using elliptic's myintegrate suite of functions, which generalize the stats package's integrate() function to the complex plane.Here, function integrate.contour() is used to integrate round the unit circle.This function takes three arguments: first, the function to be integrated; second, a function that describes the contour to be integrated along; and third, a function that describes the derivative of the contour.We may now integrate over a closed loop, using arguments u and udash which specify a contour following the unit circle: > jj <-integrate.contour(Zeta,u, udash) > maxdiff(jj, 2 * pi * (0+1i)) [1] 9.344997e-16 showing reasonable numerical accuracy.Compare Weierstrass's ℘ function, which has a second order pole at the origin: > abs(integrate.contour(WeierstrassP,u, udash)) [1] 4.002966e-16 The PARI system Perhaps the most convincing evidence for numerical accuracy and consistency of notation in the software presented here is provided by comparison of the package's results with that of PARI (Batut et al. 2000).The PARI system is an open-source project aimed at number theorists, with an emphasis on pure mathematics; it includes some elliptic function capability.Function P.pari() of package elliptic calls the pari system directly to evaluate elliptic functions from within an R session, enabling quick verification: again showing reasonable agreement, this time between two independent computational systems.

Visualization of complex functions
In the following, a Weierstrass elliptic function with invariants of 1 + i and 2 − 3i will be considered.The half periods ω 1 , ω 2 are first evaluated: and these may be visualized by using latplot(), as in figure 1. Figure 2 shows the real part of such a function, shown over part of the complex plane, and figure 3 shows the same function using the view() function.
The σ function with the same invariants is visualized in figure 4, showing that its zeros lie on the same lattice as figure 1.
Figure 5 shows the zeta function, and figure 6 shows Jacobi's "sn" function.

Potential flow
One application of complex analysis is to fluid dynamics.In particular, potential flow (steady, two-dimensional, inviscid, incompressible) may be studied with the aid of analytic complex functions.Here I show how the elliptic functions discussed in this paper may be used to simulate potential flow in a rectangular domain.Although the tenets of potential flow appear to be absurdly idealized3 , it is nevertheless a useful technique in many branches of practical fluid mechanics: it is often used to calculate a "theoretical" flowfield with which measured velocities may be compared.A short sketch of potential theory is given here but the reader is referred to Milne-Thomson (1949) for a full exposition.Briefly, we define a complex potential w(z) to be a complex function and observe that both φ and ψ obey Laplace's equation if w is differentiable.Given this, we may take the velocity vector v = (v x , v y ) of the fluid to be and observe that streamlines are given by contours of equal ψ; contours of equal φ are called equipotential lines.The two systems of lines cross at right angles (this follows from the Cauchy-Riemann conditions).
Consider, for example, the function w(z) = z 2 , whose associated flow field is shown in figure 7.This corresponds to a stagnation point, at which the speed vanishes; the streamlines (solid) intersect the equipotential lines (dashed) at right angles.
Now consider a slightly more complicated case.A point source of strength m at z 0 may be represented by the function m log(z − z 0 ) (a sink corresponds to m < 0).Any finite number of sources or sinks may be combined, as in i m i log(z − z i ) where the i th source is at z i and has strength m i , because the system is linear4 .Figure 8 shows two sources and two sinks, all of equal strength.Because the flowfield is symmetric with respect to the real axis, there is no flux across it; we may therefore ignore the flow in the lower half plane (ie {z : Im(z) < 0}) and consider the flow to be bounded below by the real axis.This is known as the method of images (Milne-Thomson 1949).Now, one may transform a potential flowfield into another form using a conformal mapping from the z-plane to the ζ-plane, traditionally denoted This technique finds application when flow is desired (in the ζ-plane) that obeys some specific boundary condition that is simple to specify in the z-plane.
In this case, we seek a conformal transformation that maps the upper half plane to a rectangle.If we consider the flowfield shown in figure 8, then the map given by takes the upper half plane of the ζ-plane to a rectangle in the z-plane (Milne-Thomson 1949).Using equation 10, this is equivalent to z = sn(ζ; m), where sn(•, •) is a Jacobian elliptic function and m a constant of integration.
Figure 9 shows the resulting flow field: observe how the flow speed, which is proportional to the spacing between the streamlines, is very small near the left-hand edge of the rectangle.

Bayesian analysis of potential flow
When considering potential flows, it is often necessary to infer the locations of singularities in the flow from sparse and imperfect data (Johnson and McDonald 2004).
Here, I apply the methods of Kennedy and O'Hagan (2001a) and Kennedy and O'Hagan (2001b) (hereafter KOH and KOHa respectively) using the BACCO package (Hankin 2005) to assess some characteristics of potential flow in a rectangle.
Kennedy and O'Hagan considered the following inference problem for a set of parameters θ ∈ R q that are inputs to a computer program.Given an independent variable x ∈ R n , and a set of scalar ("field") observations z = z(x), they assume where ρ is a constant of proportionality (notionally unity); η(•, •) a Gaussian process with unknown coefficients; θ the true, but unknown parameter values; δ(•) a model inadequacy term (also a Gaussian process with unknown coefficients); and ∼ N (0, λ 2 ) uncorrelated normal observational errors.
Inferences about η(•, •) are made from point observations of the process: Kennedy and O'Hagan call these the "code observations" on the grounds that their chief motivation was the understanding of complex computer codes.
Here, potential flow in a rectangle is considered.The source is at one corner of the rectangle, which is considered to have lower left point (−1, 0) and upper right point (1, 1).The position of the sink is unknown.
I now show how the position of the sink may be inferred from a sparse and noisy set of observed fluid speeds.Similar inference problems are encountered in practice when considering oceanographic flows such as those occurring near deep sea vents, although the geometry is generally more complex than considered here.
The independent variable x is the two-dimensional position within the rectangle, and the field observation z(x) is the fluid speed at that point, plus obervational error .The parameter set θ thus has two degrees of freedom corresponding to the x− and y− coordinates of the sink.
Field observations will be obtained numerically, using the elliptic package.The simulated flowfield has a sink at a known position-in this case the geometric centre of the rectangleand Bayesian methods will be used to infer its position using only fluid speed data.
In the terminology of KOHa, dataset y corresponds to modelled fluid speed, obtained from the appropriate simulations carried out with the sink placed at different locations within the rectangle.Dataset z corresponds to field observations, which in this case is fluid speed at several points in the rectangle, obtained from simulations with the sink at the centre of the rectangle.
The code evaluation design matrix D1 is chosen according to a random Latin hypercube design, and the observation is calculated using the elliptic package: > head(D1) So the first line shows a simulation with the sink at (-0.17,0.15); the log of the fluid speed at (0.77, 0.35) is -1.87.There are a total of 30 such observations.Figure 10 shows these points superimposed on the "true" flow field.
The field observations are similarly determined: > head(D2) showing that the first field observation, at (0.26, 0.63), is -2.5.There are a total of 31 such observations.Figure 11 shows the first code observation in the context of the "true" flow field.
Kennedy and O'Hagan give, inter alia, an expression for the likelihood of any value of θ being the true parameter set (in this case, the true position of the sink) in terms of the code evaluations and field observations.
Here, function support() calculates the log-likelihood for a pair of coordinates of the sink.This may be evaluated at the centre of the rectangle, and again at the top left corner: > support(c(0, 1/2)) [1] 0 > support(c(-1, 1)) [1] -52.29362 showing, as expected, that the support is very much larger at the centre of the rectangle than the edge (here the arbitrary additive constant is such that the support at c(0,1/2) is exactly zero).It is now possible to identify the position of the sink that corresponds to maximum support using numerical optimization techniques: (mle <-optim(c(0,1/2),support)) [1] -0.1842530 0.6190824 Thus the maximum likelihood estimate for the sink is a distance of about 0.2 from the true position.The support at this point is about 3.9 units of likelihood: > support(mle) [1] 3.908104 Discussion of Bayesian statistical analysis The above example shows the ideas of KOH being applied straightforwardly, but with the novel twist of θ being interpreted as physical characteristics of a fluid flow.In this case θ is the coordinates of the sink.
The MLE is better supported than the true position by about 3.9 units of likelihood: thus, in the language of Edwards (1992), the hypothesis of θ true = (0, 0.5) would not be rejected if one accepted Edwards's 2 units of likelihood per degree of freedom.
The discrepancy between θ and θ true (a distance of about 0.2) may be due to due to the coarseness of the form adopted for the basis functions, and better results might be obtained by using a more sophisticated system of model inadequacy than the simple linear form presented here.
The methods of KOH allow one to make statistically robust statements about the physical characteristics of an interesting flow that are difficult to make in any other way.

Conclusions
Elliptic functions are an interesting and instructive branch of complex analysis, and are frequently encountered in applied mathematics: here they were used to calculate a potential flow field in a rectangle.
This paper introduced the R package elliptic, which was then used in conjunction with Bayesian statistical methods (the BACCO bundle) to make statistically sound inferences about a flow with uncertain parameters: in this case the position of the sink was estimated from a sparse and noisy dataset.

Figure 1 :
Figure 1: The lattice generated by ℘(z; 1 + i, 2 − 3i); fundamental period parallelogram shown in light gray and a basic period parallelogram shown in darker gray

Figure 7 :Figure 9 :
Figure 7: Potential flow on the complex plane: field corresponding to the function (z) = z 2 .Solid lines represent streamlines and dotted lines represent equipotentials; these intersect at right angles.Note stagnation point at the origin

Figure 10 :
Figure 10: Streamlines of first code observation point; field observation point shown as a cross.The sink is at (-0.17,0.15)