\documentclass[nojss]{jss} \usepackage{dsfont} \usepackage{bbm} \usepackage{amsfonts} \usepackage{wasysym} \usepackage{wrapfig2} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% declarations for jss.cls %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% just as usual \author{Robin K. S. Hankin} \title{Introducing \pkg{elliptic}, an \proglang{R} package for elliptic and modular functions} %\VignetteIndexEntry{A vignette for the elliptic package} %% for pretty printing and a nice hypersummary also set: %% \Plainauthor{Achim Zeileis, Second Author} %% comma-separated \Plaintitle{Introducing elliptic, an R package for elliptic and modular functions} \Shorttitle{Elliptic functions with \proglang{R}} %% an abstract and keywords \Abstract{ To cite the package in publications, use \citet{hankin2006}. This paper introduces the \pkg{elliptic} package of \proglang{R} routines, for numerical calculation of elliptic and related functions. Elliptic functions furnish interesting and instructive examples of many ideas of complex analysis, and package \pkg{elliptic} illustrates these numerically and visually. A statistical application in fluid mechanics is presented. } \Keywords{Elliptic functions, modular functions, Weierstrass elliptic functions, visualization of complex functions} \Address{ Robin K. S. Hankin\\ The University of Stirling\\ E-mail: \email{hankin.robin@gmail.com} } %% need no \usepackage{Sweave.sty} \SweaveOpts{echo=FALSE} \begin{document} <>= require(elliptic,quietly=TRUE) @ <>= n <- 400 n_BACCO <- 40 @ \section{Introduction} \begin{wrapfigure}{r}{30mm} \includegraphics[width=1in]{\Sexpr{system.file("help/figures/elliptic.png",package="elliptic")}} \end{wrapfigure} % The elliptic functions crop up here and there in diverse areas of applied mathematics such as cosmology~\citep{kraniotis2002}, chemical engineering~\citep{koopman1991}, dynamical systems~\citep{kotus2003}, and quantum mechanics~\citep{chow2002}; here they are applied to fluid mechanics~\citep{johnson2004,johnson2005}. 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 \pkg{elliptic}, a new \proglang{R} package for numerical calculation of Weierstrass and Jacobi elliptic functions, theta functions and modular functions; it is based on~\citet{hankin2006}. The emphasis is on efficient numerical calculation, and informative visualization techniques. The package is available on CRAN, \url{https://CRAN.R-project.org/package=elliptic} \citep{rcore2005}. \section{Elliptic functions}\label{section:introduction} 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~\cite{whittaker1952}. \cite{chandrasekharan1985} approaches the field from a more modern perspective, and \cite{abramowitz1965} provide the definitive reference work for the case of real invariants. A meromorphic function~$f$ is said to be elliptic if~$\exists\,\omega_1,\omega_2\in\mathbbm{C}$ with~$\omega_2/\omega_1\in\mathbbm{C}\backslash\mathbbm{R}$ such that \begin{equation} f(z)=f(z+2m\omega_1+2n\omega_2) \end{equation} whenever~$f(z)$ is defined and~$m,n\in\mathbbm{Z}$. Notation in this paper is consistent with that of~\citet{abramowitz1965}; $\omega_1$ and~$\omega_2$ are called the {\em half periods}. In 1862, Weierstrass introduced his~$\wp$ function which is defined as \begin{equation}\label{direct.sum} \wp(z)= \frac{1}{z^2}+ \sum_{m,n\in\mathbbm{Z}\atop m,n\neq 0} \left\{ \frac{1}{\left(z-2m\omega_1-2n\omega_2\right)^2} -\frac{1}{\left( 2m\omega_1+2n\omega_2\right)^2} \right\}. \end{equation} The~$\wp$ function is, in a well defined sense, the simplest nontrivial elliptic function~\citep{whittaker1952}. Given this, we have a Laurent expansion of the form \begin{equation} \wp(z)-z^{-2}=\frac{1}{20}g_2z^2+\frac{1}{28}g_3z^4+O(z^6) \end{equation} with \begin{equation} g_2=60{\sum}'\frac{1}{\left(2m\omega_1+2n\omega_2\right)^4}, \qquad g_3=140{\sum}'\frac{1}{\left(2m\omega_1+2n\omega_2\right)^6}, \end{equation} where a prime indicates summation over~$\mathbbm{Z}^2$ excluding~$(m,n)=(0,0)$. For reasons to be made clear in section~\ref{section.unimodularity}, $g_2$ and~$g_3$ are known as the {\em invariants}. Other equivalent forms for~$\wp$ include its differential equation \begin{equation}\label{P.differential.eqn.definition} \left(\frac{d\wp}{dz}\right)^2=4\wp^3-g_2\wp-g_3 \end{equation} and the relation \begin{equation}\label{P.integral.definition} z=\int_{t=w}^\infty\frac{dt}{\sqrt{4t^3-g_2t-g_3}} \end{equation} which is equivalent to~$w=\wp(z)$. Related functions include the zeta function~$\zeta(z)$, defined by \begin{equation}\label{zeta.definition} \frac{d\zeta(z)}{dz}=-\wp(z) \end{equation} and the sigma function~$\sigma(z)$, defined by \begin{equation}\label{sigma.definition} \frac{d\log\sigma(z)}{dz}=\zeta(z),\qquad{\lim_{\mbox{\tiny $z\longrightarrow 0$}}}\left[\frac{\sigma(z)}{z}\right]=1 \end{equation} (neither~$\sigma(z)$ nor~$\zeta(z)$ is elliptic). It may be shown\label{zeta.analytic} that~$\zeta(z)$ 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~$\omega_1,\omega_2$ is a pair of basic periods\footnote{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 {\em fundamental} periods that are, in a well-defined sense, optimal~\citep{chandrasekharan1985}.}, with~$\rm{Im}(\omega_2/\omega_1)>0$, then \begin{equation} \eta_1\omega_2-\eta_2\omega_1=\pi i\label{legendre} \end{equation} where~$\eta_1=\zeta(\omega_1)$ and~$\eta_2=\zeta(\omega_2)$. \subsection{Jacobian elliptic functions} Jacobi approached the description of elliptic functions from a different perspective~\citep{weisstein2005}. Given~$m=k^2$ and~$m_1$ with~$m+m_1=1$, Jacobi showed that if \[ u=\int_{t=0}^\phi\frac{dt}{\sqrt{(1-t^2)(1-mt^2)}} \] the functions~${\rm sn}(u,k)$, ${\rm cn}(u,k)$ and~${\rm dn}(u,k)$ defined by \begin{equation}\label{sn.definition} {\rm sn} u=\sin\phi,\qquad {\rm cn} u=\cos\phi,\qquad {\rm dn} u=\sqrt{1-k^2\sin^2\phi} \end{equation} are elliptic with periods \begin{equation} K=\int_{\theta=0}^{\pi/2}\frac{d\theta}{\sqrt{1-m\sin^2\theta}} \end{equation} and \begin{equation} iK'=i\int_{\theta=0}^{\pi/2}\frac{d\theta}{\sqrt{1-m_1\sin^2\theta}}. \end{equation} The Jacobian elliptic functions are encountered in a variety of contexts and bear a simple analytical relation with the Weierstrass~$\wp$ function. \section{Numerical evaluation and Jacobi's theta functions} Although equation~\ref{direct.sum} 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: \begin{eqnarray}\label{theta.defn} \theta_1(z,q)&=&2q^{1/4}\sum_{n=0}^\infty(-1)^nq^{n(n+1)}\sin(2n+1)z\\ \theta_2(z,q)&=&2q^{1/4}\sum_{n=0}^\infty q^{n(n+1)}\cos(2n+1)z\\ \theta_3(z,q)&=&1+2\sum_{n=1}^\infty q^{n^2}\cos 2nz\\ \theta_4(z,q)&=&1+2\sum_{n=1}^\infty(-1)^n q^{n^2}\cos 2nz \end{eqnarray} It may be seen that, provided~$|q|<1$, the series converges for all~$z\in\mathbbm{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~\ref{section:introduction} may be expressed in terms of the theta functions. This is the default method in \pkg{elliptic}, although alternative algorithms are implemented where possible to provide a numerical and notational check. For example, Weierstrass's~$\wp$ function is given by \begin{equation}\label{P.in.terms.of.theta} \wp(z)= \frac{\pi^2}{4\omega_1^2}\left( \frac{\theta_1'(0,q)\theta_2(v,q)}{\theta_2(0,q)\theta_1(v,q)} \right)^2 \end{equation} where~$q=e^{i\omega_2/\omega_1}$; the other functions have similar theta function definitions. <>= <>= require(elliptic) require(emulator) require(calibrator) @ \section[Package ''elliptic'' in use]{Package \pkg{elliptic} in use} This section shows \pkg{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 \code{view()}; and finally, the package is used to calculate flows occurring in an oceanographic context. The primary function in package \pkg{elliptic} is~\code{P()}: this calculates the Weierstrass~$\wp$ function, and may take named arguments that specify either the invariants~\code{g} or half periods~\code{Omega}: <>= z <- 1.9+1.8i P(z,g=c(1,1i)) P(z,Omega=c(1,1i)) @ \subsection{Numerical verification} Work in the field of elliptic functions is very liable to mistakes\footnote{\cite{abramowitz1965} state that there is a ``bewildering'' variety of notations in use; the situation has become more confusing in the intervening 40 years.}, and package \pkg{elliptic} includes a number of numerical checks to guard against notational inexactitude. These checks generally use the convenient (trivial) function \code{maxdiff()} that shows the maximum absolute difference between its two arguments: <>= maxdiff <- function(x,y){max(abs(x-y))} @ For example, we may compare the output of \code{P()}, which uses equation~\ref{P.in.terms.of.theta}, against the straightforward Laurent expansion, used by \code{P.laurent()}: <>= g <- c(3,2+4i) z <- seq(from=1,to=0.4+1i,len=34) <>= maxdiff(P(z,g), P.laurent(z,g)) @ showing reasonable agreement; note that function \code{P()} uses the conceptually distinct theta function formula of equation~\ref{P.in.terms.of.theta}. Package \pkg{elliptic} includes a large number of such numerical verification tests in the \code{test} suite provided in the package, but perhaps more germane is the inclusion of named identities appearing in \cite{abramowitz1965}. For example, consider function~\code{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: \begin{eqnarray} 12\omega_1^2e_1 &=& \hphantom{-}\pi^2\left[\theta_3^4(0,q)+\theta_4^4(0,q)\right]\nonumber\\ 12\omega_1^2e_2 &=& \hphantom{-}\pi^2\left[\theta_2^4(0,q)-\theta_4^4(0,q)\right]\\ 12\omega_1^2e_3 &=& - \pi^2\left[\theta_3^4(0,q)+\theta_4^4(0,q)\right]\nonumber \end{eqnarray} where~$q=e^{-\pi K'/K}$. From the manpage: <>= abs( e18.10.9(parameters(g=g))) @ 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~\code{g2} and \code{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 \code{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~\ref{latplot}---but it usual to specify the unique {\em fundamental periods} as this pair usually has desirable numerical convergence properties. \begin{figure}[htbp] \begin{center} <>= jj <- parameters(g=c(1+1i,2-3i))$Omega latplot(jj,xlim=c(-4,4),ylim=c(-4,4),xlab="Re(z)", ylab="Im(z)") polygon(Re(c(jj[1],sum(jj),jj[2],0)), Im(c(jj[1],sum(jj),jj[2],0)),lwd=2,col="gray90",pch=16,cex=3) kk <- -c(3*jj[1] + 2*jj[2] , jj[1] + jj[2]) #det(matrix(c(3,2,1,1),2,2,T)==1 polygon(Re(c(kk[1],sum(kk),kk[2],0)), Im(c(kk[1],sum(kk),kk[2],0)),lwd=2,col="gray30",pch=16,cex=3) @ \caption{The\label{latplot} lattice generated by~$\wp(z;1+i,2-3i)$; fundamental period parallelogram shown in light gray and a basic period parallelogram shown in darker gray} \end{center} \end{figure} \subsubsection{Unimodularity}\label{section.unimodularity} Many functions of the package are {\em unimodular}. The invariants~$g_2$ and~$g_3$ are defined in terms of a pair of basic periods~$\omega_1$ and~$\omega_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~$\omega_1,\omega_2$ and~$\tilde{\omega}_1,\tilde{\omega}_2$ are two pairs of basic periods then there exist integers~$a,b,c,d$ with~$ad-bc=1$ and \[ \left( \begin{array}{cc} a&b\\ c&d \end{array} \right) \left( \!\! \begin{array}{c} \omega_1\\ \omega_2 \end{array} \!\! \right)= \left(\!\! \begin{array}{c} \tilde{\omega}_1\\ \tilde{\omega}_2 \end{array} \!\! \right) \] Formally, a unimodular function~$f(\cdot,\cdot)$ is one with arity~2---it is conventional to write~$f(\mathbf{v})=f(a,b)$---and for which \begin{equation} f(\mathbf{v})=f(\mathbf{M}\mathbf{v})\end{equation} where~$\mathbf{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. The package includes functions that generate unimodular matrices. The underlying function is \code{congruence()}, which generates~$2\times 2$ integer matrices with a determinant of~1, given the first row. For example: <>= M <- congruence(c(4,9)) @ (observe that the determinant of~$\mathbf{M}$ is unity) and thus we may verify the unimodularity of, for example, \code{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~$\mathbf{M}$). In \proglang{R} idiom: <>= o <- c(1,1i) <>= maxdiff(g.fun(o), g.fun(M %*% o,maxiter=800)) @ showing that the invariants for period pair~$o=(1,i)^T$ are almost identical to those for period pair~$o'=\mathbf{M}o=(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 \code{unimodular()} generates unimodular matrices systematically, and function \code{unimodularity()} checks for a function's being unimodular. \subsubsection{Contour integration and the residue theorem} As noted in section~\ref{zeta.analytic}, the zeta function~$\zeta(z)$ possesses a simple pole of residue~1 at the origin. The residue theorem would imply that \[ \varoint\zeta(z)\,dz=2\pi i \] when the contour is taken round a path that encircles the origin but no other poles. This may be verified numerically using \pkg{elliptic}'s \code{myintegrate} suite of functions, which generalize the \pkg{stats} package's \code{integrate()} function to the complex plane. Here, function \code{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~\code{u} and~\code{udash} which specify a contour following the unit circle: <>= u <- function(x){exp(pi*2i*x)} udash <- function(x){pi*2i*exp(pi*2i*x)} Zeta <- function(z){zeta(z,g)} Sigma <- function(z){sigma(z,g)} WeierstrassP <- function(z){P(z,g)} @ <>= jj <- integrate.contour(Zeta,u,udash) <>= maxdiff(jj, 2*pi*1i) @ showing reasonable numerical accuracy. Compare Weierstrass's~$\wp$ function, which has a second order pole at the origin: <>= abs(integrate.contour(WeierstrassP,u,udash)) @ \subsubsection[The PARI system]{The \proglang{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 \proglang{PARI}~\citep{batut2000}. The \proglang{PARI} system is an open-source project aimed at number theorists, with an emphasis on pure mathematics; it includes some elliptic function capability. Function \code{P.pari()} of package \pkg{elliptic} calls the \code{pari} system directly to evaluate elliptic functions from within an \proglang{R} session, enabling quick verification: \begin{Schunk} \begin{Sinput} > omega <- c(1,1i) \end{Sinput} \end{Schunk} \begin{Schunk} \begin{Sinput} > z <- seq(from=pi,to=pi*1i,len=10) \end{Sinput} \end{Schunk} \begin{Schunk} \begin{Sinput} > maxdiff(P.pari(z,Omega=omega), P(z,params=parameters(Omega=omega))) \end{Sinput} \begin{Soutput} [1] 2.760239e-14 \end{Soutput} \end{Schunk} again showing reasonable agreement, this time between two independent computational systems. \subsection{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 $\omega_1,\omega_2$ are first evaluated: <>= jj.omega <- half.periods(g=c(1+1i,2-3i)) @ and these may be visualized by using \code{latplot()}, as in figure~\ref{latplot}. Figure~\ref{P.persp.re} shows the real part of such a function, shown over part of the complex plane, and figure~\ref{P.view} shows the same function using the \code{view()} function. <>= x <- seq(from=-4, to=4, len=n) y <- x z <- outer(x,1i*x, "+") f <- P(z, c(1+1i,2-3i)) @ %% Thanks to Dario Strbenac for the following structure <>= png("wp_figure.png",width=800,height=800) @ <>= persp(x, y, limit(Re(f)), xlab="Re(z)",ylab="Im(z)",zlab="Re(P(z))", theta=30, phi=30, r=1e9, border=NA, shade=0.8, expand=0.6) @ <>= null <- dev.off() @ \begin{figure}[htbp] \begin{center} \includegraphics{wp_figure.png} \caption{Real part of~$\wp(z,1,1+2i)$. Note \label{P.persp.re} the second order poles at each lattice point} \end{center} \end{figure} <>= png("thallerfig.png",width=800,height=800) @ <>= par(pty="s") view(x,y,f,code=0,real.contour=FALSE, imag.contour=FALSE,drawlabel=FALSE,col="red",axes=FALSE,xlab="Re(z)",ylab="Im(z)") axis(1,pos = -4) axis(2,pos = -4) lines(x=c(-4,4),y=c(4,4)) lines(y=c(-4,4),x=c(4,4)) @ <>= null <- dev.off() @ \begin{figure}[htbp] \begin{center} \includegraphics{thallerfig.png} \caption{Visualization of~$\wp(z,1,1+2i)$\label{P.view} using the scheme of \cite{thaller1998}: white corresponds to a pole, black to a zero, and full saturation to~$|\wp(z)|=1$. The poles of~$\wp(z)$ occur on a regular lattice, and the zeros on two shifted lattices. Note how each of the poles is surrounded by two cycles of hue, indicating that they are of second order; and each of the zeros is surrounded by one cycle of hue, indicating that they are simple roots} \end{center} \end{figure} The~$\sigma$ function with the same invariants is visualized in figure~\ref{sigma.green}, showing that its zeros lie on the same lattice as figure~\ref{latplot}. <>= x <- seq(from= -12, to=12, len=n) y <- x z <- outer(x, 1i*y, "+") f <- sigma(z, c(1+1i,2-3i)) @ <>= png("sigma_green.png",width=800,height=800) @ <>= par(pty="s") view(x,y,f,scheme=4,real.contour=FALSE,drawlabels=FALSE,axes=FALSE, xlab="Re(z)",ylab="Im(z)") axis(1,pos= -12) axis(2,pos= -12) lines(x=c(-12,12),y=c(12,12)) lines(y=c(-12,12),x=c(12,12)) lines(x=c(-12,12),y=-c(12,12)) lines(y=c(-12,12),x=-c(12,12)) @ <>= null <- dev.off() @ \begin{figure}[htbp] \begin{center} \includegraphics{sigma_green.png} \caption{Visualization of~$f=\sigma(z,1,1+2i)$ using \code{view()}; colour indicates~${\rm Arg}(f)$. Thus points at which~$f(z)$ is on the negative real axis, that is $\{z\colon f(z)\in\mathbbm{R}^-\}$, are visible as discontinuities of (colourimetric) value. These discontinuities are semi-infinite; note that the zeros of~$f$ occur, \label{sigma.green} at the (finite) end of each line, on a regular lattice. As~$|z|$ increases, each discontinuity threads its way through an increasing number of other discontinuities and zeros, and the spacing between the discontinuities becomes less and less} \end{center} \end{figure} Figure~\ref{zeta.thaller} shows the zeta function, and figure~\ref{sn.thaller} shows Jacobi's ``sn'' function. <>= zeta.z <- zeta(z, c(1+1i,2-3i)) @ <>= png("zetafig.png",width=800,height=800) @ <>= par(pty="s") view(x,y,zeta.z,scheme=0,real.contour=FALSE,drawlabels=FALSE,xlab="Re(z)",ylab="Im(z)") @ <>= null <- dev.off() @ \begin{figure}[htbp] \begin{center} \includegraphics{zetafig.png} \caption{Visualization of~$\zeta(z,1,1+2i)$ using \code{view()} and the colouring scheme of Thaller. Poles appear as white regions, and zeros as black regions. \label{zeta.thaller} Each pole is of single order, each zero is a simple root (one cycle of hue). The poles occur on a lattice; there is no simple pattern to the zeros. Note the overall tendency towards the edges of the square to whiteness: $|f|$ increases with~$|z|$ as per equation~\ref{zeta.definition}} \end{center} \end{figure} <>= jj <- seq(from=-40,to=40,len=n) m <- outer(jj,1i*jj,"+") f <- sn(u=5-2i,m=m,maxiter=1000) @ <>= png("sn_figure.png",width=800,height=800) @ <>= par(pty="s") view(jj,jj,f,scheme=0,r0=1/5,real=T,imag=F,levels=c(0,-0.4,-1),drawlabels=F,xlab="Re(m)",ylab="Im(m)") @ <>= null <- dev.off() @ \begin{figure}[htbp] \begin{center} \includegraphics{sn_figure.png} \caption{Jacobi's ``sn'' function\label{sn.thaller} using the elliptic package. Here, $f={\rm sn}(5-2i,m)$ is visualized, with background utilizing Thaller's scheme, and contours of equal~${\rm Re}(f)$ at three selected values shown as black lines. Note the aperiodic arrangement of poles (white areas) and zeros (black areas)} \end{center} \end{figure} \subsection{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 idealized\footnote{\cite{feynman1966} famously described potential flow as the study of ``dry water''}, 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~\cite{milne1949} for a full exposition. Briefly, we define a {\em complex potential} $w(z)$ to be a complex function \[ w(z)=\phi+i\psi\] and observe that both~$\phi$ and~$\psi$ obey Laplace's equation if~$w$ is differentiable. Given this, we may take the velocity vector~$\mathbf{v}=(v_x,v_y)$ of the fluid to be \[ v_x=\frac{\partial\phi}{\partial x},\qquad v_y=\frac{\partial\phi}{\partial y},\qquad \] and observe that streamlines are given by contours of equal~$\psi$; contours of equal~$\phi$ 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~\ref{z.squared.pot.flow}. This corresponds to a stagnation point, at which the speed vanishes; the streamlines (solid) intersect the equipotential lines (dashed) at right angles. <>= f <- function(z){1i*z^2} x <- seq(from=-6,to=6,len=n) y <- seq(from=-6,to=6,len=n) z <- outer(x,1i*y,"+") @ <>= png("stag_point.png",width=800,height=800) @ <>= par(pty="s") view(x,y,f(z),nlevels=14,imag.contour=TRUE,real.cont=TRUE,scheme=-1, drawlabels=FALSE, axes=FALSE,xlab="Re(z)",ylab="Im(z)") axis(1,pos=-6) axis(2,pos=-6) lines(x=c(-6,6),y=c(6,6)) lines(y=c(-6,6),x=c(6,6)) d1 <- c(-0.1,0,0.1) d2 <- c( 0.1,0,0.1) lines(x=d1,y=1+d2) lines(x=d1,y=-1-d2) lines(x=1-d2,y=d1) lines(x=-1+d2,y=d1) @ <>= null <- dev.off() @ \begin{figure}[htbp] \begin{center} \includegraphics{stag_point.png} \caption{Potential flow on the complex plane: field\label{z.squared.pot.flow} 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} \end{center} \end{figure} 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~$\sum_i m_i\log(z-z_i)$ where the~$i^{\rm th}$ source is at~$z_i$ and has strength~$m_i$, because the system is linear\footnote{It is often more convenient to work with the algebraically equivalent form~$\log\left(\prod (z-z_i)^{m_i}\right)$, as there are fewer branch cuts to deal with.}. Figure~\ref{upper.halfplane.flow} 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\colon \rm{Im}(z)<0\}$) and consider the flow to be bounded below by the real axis. This is known as {\em the method of images}~\citep{milne1949}. <>= f <- function(z){1i*log((z-1.7+3i)*(z-1.7-3i)/(z+1-0.6i)/(z+1+0.6i))} x <- seq(from=-6,to=6,len=n) y <- seq(from=-6,to=6,len=n) z <- outer(x,1i*y,"+") @ <>= png("two_sources_two_sinks.png",width=800,height=800) @ <>= par(pty="s") view(x,y,f(z),nlevels=24,imag.contour=TRUE,real.cont=TRUE,scheme=17,power=0.1,drawlabels=FALSE,axes=FALSE,xlab="Re(z)",ylab="Im(z)") axis(1,pos=-6) axis(2,pos=-6) lines(x=c(-6,6),y=c(6,6)) lines(y=c(-6,6),x=c(6,6)) @ <>= null <- dev.off() @ \begin{figure}[htbp] \begin{center} \includegraphics{two_sources_two_sinks.png} \caption{Potential flow in on the complex plane: two sources and two sinks, all of equal strength. Solid \label{upper.halfplane.flow} lines denote streamlines, dotted lines equipotentials; colour scheme uses the \code{hcl()} system: hue varies with the argument, and chrominance varies with the modulus, of the potential. There is no flux between the lower and the upper half plane, but there is flux out of, and in to, the diagram. Note the stagnation point at approximately $5+0i$} \end{center} \end{figure} Now, one may transform a potential flowfield into another form using a conformal mapping from the~$z$- plane to the~$\zeta$- plane, traditionally denoted \[ \zeta=f(z). \] This technique finds application when flow is desired (in the~$\zeta$- plane) that obeys some specific boundary condition that is simple to specify in the~$z$- plane. %In the present case, we make use of the Schwartz-Christoffel theorem, %which states that if~$a,b,c,\ldots$ are~$n$ points on the real axis of %the~$\zeta$- plane with~$a>= m <- 0.5 K <- K.fun(m) iK <- K.fun(1-m) #b <- sn(1.8 + 0.8i, m=m) # 1.8 to the right and 0.8 up. #b <- 0 # middle bottom b <- sn(0 + 1i*iK/2,m=m) #dead centre of the rectangle. #b <- -1 # lower left #b <- 1/sqrt(m) # top right #b <- -1/sqrt(m) # top left #b <- 1e9*1i # top centre a <- 1 #bottom right hand side corner f <- function(z){1i*log((z-a)*(z-Conj(a))/(z-b)/(z-Conj(b)))} x <- seq(from=-K,to=K,len=n) y <- seq(from=0,to=iK,len=n) z <- outer(x,1i*y,"+") fsn <- f(sn(z,m=m)) @ <>= png("rectangle_pot_flow.png",width=800,height=800) @ <>= view(x,y,fsn,nlevels=44,imag.contour=FALSE,real.cont=TRUE,scheme=17,power=0.1,drawlabels=FALSE,axes=FALSE,xlab="",ylab="") rect(-K,0,K,iK,lwd=3) @ <>= null <- dev.off() @ \begin{figure}[htbp] \begin{center} \includegraphics{rectangle_pot_flow.png} \caption{Potential flow in a rectangle of aspect ratio~2: source and sink of equal \label{box.flow} strength. Colour scheme as in figure~\ref{upper.halfplane.flow}. Note the dividing streamline which terminates in a stagnation point on the rectangle boundary} \end{center} \end{figure} \subsection{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~\citep{johnson2004}. Here, I apply the methods of~\cite{kennedy2001} and~\cite{kennedy2001a} (hereafter KOH and KOHa respectively) using the~\pkg{BACCO} package~\citep{hankin2005} to assess some characteristics of potential flow in a rectangle. Kennedy and O'Hagan considered the following inference problem for a set of parameters~$\theta\in{\mathcal R}^q$ that are inputs to a computer program. Given an independent variable~$x\in{\mathcal R}^n$, and a set of scalar (``field'') observations~$z=z(x)$, they assume \begin{equation} z(x)=\rho\cdot\eta\left(x,\theta\right)+ \delta(x)+\epsilon \end{equation} where~$\rho$ is a constant of proportionality (notionally unity); $\eta(\cdot,\cdot)$ a Gaussian process with unknown coefficients; $\theta$ the true, but unknown parameter values; $\delta(\cdot)$ a model inadequacy term (also a Gaussian process with unknown coefficients); and~$\epsilon\sim N(0,\lambda^2)$ uncorrelated normal observational errors. Inferences about~$\eta(\cdot,\cdot)$ 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~$\mathbf{x}$ is the two-dimensional position within the rectangle, and the field observation~$z(\mathbf{x})$ is the fluid speed at that point, plus obervational error~$\epsilon$. The parameter set~$\theta$ thus has two degrees of freedom corresponding to the $x-$ and $y-$ coordinates of the sink. Field observations will be obtained numerically, using the \pkg{elliptic} package. The simulated flowfield has a sink at a {\em known} position---in this case the geometric centre of the rectangle---and Bayesian methods will be used to infer its position using only fluid speed data. In the terminology of KOHa, dataset~\code{y} corresponds to modelled fluid speed, obtained from the appropriate simulations carried out with the sink placed at different locations within the rectangle. Dataset~\code{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. <>= # Choose the size of the computational mesh: n <- n_BACCO # Choose the number of code observations for D1: n.code.obs <- 30 # And the number of field observations for D2: n.field.obs <- 31 # First, up the D1 design matrix. Recall that D1 is the set of code # observations, which here means the observations of fluid speed when # the sink is at a known, specified, position. suppressWarnings(RNGversion("3.5.0")) set.seed(0) latin.hypercube <- function (n, d){ sapply(seq_len(d), function(...) { (sample(1:n) - 0.5)/n }) } D1.elliptic <- latin.hypercube(n.code.obs , 4) colnames(D1.elliptic) <- c("x","y","x.sink","y.sink") D1.elliptic[,c(1,3)] <- (D1.elliptic[,c(1,3)] -0.5)*2 #D1.elliptic[,c(2,4)] <- D1.elliptic[,c(2,4)] *iK # now a D2 design matrix. This is field observations: observations of # fluid speed when the sink is at the true, unknown, specified, # position. D2.elliptic <- latin.hypercube(n.field.obs , 2) colnames(D2.elliptic) <- c("x","y") D2.elliptic[,1] <- (D2.elliptic[,1] -0.5)*2 # Now a function that, given x and y and x.sink and y.sink, returns # the log of the fluid speed at x,y: fluid.speed <- function(x.pos, y.pos, x.sink, y.sink){ a <- 1 #bottom right hand side corner b <- sn(x.pos/K + 1i*iK*y.pos,m=m) #position (x.pos , y.pos) f <- function(z){1i*log((z-a)*(z-Conj(a))/(z-b)/(z-Conj(b)))} x <- seq(from=-K,to=K,len=n) y <- seq(from=0,to=iK,len=n) z <- outer(x,1i*y,"+") potential <- f(sn(z,m=m)) get.log.ke <- function(x,y,potential){ jj <- Re(potential) jj.x <- cbind(jj[,-1]-jj[,-ncol(jj)],0) jj.y <- rbind(jj[-1,]-jj[-nrow(jj),],0) kinetic.energy <- jj.x^2 + jj.y^2 n.x <- round(n * (x-(-1))/2) n.y <- round(n * y) return(log(kinetic.energy[n.x , n.y]+0.01)) } return(get.log.ke(x.pos,y.pos,potential)) } # now fill in code outputs y: y.elliptic <- rep(NA,nrow(D1.elliptic)) for(i in 1:length(y.elliptic)){ jj <- D1.elliptic[i,,drop=TRUE] y.elliptic[i] <- fluid.speed(jj[1],jj[2],jj[3],jj[4]) } # Now do the field observations; here the source is known to be at the # centre of the rectangle: z.elliptic <- rep(NA,nrow(D2.elliptic)) for(i in 1:length(z.elliptic)){ jj <- D2.elliptic[i,,drop=TRUE] z.elliptic[i] <- fluid.speed(jj[1],jj[2],0,0.5) } # Create design matrix plus observations for didactic purposes: D1 <- round(cbind(D1.elliptic,observation=y.elliptic),2) D2 <- round(cbind(D2.elliptic,observation=z.elliptic),2) # create a data vector: d.elliptic <- c(y.elliptic , z.elliptic) #now a h1.toy() equivalent: h1.elliptic <- function(x){ out <- c(1,x[1]) } #now a H1.toy() equivalent: H1.elliptic <- function (D1) { if (is.vector(D1)) { D1 <- t(D1) } out <- t(apply(D1, 1, h1.elliptic)) colnames(out)[1] <- "h1.const" return(out) } h2.elliptic <- function(x){ c(1,x[1]) } H2.elliptic <- function(D2){ if (is.vector(D2)) { D2 <- t(D2) } out <- t(apply(D2, 1, h2.elliptic)) colnames(out)[1] <- "h2.const" return(out) } #Now an extractor function: extractor.elliptic <- function (D1) { return(list(x.star = D1[, 1:2, drop = FALSE], t.vec = D1[, 3:4, drop = FALSE])) } # Now a whole bunch of stuff to define a phi.fun.elliptic() # and, after that, to call it: phi.fun.elliptic <- function (rho, lambda, psi1, psi1.apriori, psi2, psi2.apriori, theta.apriori, power) { "pdm.maker.psi1" <- function(psi1) { jj.omega_x <- diag(psi1[1:2]) rownames(jj.omega_x) <- names(psi1[1:2]) colnames(jj.omega_x) <- names(psi1[1:2]) jj.omega_t <- diag(psi1[3:4]) rownames(jj.omega_t) <- names(psi1[3:4]) colnames(jj.omega_t) <- names(psi1[3:4]) sigma1squared <- psi1[5] return(list(omega_x = jj.omega_x, omega_t = jj.omega_t, sigma1squared = sigma1squared)) } "pdm.maker.psi2" <- function(psi1) { jj.omegastar_x <- diag(psi2[1:2]) sigma2squared <- psi2[3] return(list(omegastar_x = jj.omegastar_x, sigma2squared = sigma2squared)) } jj.mean <- theta.apriori$mean jj.V_theta <- theta.apriori$sigma jj.discard.psi1 <- pdm.maker.psi1(psi1) jj.omega_t <- jj.discard.psi1$omega_t jj.omega_x <- jj.discard.psi1$omega_x jj.sigma1squared <- jj.discard.psi1$sigma1squared jj.discard.psi2 <- pdm.maker.psi2(psi2) jj.omegastar_x <- jj.discard.psi2$omegastar_x jj.sigma2squared <- jj.discard.psi2$sigma2squared jj.omega_t.upper <- chol(jj.omega_t) jj.omega_t.lower <- t(jj.omega_t.upper) jj.omega_x.upper <- chol(jj.omega_x) jj.omega_x.lower <- t(jj.omega_x.upper) jj.a <- solve(solve(jj.V_theta) + 2 * jj.omega_t, solve(jj.V_theta, jj.mean)) jj.b <- t(2 * solve(solve(jj.V_theta) + 2 * jj.omega_t) %*% jj.omega_t) jj.c <- jj.sigma1squared/sqrt(det(diag(nrow = nrow(jj.V_theta)) + 2 * jj.V_theta %*% jj.omega_t)) jj.A <- solve(jj.V_theta + solve(jj.omega_t)/4) jj.A.upper <- chol(jj.A) jj.A.lower <- t(jj.A.upper) list(rho = rho, lambda = lambda, psi1 = psi1, psi1.apriori = psi1.apriori, psi2 = psi2, psi2.apriori = psi2.apriori, theta.apriori = theta.apriori, power = power, omega_x = jj.omega_x, omega_t = jj.omega_t, omegastar_x = jj.omegastar_x, sigma1squared = jj.sigma1squared, sigma2squared = jj.sigma2squared, omega_x.upper = jj.omega_x.upper, omega_x.lower = jj.omega_x.lower, omega_t.upper = jj.omega_t.upper, omega_t.lower = jj.omega_t.lower, a = jj.a, b = jj.b, c = jj.c, A = jj.A, A.upper = jj.A.upper, A.lower = jj.A.lower) } # OK, that's the function defined. Now to create some jj.* variables # to call it: jj.psi1 <- c(rep(1,4),0.3) names(jj.psi1)[1:4] <- colnames(D1.elliptic) names(jj.psi1)[5] <- "sigma1squared" jj.mean.psi1 <- rep(1,5) names(jj.mean.psi1) <- names(jj.psi1) jj.sigma.psi1 <- diag(0.1,nrow=5) rownames(jj.sigma.psi1) <- names(jj.psi1) colnames(jj.sigma.psi1) <- names(jj.psi1) jj.psi2 <- c(1,1,0.3) names(jj.psi2)[1:2] <- colnames(D2.elliptic) names(jj.psi2)[3] <- "sigma2squared" jj.mean.psi2 <- rep(1,4) names(jj.mean.psi2) <- c("x.sink", "y.sink","rho","lambda") jj.sigma.psi2 <- diag(0.1,4) rownames(jj.sigma.psi2) <- names(jj.mean.psi2) colnames(jj.sigma.psi2) <- names(jj.mean.psi2) jj.mean.th <- c(1,0.5) names(jj.mean.th) <- colnames(D1.elliptic)[3:4] jj.sigma.th <- diag(rep(1,2)) rownames(jj.sigma.th) <- colnames(D1.elliptic)[3:4] colnames(jj.sigma.th) <- colnames(D1.elliptic)[3:4] # Now call phi.fun.elliptic(): phi.elliptic <- phi.fun.elliptic( rho=1, lambda=0.1, psi1=jj.psi1, psi2=jj.psi2, psi1.apriori=list(mean=jj.mean.psi1, sigma=jj.sigma.psi1), psi2.apriori=list(mean=jj.mean.psi2, sigma=jj.sigma.psi2), theta.apriori=list(mean=jj.mean.th, sigma=jj.sigma.th), power=2 ) # Now an E.theta.elliptic(): E.theta.elliptic <- function (D2 = NULL, H1 = NULL, x1 = NULL, x2 = NULL, phi, give.mean = TRUE) { if (give.mean) { m_theta <- phi$theta.apriori$mean return(H1(D1.fun(D2, t.vec = m_theta))) } else { out <- matrix(0, 2,2) rownames(out) <- c("h1.const","x") colnames(out) <- c("h1.const","x") return(out) } } #Now an Edash.theta.elliptic(). Because the basis vector is not a #function of theta, this is a bit academic as we can use a function #that is identical to Edash.theta.toy(): Edash.theta.elliptic <- function (x, t.vec, k, H1, fast.but.opaque = FALSE, a = NULL, b = NULL, phi = NULL) { if (fast.but.opaque) { edash.mean <- a + crossprod(b, t.vec[k, ]) } else { V_theta <- phi$theta.apriori$sigma m_theta <- phi$theta.apriori$mean omega_t <- phi$omega_t edash.mean <- solve(solve(V_theta) + 2 * omega_t, solve(V_theta, m_theta) + 2 * crossprod(omega_t, t.vec[k, ])) } jj <- as.vector(edash.mean) names(jj) <- rownames(edash.mean) edash.mean <- jj return(H1(D1.fun(x, edash.mean))) } # Define a wrapper for equation 8: # First, calculate the constant to subtract to ensure that # the support has a maximum of about zero: maximum.likelihood.support <- p.eqn8.supp(theta=c(0,1/2), D1=D1.elliptic, D2=D2.elliptic, H1=H1.elliptic, H2=H2.elliptic, d=d.elliptic, include.prior=FALSE, lognormally.distributed=FALSE, return.log=TRUE, phi=phi.elliptic) support <- function(x){ p.eqn8.supp(theta=x, D1=D1.elliptic, D2=D2.elliptic, H1=H1.elliptic, H2=H2.elliptic, d=d.elliptic, include.prior=FALSE, lognormally.distributed=FALSE, return.log=TRUE, phi=phi.elliptic) - maximum.likelihood.support } #define a local function called optim() for aesthetic reasons (ie it # improves the appearance of the call to optim(): optim <- function(par,fn){ stats::optim(par=par,fn=fn,control=list(fnscale = -1))$par } @ The code evaluation design matrix~\code{D1} is chosen according to a random Latin hypercube design, and the observation is calculated using the \pkg{elliptic} package: <>= head(D1) @ So the first line shows a simulation with the sink at~(\Sexpr{D1[1,3]},\Sexpr{D1[1,4]}); the log of the fluid speed at~(\Sexpr{D1[1,1]}, \Sexpr{D1[1,2]}) is~\Sexpr{D1[1,5]}. There are a total of~\Sexpr{n.code.obs} such observations. Figure~\ref{code.obs} shows these points superimposed on the ``true'' flow field. The field observations are similarly determined: <>= head(D2) @ showing that the first field observation, at~(\Sexpr{D2[1,1]}, \Sexpr{D2[1,2]}), is~\Sexpr{D2[1,3]}. There are a total of~\Sexpr{n.field.obs} such observations. Figure~\ref{field.obs} shows the first code observation in the context of the ``true'' flow field. <>= b <- sn(D1[1,3] + 1i*D1[1,4],m=m) #point corresponding to first line of D1 fsnz2 <- f(sn(z,m=m)) @ <>= png("code_obs.png",width=800,height=800) @ <>= view(x,y,fsnz2,nlevels=44,imag.contour=FALSE,real.cont=TRUE,scheme=-1,drawlabels=FALSE,axes=FALSE,xlab="",ylab="") points(x=K*D1[1,1],y=D1[1,2]*iK,pch=4) rect(-K,0,K,iK,lwd=3) @ <>= null <- dev.off() @ \begin{figure}[htbp] \begin{center} \includegraphics{code_obs.png} \caption{Streamlines\label{code.obs} of first code observation point; field observation point shown as a cross. The sink is at~(\Sexpr{D1[1,3]},\Sexpr{D1[1,4]})} \end{center} \end{figure} <>= b <- sn(0 + 1i*iK/2,m=m) fsnzz <- f(sn(z,m=m)) @ <>= png("true_flow.png",width=800,height=800) @ <>= view(x,y,fsnzz,nlevels=44,imag.contour=FALSE,real.cont=TRUE,scheme=-1,drawlabels=FALSE,axes=FALSE,xlab="",ylab="") points(x=K*D2[,1],y=D2[,2]*iK,pch=4) rect(-K,0,K,iK,lwd=3) @ <>= null <- dev.off() @ \begin{figure}[htbp] \begin{center} \includegraphics{true_flow.png} \caption{Streamlines\label{field.obs} of ``true'' flow; field observation points shown as crosses} \end{center} \end{figure} Kennedy and O'Hagan give, {\em inter alia,} an expression for the likelihood of any value of $\theta$ 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 \code{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)) #centre of the rectangle support(c(-1,1)) #top left corner @ 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 \code{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) @ \begin{Schunk} \begin{Sinput} (mle <- optim(c(0,1/2),support)) \end{Sinput} \end{Schunk} <>= mle @ 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) @ \subsubsection{Discussion of Bayesian statistical analysis} The above example shows the ideas of KOH being applied straightforwardly, but with the novel twist of $\theta$ being interpreted as physical characteristics of a fluid flow. In this case~$\theta$ 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~\cite{edwards1992}, the hypothesis of $\theta_\mathrm{true}=(0,0.5)$ would not be rejected if one accepted Edwards's 2 units of likelihood per degree of freedom. The discrepancy between~$\hat{\theta}$ and~$\theta_\mathrm{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. \section{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 \proglang{R} package \pkg{elliptic}, which was then used in conjunction with Bayesian statistical methods (the \pkg{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. \subsection*{Acknowledgements} I would like to acknowledge the many stimulating and helpful comments made by the \proglang{R}-help list over the years. \bibliography{elliptic} \end{document}