20.20 ELLIPFN: A Package for Elliptic Functions and
Integrals
Additional documentation and examples can be found in the files efellip.tex,
eftheta.tst and efweier.tst in the packages/ellipfn directory.
Author: Lisa Temme and Alan Barnes, with contributions from Winfried Neun and
several others
20.20.1 Elliptic Functions: Introduction
\(\newcommand {\f }[1]{\texttt {#1}}\)The package ELLIPFN is designed to provide algebraic and numeric manipulations of
many elliptic functions, namely:
The implementation of the functions in this and the next two subsections have been
substantially revised by Alan Barnes in 2019. This is to bring the notation more into line
with standard (British) texts such as Whittaker & Watson [WW69] and Lawden
[Law89] and also to correct a number of errors and omissions. These changes
and additions will be itemised in the relevant sections below. New subsections
has been added starting in 2021 to support Weierstrassian Elliptic functions,
Sigma functions, inverse Jacobi elliptic functions and finally symmetric elliptic
integrals.
The functions in this subsection are for the most part autoloading; the exceptions being
the subsidiary utility functions such as the AGM function, the quasi-period
factors, lattice functions, derivatives of the theta functions and symmetric elliptic
integrals.
20.20.2 Jacobi Elliptic Functions
The following functions have been implemented:
- The 12 Jacobi Elliptic Functions
- The Jacobi Amplitude Function
- Arithmetic Geometric Mean
- Descending Landen Transformation
The following Jacobi elliptic functions are available:-
- jacobisn(u,k)
- jacobidn(u,k)
- jacobicn(u,k)
- jacobicd(u,k)
- jacobisd(u,k)
- jacobind(u,k)
- jacobidc(u,k)
- jacobinc(u,k)
- jacobisc(u,k)
- jacobins(u,k)
- jacobids(u,k)
- jacobics(u,k)
These differ somewhat from the originals implemented by Lisa Temme in that the second
argument is now the modulus (usually denoted by \(k\) in most texts rather than
its square \(m\)). The notation for the most part follows Lawden [Law89]. The last
nine Jacobi functions are related to the three basic ones: jacobisn(u,k),
jacobicn(u,k) and jacobidn(u,k) and use Glaisher’s notation. For example
\[ \mathrm {ns}(x,k) = \frac {1}{\mathrm {sn}(u,k)}, \qquad \mathrm {cs}(x,k) = \frac {\mathrm {cn}(u,k)}{\mathrm {sn}(u,k)}, \qquad \mathrm {cd}(x,k) = \frac {\mathrm {cn}(u,k)}{\mathrm {dn}(u,k)}. \]
All twelve functions are doubly periodic in the complex plane. The primitive
periods and the positions of the zeros and poles of the functions are conveniently
expressed in terms of the so-called quarter periods \(\mathrm {K}(k)\) and \(i\mathrm {K}'(k)\) which are complete
elliptic integrals of the first kind. The details are displayed in the table below; in
all cases the pole is single and the corresponding residue is given in the last
column of the table. As the functions are doubly-periodic, there are, of course,
infinitely many other zeros and poles. These occur at points ‘congruent’ to the
points given in the table obtained by translating by \(2m\mathrm {K} +2i n\mathrm {K}'\) where \(m\) and \(n\) are arbitrary
integers.
| Function | Period1 | Period2 | Zero | Pole | Residue |
| \(\mathop {\mathrm {sn}}\) | \(4\mathrm {K}\) | \(2i\mathrm {K}'\) | \(0\) | \(i\mathrm {K}'\) | \(1/k\) |
| \(\mathop {\mathrm {ns}}\) | \(4\mathrm {K}\) | \(2i\mathrm {K}'\) | \(i\mathrm {K}'\) | \(0\) | \(1\) |
| \(\mathop {\mathrm {cd}}\) | \(4\mathrm {K}\) | \(2i\mathrm {K}'\) | \(\mathrm {K}\) | \(\mathrm {K}+i\mathrm {K}'\) | \(-1/k\) |
| \(\mathop {\mathrm {dc}}\) | \(4\mathrm {K}\) | \(2i\mathrm {K}'\) | \(\mathrm {K}+i\mathrm {K}'\) | \(\mathrm {K}\) | \(-1\) |
| \(\mathop {\mathrm {cn}}\) | \(4\mathrm {K}\) | \(2(\mathrm {K}+i\mathrm {K}')\) | \(\mathrm {K}\) | \(i\mathrm {K}'\) | \(-i/k\) |
| \(\mathop {\mathrm {nc}}\) | \(4\mathrm {K}\) | \(2(\mathrm {K}+i\mathrm {K}')\) | \(i\mathrm {K}'\) | \(\mathrm {K}\) | \(-1/k'\) |
| \(\mathop {\mathrm {sd}}\) | \(4\mathrm {K}\) | \(2(\mathrm {K}+i\mathrm {K}')\) | \(0\) | \(\mathrm {K}+i\mathrm {K}'\) | \(-i/(k k')\) |
| \(\mathop {\mathrm {ds}}\) | \(4\mathrm {K}\) | \(2(\mathrm {K}+i\mathrm {K}')\) | \(\mathrm {K}+i\mathrm {K}'\) | \(0\) | \(1\) |
| \(\mathop {\mathrm {dn}}\) | \(4i\mathrm {K}'\) | \(2\mathrm {K}\) | \(\mathrm {K}+i\mathrm {K}'\) | \(i\mathrm {K}'\) | \(-i\) |
| \(\mathop {\mathrm {nd}}\) | \(4i\mathrm {K}'\) | \(2\mathrm {K}\) | \(i\mathrm {K}'\) | \(\mathrm {K}+i\mathrm {K}'\) | \(-i/k'\) |
| \(\mathop {\mathrm {sc}}\) | \(4i\mathrm {K}'\) | \(2\mathrm {K}\) | \(0\) | \(\mathrm {K}\) | \(-1/k'\) |
| \(\mathop {\mathrm {cs}}\) | \(4i\mathrm {K}'\) | \(2\mathrm {K}\) | \(\mathrm {K}\) | \(0\) | \(1\) |
All other periods of the Jacobi functions can be expressed as linear combinations of the
two primitive periods with integer coefficients. Thus, for example, \(4i\mathrm {K}'\) is a period of \(\mathrm {cn}\) as
\[4i\mathrm {K}' = 2(\mathrm {K}+\mathrm {K}') -4\mathrm {K}\]
Extensive rule lists are provided for differentiation of these functions with respect to
either argument, for argument shifts by integer multiples of the two quarter
periods \(\mathrm {K}(k)\) and \(i\mathrm {K}'(k)\) and finally Jacobi’s transformations for a purely imaginary first
argument.
Rules are also provided for the values of the twelve Jacobi functions at the ‘eighth’
period values: \(\mathrm {K}/2\), \(i\mathrm {K}'/2\) and \((\mathrm {K}+i\mathrm {K}')/2\). For these rules to yield correct values it is essential that the
switch precise_complex is ON, otherwise Reduce will often incorrectly
simplify the results if they involve complex values of \(k\) or \(k'\). It should also be
noted that the rule for \(\mathrm {dn}((\mathrm {K}(k)+i\mathrm {K}'(k))/2, k)\) differs from that on DLMF:NIST which appears to give
incorrect values for some values of \(k\). The rule used in REDUCE is not only
simpler, but appears to give correct values in all cases as do the derived rules
for \(\mathrm {cd}\), \(\mathrm {sd}\), \(\mathrm {nd}\), \(\mathrm {ds}\) and \(\mathrm {dc}\). It is derived from the third of the identities (2.2.19) of Lawden
[Law89] with \(u =-(\mathrm {K}(k)+i\mathrm {K}'(k))/2\) and the corresponding identities for \(\mathrm {sn}\) and \(\mathrm {cn}\) on the DLMF NIST
site.
Some care must be exercised when applying the ‘eighth’ period rules when the modulus \(k\)
belongs to the negative real or imaginary axes as \(\mathrm {K}(k)\) and all twelve Jacobi functions \(\mathrm {sn}(x,k)\) etc. are
even functions of the modulus \(k\). The rule should be applied when the modulus has no
assigned value and then resimplifying the result after assigning the modulus its required
value.
Complex split functions have recently been implemented for the 12 Jacobi functions for
the cases when the modulus \(k\) is either real or imaginary. Thus REPART and IMPART will
now return rational expressions involving the twelve Jacobi functions with real
arguments and a real modulus \(|k|<1\). As far as I (AB) am aware there are no known
methods of implementing split functions for general complex values of the
modulus.
Four useful rule lists TO_SN, TO_CN, TO_DN and NO_GLAISHER are
provided for user convenience. The first TO_SN replaces occurences of the
squares of \(\mathrm {cn}\) and \(\mathrm {dn}\) in favour of the square of \(\mathrm {sn}\) using the ‘Pythagorean’ identities.
TO_CN and TO_DN apply these identities to replace squares in favour of \(\mathrm {cn}\) and \(\mathrm {dn}\)
respectively. At most one of these rule sets should be active at any one time to
avoid a potential infinite recursion in the simplification. The fourth rule list
NO_GLAISHER replaces all occurences of the nine subsidiary Jacobi elliptic functions \(\mathrm {ns}\),
\(\mathrm {sc}\), \(\mathrm {sd}\), …by reciprocals and quotients of the ‘basic’ Jacobi elliptic functions \(\mathrm {sn}\), \(\mathrm {cn}\) and
\(\mathrm {dn}\).
A fifth rule list JACOBIADDITIONRULES applies addition rules when the
first argument of any Jacobi function is a sum. Note that this rule list is NO
LONGER active by default. There are many equivalent ways of expressing the
right-hand sides of these rules. Previously all the rules were expressed using
only the ‘basic’ Jacobi functions: \(\mathrm {sn}\), \(\mathrm {cn}\) and \(\mathrm {dn}\). Currently, however, the right-hand
sides of the rules for the reciprocal and quotient Glaisher functions \(\mathrm {ns}\), \(\mathrm {cs}\), \(\mathrm {sd}\) etc are
expressed using the three Glaisher functions with the same ‘denominator’ as
on the left-hand side. Thus the rule for \(\mathrm {ns}\) is expressed in terms of \(\mathrm {ns}\), \(\mathrm {cs}\) and \(\mathrm {ds}\) whilst
that for \(\mathrm {cd}\) is expressed in terms of \(\mathrm {nd}\), \(\mathrm {cd}\) and \(\mathrm {sd}\) whereas the rules for the three ‘basic’
functions are unchanged and use only \(\mathrm {sn}\), \(\mathrm {cn}\) and \(\mathrm {dn}\). Note, however, that the previous
form of these rules will be obtained if the rule-list NO_GLAISHER is also
active.
If the switch rounded is ON and both arguments of a Jacobi function are purely
numerical, it will be evaluated numerically (previously it was necessary for both the
switches rounded and complex to be ON). Of course, if either argument is a
complex number, the switch complex must also be ON for numerical evaluation to be
triggered.
The numerical evaluation of the Jacobi functions uses their definitions in terms of theta
functions which are valid for all complex values of the argument and modulus.
Since theta functions are inherently complex-valued, it is necessary to turn the
switch complex ON during the evaluation even if both the argument and the
modulus of the Jacobi function are real (for example when the modulus \(|k|>1\)). The
switch complex is, however, returned to its old value as the computation
completes.
The traditional AGM and \(\phi \)-function based algorithms which were used when the
argument and modulus are both real and \(|k| < 1\) have been corrected. There was a problem with
the numerical evaluation of dn, nd, sd, ds, cd, dc resulting in largish
rounding errors for some values of the arguments. There was no problem with previous
evaluation of sn, ns, cn, nc, sc, cs and these produced the same results
(modulo acceptable rounding errors) as the theta function based methods. In fact
in all cases when both the arguments are real and \(|k|<1\), the traditional algorithms
have been reinstated as they are somewhat faster than the theta function based
methods.
Jacobi Amplitude Function
This function is defined as
\[ \mathrm {am}(u,k) = \int _0^u\mathrm {dn}(z,k)\mathrm {d}z.\]
As \(\mathrm {dn}(z,k)\) has poles with residue \(-i\) at points \(z=(4m+1)i\mathrm {K}'(k)+2n\mathrm {K}(k)\) and others with residue \(+i\) when \(z=(4m+3)i\mathrm {K}'(k)+2n\mathrm {K}(k)\) where
\(m\) and \(n\) are arbitrary integers, \(\mathrm {am}(z,k)\) has logarithmic singularities at these points. In
fact the amplitude function is multivalued with its principal value \(v\), say, being
given by choosing the contour in the defining integral to be the straight line
segment between \(0\) and \(u\). Other values are given by \(v+2n\pi \) where \(n\) is an arbitrary integer
and depend on how the chosen contour winds around the logarithmic branch
points. To obtain a single valued function it is necessary to introduce branch
cuts between the points \(i\mathrm {K}'(k)\) and \(3i\mathrm {K}'(k)\) and between corresponding congruent branch
points.
A rule list is provided for to simplify this function for special values of the
arguments and for differentiation with respect to either argument. When the switch
rounded is ON and both arguments are numeric, REDUCE will attempt to evaluate
jacobiam(z,k) numerically. Several possible methods are available to evalaute the
function, for example summation of its Fourier series and an AGM-based method due to
Sala (see DLMF:NIST for more details. Currently the AGM-based method is used and
this certainly produces reliable results if the modulus \(k\) is real and \(|\Im (z)| < \mathrm {K}'(k)\). These agree with those
produced by the Fourier series method although rounding errors become significant as \(|\Im (z)|\)
approaches \(\mathrm {K}'(k)\). If \(k\) is complex both methods produce the same results when the imaginary
parts of \(k\) and \(z\) are not too large, but it is difficult to give precise bounds on the size of
these.
For general values of \(z\) and \(k\) the AGM-based method always produces a value except
perhaps at the logarithmic branch points but these, I believe, are not correct as they fail to
satisfy the identity
\[\mathrm {sn}(z,k) = \sin (\mathrm {am}(z,k)) \]
where \(\mathrm {sn}\) is calculated using the theta function method. The Fourier
series method fails to converge in these cases and so is not used except for comparison
purposes. Currently an alternative method of evaluation due to Sala is being investigated;
it uses the Poisson summation formula and is claimed to be valid for all \(z\) and \(k\) except at
branch points.
20.20.3 Some Numerical Procedures
This section briefly describes several procedures which are primarily intended for use in
the numerical evaluation of the various elliptic functions and integrals rather than for
direct use by users.
Arithmetic Geometric Mean (AGM)
A procedure to evaluate the AGM of initial values \(a_0,b_0,c_0\) exists as
AGM_function(\(a_0,b_0,c_0\)) and will return
\(\{ N, AGM, \{ a_N, \ldots ,a_0\}, \{ b_N, \ldots ,b_0\}, \{c_N, \ldots ,c_0\}\}\), where \(N\) is the number of steps to compute the AGM to the desired accuracy.
To determine the Elliptic Integrals K(\(m\)), E(\(m\)) we use initial values \(a_0 = 1\); \(b_0 = \sqrt {1-k^2}\) ; \(c_0 = k\).
The procedure to evaluate the Descending Landen Transformation of \(\phi \) and \(\alpha \) uses the
following equations:
\begin{align*} (1+\sin \alpha _{n+1})(1+\cos \alpha _n)=2 &\text { where } \alpha _{n+1}<\alpha _n, \\ \tan (\phi _{n+1}-\phi _n)=\cos \alpha _n \tan \phi _n & \text { where } \phi _{n+1}>\phi _n. \end{align*}
It can be called using landentrans(\(\phi _0\), \(\alpha _0\)) and will return
\(\{\{\phi _0, \ldots ,\phi _n\},\{\alpha _0, \ldots ,\alpha _n\}\}\).
20.20.4 Legendre’s Elliptic Integrals
The following functions have been implemented:
- Complete & Incomplete Elliptic Integrals of the First Kind
- Complete & Incomplete Elliptic Integrals of the Second Kind
- Jacobi’s Epsilon & Zeta Functions
- Complete & Incomplete Ellpitic Integrals of the Third Kind
- Some Utility Functions
These again differ somewhat from the originals implemented by Lisa Temme as the
second argument is now the modulus \(k\) rather that its square. Also in the original
implementation there was some confusion between Legendre’s form and Jacobi’s form of
the incomplete elliptic integrals of the second kind; \(E(u,k)\) denoted the first in numerical
evaluations and the second in the derivative formulae for the Jacobi elliptic functions
with respect to their second argument. This confusion was perhaps understandable as in
the literature some authors use the notation \(\mathrm {E}(u, k)\) for the Legendre form and others for Jacobi’s
form.
To bring the notation more into line with that in the NIST Digital Library of
Mathematical Functions and avoid any possible confusion, \(\mathrm {E}(u, k)\) is used for the Legendre form
and \(\mathcal {E}(u, k)\) for Jacobi’s form. This differs from the 2019 version of this section which followed
Lawden [Law89] chapter 3, where the notation \(\mathrm {D}(\phi , k)\) and \(\mathrm {E}(u, k)\) were used for the Legendre and
Jacobi forms respectively.
A number of rule lists have been provided to implement, where appropriate, derivatives
of these functions, addition rules and periodicity and quasi-periodicity properties and to
provide simplifications for special values of the arguments.
Elliptic F
The Elliptic F function can be used as EllipticF(phi,k) and will return the value
of the Incomplete Elliptic Integral of the First Kind:
\[\mathrm {F}(\phi , k)=\int _0^\phi (1-k^2 \sin ^2 \theta )^{-1/2} \mathrm {d}\theta .\]
This is actually closely related to
the inverse Jacobi function \(\mathrm {arcsn}\); in fact for \(-\pi /2 <= \Re \phi <=\pi /2\):
\[ \mathrm {F}(\phi , k) = \mathrm {arcsn}(\sin \phi , k) \]
.
Elliptic K
The Elliptic K function can be used as EllipticK(k) and will return the value
of the Complete Elliptic Integral of the First Kind:
\[\mathrm {K}(k)=\mathrm {F}(\pi /2, k) =\int _0^{\pi /2}(1-k^2 \sin ^2 \theta )^{-1/2}\mathrm {d}\theta .\]
This is one of the quarter
periods of the Jacobi elliptic functions and is often used in the calculation of
other elliptic functions. The complementary Elliptic K\('\) function can be used as
EllipticK!\('\)(k). Note that \(i\mathrm {K}'(k)\) is the other quarter period of the Jacobi elliptic
functions.
Elliptic E
The Elliptic E function comes with either one or two arguments; used with two
arguments as EllipticE(u,k) it will return the value of Legendre’s form of
the Incomplete Elliptic Integral of the Second Kind:
\[\mathrm {E}(\phi , k)=\int _0^\phi \sqrt {1-k^2 \sin ^2 \theta } \,\mathrm {d}\theta .\]
When called with one
argument EllipticE(k) will return the value of the Complete Elliptic Integral
of the Second Kind:
\[\mathrm {E}(k)=\mathrm {E}(\pi /2, k) = \int _0^{\pi /2} \sqrt {1-k^2 \sin ^2 \theta } \,\mathrm {d}\theta .\]
The complementary Elliptic E\('\) function can be used as
EllipticE!\('\)(k).
The complete integrals are actually multi-valued; to obtain single valued functions it is
necessary to introduce branch cuts. For \(\mathrm {K}(k)\) and \(\mathrm {E}(k)\) these are \((-\infty , -1] \cup [1, +\infty )\). The functions are continuous if
the two components of the branch cut are approached from the second and
fourth quadrants respectively. For \(\mathrm {K}'(k)\) and \(\mathrm {E}'(k)\) the branch cut is \((-\infty , 0]\) with continuity if
the cut is approached from the second quadrant. For more details, see Lawden
[Law89] sections 8.12 to 8.14. The principal values of \(\mathrm {K}(k)\) and \(\mathrm {E}(k)\) are even functions of
\(k\).
The numerical evaluation of the complete integrals is more robust and uses symmetric
elliptic integrals. They should now work for all complex values of the modulus and
return the principal value of the integral concerned. Note that for all complex
values of \(k\), the well-known identities:
\[\mathrm {K}'(k)=\mathrm {K}(\sqrt {1-k^2})\qquad \mathrm {E}'(k) = \mathrm {E}(\sqrt {1-k^2}).\]
are not actually valid for the principal
values of the functions concerned. It is necessary that \(\Re k >=0\) with in addition if \(\Re k=0\), \(\Im k>=0\). One
of the consequences of this is that the principal values of \(\mathrm {K}'(k)\) and \(\mathrm {E}'(k)\) are not even
functions.
General values for the four complete integrals together with the principal value when \(\Re {k} <0\) are
given in the table below, where \(n\) is an arbitrary integer, \(k' = \sqrt {1-k^2}\) and the expressions in the second
and third columns refer always to principal values:
| Function | General Value | Principal Value when \(\Re (k) <0\) |
| | | and when \(\Re (k)=0\) and \(\Im (k) <0\) |
| \(\mathrm {K}(k)\) | \(\mathrm {K}(k)-2in\mathrm {K}(k')\) | \(\mathrm {K}(-k)\) |
| \(\mathrm {K}'(k)\) | \(\mathrm {K}(k')-4in\mathrm {K}(k)\) | \(\mathrm {K}(k')\mp 2i\mathrm {K}(-k)\) |
| \(\mathrm {E}(k)\) | \(\mathrm {E}(k)-2in(\mathrm {K}(k')-\mathrm {E}(k'))\) | \(\mathrm {E}(-k)\) |
| \(\mathrm {E}'(k)\) | \(\mathrm {E}(k')-4in(\mathrm {K}(k)-\mathrm {E}(k))\) | \(\mathrm {E}(k')\mp 2i(\mathrm {K}(-k)-\mathrm {E}(-k))\) |
In the third column the upper and lower alternative signs are taken when \(k\) lies in the
second and third quadrants respectively.
One quite subtle point arises: although the twelve Jacobi elliptic functions and \(\mathrm {K}(k)\) are even
functions of the modulus \(k\), the quarter period \(i\mathrm {K}'(k)\) is not (see the second row of the table
above). If \(\Re (k)>0\), then for example \(4\mathrm {K}(k)\) and \(2i\mathrm {K}'(k)=2i\mathrm {K}(k')\) are primitive periods of \(\mathrm {sn}(x,k)\). However, if we use \(-k\) as
the modulus, the first primitive period obtained is the same (since \(\mathrm {K}(k)\) is an even
function) but the second is different namely: \(2i\mathrm {K}'(-k)=2i\mathrm {K}(k') \pm 4\mathrm {K}(k)\). This explains why some of the
values returned by the ‘eighth’ period rules for \(\mathrm {sn}\) etc. are not even functions of
k.
Elliptic D
The Elliptic D function also comes with either one or two arguments; used with two
arguments as EllipticD(u,k) it will return the value of an alternative form of
Legendre’s Incomplete Elliptic Integral of the Second Kind:
\[\mathrm {D}(\phi , k)=\int _0^\phi \frac {\sin ^2 \theta }{\sqrt {1-k^2\sin ^2\theta }} \,\mathrm {d}\theta .\]
When called with
one argument EllipticD(k) will return the value of the Complete Elliptic
Integral of the Second Kind:
\[\mathrm {D}(k)=\mathrm {D}(\pi /2, k) = \int _0^{\pi /2} \frac {\sin ^2 \theta }{\sqrt {1-k^2 \sin ^2 \theta }} \,\mathrm {d}\theta .\]
The integrals of the first and second kind are related:
\[\mathrm {E}(k) = \mathrm {K}(k)-k^2\mathrm {D}(k), \qquad \mathrm {E}(\phi , k) = \mathrm {F}(\phi , k)-k^2\mathrm {D}(\phi , k).\]
For all the elliptic integrals, rule lists are provided for simplification for special values
of the argument(s), differentiation with respect to either argument and shift
rules for the incomplete integrals so that their first argument lies in the range
\(0 <= \phi <= \pi /2\).
If the arguments of the elliptic integral function are numeric and the switch rounded is
ON, the numerical value will be returned. As with Jacobi elliptic functions, if any
argument is complex then numerical evaluation will only occur if the switch
complex is also ON.
Numerical evaluation should now succeed for all complex values of the arguments
provided, of course, the corresponding integral exists. The incomplete integrals \(\mathrm {F}(\phi , k)\), \(\mathrm {E}(\phi , k)\), \(\mathrm {D}(\phi , k)\) and \(\Pi (\phi , a, k)\)
are all multi-valued as there are branch points in the defining integrands when \(\sin (\phi ) =\pm 1\) and
\(k\sin (\phi ) = \pm 1\).
For numerical evaluation when the first argument \(\phi \) is complex there are issues that
arise which are absent when \(\phi \) is real. One might use the the straight line contour
from \(\theta =0\) to \(\theta =\phi \), but this is leads to integrals that are difficult to evaluate. When \(|\Re {\phi }| \leq \pi /2\), the
change of variable \(x=\sin \theta \) is applied and the value returned uses the straight line contour
between \(x=0\) and \(x=\sin (\phi )\). This will produce a different value if the loop formed by the two
contours encloses a branch point. When \(|\Re {\phi }| > \pi /2\) the contour used is the straight line
segment along the real axis from 0 to the multiple of \(\pi \) nearest to \(\phi \) followed by the
straight line segment from between \(x=0\) and \(x=\sin (\phi )\). When \(\phi \) is real the contours in the
\(\theta \) and \(x\) planes coincide. Other contours between the two endpoints will yield
different values depending on how the contour winds around the two branch
points.
Note also that if the contour contains a branch point, there is sign ambiguity as the
branch point is crossed depending on which branch of the square root in the integrand is
chosen. It seems logical to choose the principal branch.
For example when \(k>1\) & \(1/k < y <1\)
\[ \int _0^y \frac {\mathrm {d}x}{\sqrt {1-x^2}\sqrt {1-k^2x^2}}\]
is taken to be
\[\int _0^{1/k} \frac {\mathrm {d}x}{\sqrt {1-x^2}\sqrt {1-k^2x^2}} -i\int _{1/k}^y \frac {\mathrm {d}x}{\sqrt {1-x^2}\sqrt {k^2x^2-1}}.\]
Elliptic \(\Pi \)
The Elliptic \(\Pi \) function can be used as EllipticPi( ) and will return the value of the
Elliptic Integral of the Third Kind.
The Elliptic \(\Pi \) function comes with either two or three arguments; when called with three
arguments EllipticPi(phi,a,k) will return the value of the Incomplete Elliptic
Integral of the Third Kind:
\[\Pi (\phi ,a,k)=\int _0^{\phi } \left ((1-a^2\sin ^2 \theta )\sqrt {1-k^2 \sin ^2 \theta }\right )^{-1} \,\mathrm {d}\theta .\]
For \(-\pi /2 \leq \Re {\phi } \leq \pi /2\), the incomplete integral may be written in the form:
\[\Pi (\phi ,a,k)=\int _0^{sin \phi } \left ((1-a^2t^2)\sqrt {(1-t^2)(1-k^2t^2)}\right )^{-1}\,\mathrm {d}t.\]
When called with two arguments as EllipticPi(a,k) it will return the value of the
Complete Elliptic Integral of the Third Kind:
\begin{eqnarray*} \Pi (a, k)= \Pi (\pi /2,a,k) & = & \int _0^{\pi /2} \left ((1-a^2\sin ^2 \theta )\sqrt {1-k^2 \sin ^2 \theta }\right )^{-1} \,\mathrm {d}\theta \\ & =& \int _0^1 \left ((1-a^2t^2)\sqrt {(1-t^2)(1-k^2t^2)}\right )^{-1}\,\mathrm {d}t. \end{eqnarray*}
For certain values of \(a\) namely \(0\), \(\pm 1\) and \(\pm k\) the integrals reduce to elliptic integrals of the first
and second kinds.
20.20.5 Jacobi’s Elliptic Integrals
Jacobi E
The Jacobi E function can be used as jacobiE(u,k); it will return the
value of Jacobi’s form of the Incomplete Elliptic Integral of the Second Kind:
\[\mathcal {E}(u, k)=\int _0^u \mathrm {dn}^2 (v, k) \,\mathrm {d}v.\]
The relationship between the two forms of incomplete elliptic integrals can be expressed
as
\[\mathcal {E}(u, k) = \mathrm {E}(\mathrm {am}(u), k).\]
Note that
\[\mathrm {E}(k)=\mathcal {E}(\mathrm {K}(k), k) =\int _0^{\mathrm {K}(k)} \mathrm {dn}^2(v, k) \,\mathrm {d}v.\]
On a GUI that supports calligraphic characters (NB. this is now the case with the CSL
GUI), there is no problem and it is rendered as \(\mathcal {E}(u,k)\) in accordance with NIST usage. On
non-GUI interfaces the Jacobi E function is rendered as E_j.
Jacobi’s Zeta Function
This can be called as jacobiZeta(u,k) and refers to Jacobi’s (elliptic) Zeta function
\(\mathrm {Z}(u,k)\) whereas the operator Zeta will invoke Riemann’s \(\zeta \) function. It is closely related to
Jacobi’s Epsilon function jacobie; in fact
\[\mathrm {Z}(u,k) = \mathcal {E}(u,k)-u\mathrm {E}(k)/\mathrm {K}(k).\]
20.20.6 Symmetric Elliptic Integrals
The functions in this section are currently only intended for use for in the numerical
evaluation of the ‘classical’ elliptic integrals discussed above and in certain other ‘basic’
elliptic integrals. The switch rounded should be ON.
Symmetric elliptic integrals are a relatively new development in the theory of elliptic
functions mainly due to Carlson and his collaborators; an extensive bibliography is
available on the NIST Digital Library of Mathematical Functions starting at the link
DLMF:NIST, Bibliography. They have a number of advantages over more traditional
approaches; see for example Advantages of Symmetry.
The fundamental symmetric integrals are all integrals over the positive real line and
involve the function \(s(t)=\sqrt {t+x}\sqrt {t+y}\sqrt {t+z}\) where \(x,y,z\in \mathbb {C} \setminus (-\infty , 0]\) but, except where otherwise stated, at most one of \(x,y,z\) may be
zero.
\begin{eqnarray*} \mathrm {RF}(x,y,z) & = & \frac {1}{2}\int _0^\infty \frac {\mathrm {d}t}{s(t)}\\ \mathrm {RG}(x,y,z) & = & \frac {1}{4}\int _0^\infty \frac {1}{s(t)}\left ( \frac {x}{t+x}+\frac {y}{t+y}+\frac {z}{t+z}\right )t\mathrm {d}t\\ \mathrm {RJ}(x,y,z,p) & = & \frac {3}{2}\int _0^\infty \frac {\mathrm {d}t}{s(t)(t+p)}\\ \mathrm {RD}(x,y,z) & = & \frac {3}{2}\int _0^\infty \frac {\mathrm {d}t}{s(t)(t+z)} \qquad z \neq 0, \quad x+y \neq 0\\ \mathrm {RC}(x,y) & = & \frac {1}{2}\int _0^\infty \frac {\mathrm {d}t}{\sqrt {t+x)}(t+y)}\qquad y \neq 0 \end{eqnarray*}
The first three integrals defined above are symmetric in \(x,y,z\) and are termed symmetric
elliptic integrals of the first, second and third kinds respectively. \(\mathrm {RD}\) and \(\mathrm {RC}\) are degnerate
versions of \(\mathrm {RJ}\) and \(\mathrm {RF}\) respectively as
\[\mathrm {RD}(x,y,z) =\mathrm {RJ}(x,y,z,z) \qquad \mathrm {RC}(x,y) =\mathrm {RF}(x,y,y) \]
\(\mathrm {RD}\) is an elliptic integral of the second kind and is only
symmetric in \(x,y\) whilst \(\mathrm {RC}\) is an elementary integral, but can be numerically evaluated
efficiently without needing to distinguish between the trigonometric and hyperbolic
cases. For certain special values of \(p\) the integral \(\mathrm {RJ}\) of the third kind degenerates into
integrals of the first and second kinds. For numeric evaluations it turns out that the use of \(\mathrm {RD}\)
is more convenient than that of \(\mathrm {RG}\) (which is not currently implemented in REDUCE). For
more information see §19.15-19.29 of the DLMF:NIST.
The Legendre elliptic integrals may all be expressed in terms of the symmetric integrals;
the complete integrals satisfy
\begin{eqnarray*} \mathrm {K}(k) &=& \mathrm {RF}(0,k^{\prime 2},1) \\ \mathrm {D}(k) &=& \mathrm {RD}(0,k^{\prime 2},1)/3 \\ \mathrm {E}(k) &=& \mathrm {RF}(0,k^{\prime 2},1)-k^2\mathrm {RD}(0,k^{\prime 2},1)/3\\ \Pi (a,k) &=&\mathrm {RF}(0,k^{\prime 2},1)+a^2\mathrm {RJ}(0,k^{\prime 2},1,1-a^2)/3 \end{eqnarray*}
where here and below \(k^\prime =\sqrt {1-k^2}\). For the incomplete integrals, defining \(s=\sin \phi \) with \(-\pi /2 \leq \Re {\phi } \leq \pi /2\), we have
\begin{eqnarray*} \mathrm {F}(\phi ,k) &=& s \mathrm {RF}(1-s^2, 1-k^2s^2, 1)\\ \mathrm {D}(\phi ,k) &=& s^3 \mathrm {RD}(1-s^2, 1-k^2s^2, 1)/3\\ \mathrm {E}(\phi ,k) &=&\mathrm {F}(\phi ,k)-k^2s^3\mathrm {RD}(1-s^2,1-k^2s^2,1)/3\\ \Pi (\phi ,a,k) &=&\mathrm {F}(\phi ,k) +a^2s^3\mathrm {RJ}(1-s^2,1-k^2s^2,1,1-a^2s^2)/3 \end{eqnarray*}
The above formulae are only valid when the range of integration does not include one or
more branch points of the integrand (when \(s\) is real and \(s^2>1\) or when \(k s\) is real and \((k s)^2 >1\)). In these
cases it is necessary to split the range of integration at the branch points resulting in two
or three elliptic integrals.
Currently the Carlson’s duplication method is used to evaluate the symmetric elliptic
integrals of the first, second and third kinds \(\mathrm {RF}\), \(\mathrm {RD}\) and \(\mathrm {RJ}\) (and also the related elementary
integral \(\mathrm {RC}\)). For more details see the DLMF website: Duplication Method. In particular the
REDUCE code is essentially a translation from Fortran of the code of Carlson &
Notis [CN81] generalised for complex arguments and structured to avoid GO
TO.
Alternatively the symmetric integral \(\mathrm {RF}\) may be evaluated using a sequence of quadratic
transformations which converge rapidly to the elementary hyperbolic integral:
\[\mathrm {RC}(X^2+Y^2,X^2) = \mathrm {RF}(X^2,X^2,X^2+Y^2) = \mathrm {arctanh} (Y/(X^2+Y^2))/Y.\]
More
information on this method due to Carlson in the 1990’s may be found on the DLMF
website: Quadratic Transformations.
Basic Elliptic Integrals
Elliptic integrals involving the square root \(s(t)\) of a cubic or quartic polynomial in which the
range of integration is an arbitrary interval \([y,x]\) of the real line are considered. Again there
are three basic kinds: first, second and third. Carlson [Car88] & [Car99] showed how
each of these can be expressed as a fundamental symmetric elliptic integral of the
corresponding kind over the range \([0,\infty )\). Let
\begin{eqnarray*} X_\alpha = & \sqrt {a_\alpha +b_\alpha x} \quad & 1 \leq \alpha \leq 5\\ Y_\beta = & \sqrt {a_\beta +b_\beta y} \quad & 1 \leq \beta \leq 5\\ d_{\alpha \beta } = & a_\alpha b_\beta - a_\beta b_\alpha \quad & d_{\alpha \beta } \neq 0\mbox { if } \alpha \neq \beta \\ s(t) = & \prod ^4_{\alpha =1} \sqrt {a_\alpha +b_\alpha t} \end{eqnarray*}
and where the four line segments with end points \(a_\alpha +b_\alpha y\) and \(a_\alpha +b_\alpha x\) for \(1\leq \alpha \leq 4\) lie in \( \mathbb {C} \setminus (-\infty , 0)\). Note that if
\(s(t)\) is a cubic then one simply chooses supplies unity i.e. \(a_\alpha =1, b_\alpha =0\) as a fourth factor of
\(s(t)\).
Then integrals of the first and second kinds take the form:
\begin{eqnarray*} \int _y^x \frac {\mathrm {d}t}{s(t)} &=& 2\mathrm {RF}(U_{12}^2, U_{13}^2, U_{23}^2)\\ \int _y^x\frac {(a_1+b_1t)\mathrm {d}t}{(a_4+b_4t)s(t)} &=& \frac {2}{3}d_{12}d_{14} \mathrm {RD}(U_{12}^2, U_{13}^2, U_{23}^2) +\frac {2X_1Y_1}{X_4Y_4U_{23}}\qquad \mbox {if }U_{23}\neq 0 \end{eqnarray*}
where in the second equation we have chosen (wlog) the distinguished factors to have
indices 1 & 4.
For any permutation \(\alpha ,\beta ,\gamma ,\delta \) of 1,2,3,4, let
\begin{eqnarray*} U_{\alpha \beta } & = &(X_\alpha X_\beta Y_\gamma Y_\delta + X_\gamma X_\delta Y_\alpha Y_\beta ) /(x-y) \qquad \mbox {for } x,y \mbox { finite}\\ & = & \sqrt {b_\alpha }\sqrt {b_\beta }Y_\gamma Y_\delta + Y_\alpha Y_\beta \sqrt {b_\gamma }\sqrt {b_\delta } \qquad \qquad \mbox {for } x = \infty \\ & = & X_\alpha X_\beta \sqrt {-b_\gamma }\sqrt {-b_\delta }+X_\gamma X_\delta \sqrt {-b_\alpha }\sqrt {-b_\beta }X_\gamma \qquad \mbox {for } y = -\infty \end{eqnarray*}
Clearly \(U_{\alpha \beta }=U_{\beta \alpha }=U_{\gamma \delta }=U_{\delta \gamma }\) and at most one of \(U_{12}, U_{13}, U_{23}\) is zero since \(U_{\alpha \beta }^2 - U_{\alpha \gamma }^2 = d_{\alpha \delta }d_{\beta \gamma } \neq 0\), thus at most one of the parameters of \(\mathrm {RF}\) and \(\mathrm {RD}\) is
zero. If \(X_4\) or \(Y_4\) is zero,that is if \(a_4+b_4t\) vanishes at one end of the range of integration, the integral
of the second kind diverges.
In the ‘awkward’ case when \(U_{23}=0\) or \(U_{12}^2+U_{13}^2=0\), the above method breaks down. However it is still
possible to calculate the required integral; choose \(\beta =2\) or \(\beta =3\) so that \(b_\beta \neq 0\) and note that
\begin{eqnarray*} b_\beta \int \frac {(a_1+b_1t)\mathrm {d}t} {(a_4+b_4t)s(t)} -b_1\int \frac {(a_\beta +b_\beta t)\mathrm {d}t}{(a_4+b_4t)s(t)} &=& d_{1\beta }\int \frac {\mathrm {d}t}{(a_4+b_4t)s(t)} \\ b_4\int \frac {(a_\beta +b_\beta t)\mathrm {d}t}{(a_4+b_4t)s(t)} -b_\beta \int \frac {(a_4+b_4 t)\mathrm {d}t}{(a_4+b_4t)s(t)} &=& d_{\beta 4}\int \frac {\mathrm {d}t}{(a_4+b_4t)s(t)}. \end{eqnarray*}
The second term on the lhs of the second equation reduces to \(b_\beta \int 1/s(t)\mathrm {d}t\) and so is an integral of the
first kind. Thus the required integral can be expressed in terms of an integral of the
second kind and one of the first kind.
Integrals of the third kind take the form
\[\int _y^x\frac {(a_1+b_1t)\mathrm {d}t}{(a_5+b_5t)s(t)} = \frac {2d_{12}d_{13} d_{14}}{3d_{15}}\mathrm {RJ}(U_{12}^2, U_{13}^2, U_{23}^2,U_{15}^2) + \mathrm {RC}(S_{15}^2,Q_{15}^2)\]
where
\begin{eqnarray*} S_{15} &=& \left (\frac {X_2X_3X_4}{X_1}Y_5^2+\frac {Y_2Y_3Y_4}{Y_1}X_5^2\right )/(x-y) \qquad \mbox {for }x,y\mbox { finite}\\ &=& \frac {X_2X_3X_4}{X_1}Y_5^2+\frac {Y_2Y_3Y_4}{Y_1}X_5^2 \quad \qquad \mbox {for } x = \infty \mbox { or } y = -\infty \\ U_{15}^2 &=& U_{1\beta }^2-\frac {d_{1\gamma }d_{1\delta }d_{\beta 5}}{d_{15}} = U_{\beta \gamma }^2-\frac {d_{1\beta }d_{1\gamma }d_{\delta 5}}{d_{15}} \neq 0\\ Q_{15}^2 &=& \frac {(X_5Y_5)^2}{(X_1Y_1)^2}U_{15}^2\qquad \mbox { where } S_{15}^2-Q_{15}^2 = \frac {d_{25}d_{35}d_{45}}{d_{15}} \neq 0 \end{eqnarray*}
where \(\beta ,\gamma ,\delta \) is a permutation of 2,3,4. Again at most one of the parameters of \(\mathrm {RJ}\) is zero. This
method will break down when calculating \(S_{15}\) if \(a_1+b_1t\) vanishes at either the upper or lower limit
of integration as either \(X_1\) or \(Y_1\) is zero. However the situation can be remedied by choosing \(\beta \)
so that \(X_\beta \neq 0\), \(Y_\beta \neq 0\) and \(b_\beta \neq 0\) and using similar methods to that used for the awkward case for integrals
of the second kind.
Note if the integration range is \((-\infty , +\infty )\), the above formulae for integrals of all three kinds are not
valid and the integral needs to be split into two at, say, zero. Similarly if the integration
range is such that the integrand has branch points, the integration range will need to be
split into two at each of these branch points.
In REDUCE elliptic integrals of the three kinds may be evaluated numerically when the
switch rounded is ON by calling the functions ellint_1st, ellint_2nd and
ellint_3rd. The first two parameters are the lower and upper limits of the
integration range; these are followed by 4 (or 5 for ellint_3rd) two-element lists
\(\{a_1,b_1\}, \{a_2,b_2\} \ldots \).
Known bug: As pointed out by Carlson & FitzSimmons [CF00], the algorithms for \(\mathrm {RF}\) & \(\mathrm {RJ}\)
may break down when \(s(t)\) is the product of two quadratic factors each with complex
conjugate zeros \(x_1\pm i y_1\) & \(x_2\pm i y_2\). One of the three quantities \(U_{12}, U_{13}, U_{23}\) may be negative and this causes the
algorithm to return an incorrect result. This error only arises when the crossing point
(neccessarily real) of the diagonals of the parallellogram in the complex plane formed by
the four zeros of \(s(t)\) lies inside the range of integration. The algorithm also appears to
produce the correct result when the crossing point is inside the range of integration but is
‘sufficiently close’ to either end of the range – the precise condition appears to be
unknown. In the same paper Carlson & FitzSimmons propose a modified algorithm of
the same ilk to overcome these difficulties, however it is yet to be implemented in
REDUCE.
Some Examples
In this subsection some of the advantages of symmetry are illustrated; the formula for the
integral of the first kind replaces the 28 formulae in Gradshteyn and Ryzhik. For the
cubic case one of the factors in \(s(t)\) is simply taken to be unity whilst for cases of the form
\[\int _\alpha ^\beta \frac {\mathrm {d}x}{\sqrt {(a+b x^2)(c+d x^2)}}\]
one must use the substitution \(t=x^2\) to obtain
\[\frac {1}{2}\int _{\alpha ^2}^{\beta ^2} \frac {\mathrm {d}t}{\sqrt {t(a+b t)(c+d t)}}\]
Moreover the restriction that one limit of
integration be a branch point of the integrand is eliminated without doubling the number
of standard integrals in the result. Similarly the formula for the integral of the second
kind replaces no less than 144 cases in Gradshteyn and Ryzhik (72 each for the quartic
and cubic cases). For cases where the numerator in the integrand is unity one simply
takes \(a_1=0, b_1=1\) and if there is to be no apparent term of the form \((a_4+b_4t)^{3/2}\) in its denominator one takes
\(a_4=0, b_4=1\).
When \(0 \leq \phi \leq \pi /2\), the fundamental integrals of the first, second and third kinds \(\mathrm {F}(\phi ,k)\), \(\mathrm {E}(\phi ,k)\), \(\mathrm {D}(\phi ,k)\) and \(\Pi (\phi ,a,k)\) may be
written as
\begin{eqnarray*} \mathrm {F}(\phi ,k) &=& \frac {1}{2}\int _0^{s^2} \frac {\mathrm {d}t}{\sqrt {t(1-t)(1-k^2t)}}\\ \mathrm {E}(\phi ,k) &=& \frac {1}{2}\int _0^{s^2} \frac {\sqrt {1-k^2t}\mathrm {d}t}{\sqrt {t(1-t)}}\\ \mathrm {D}(\phi ,k) &=& \frac {1}{2}\int _0^{s^2} \frac {\sqrt {t}\mathrm {d}t}{\sqrt {(1-t)(1-k^2t)}}\\ \Pi (\phi ,a,k) &=& \frac {1}{2}\int _0^{s^2} \frac {\mathrm {d}t}{\sqrt {t(1-t)(1-k^2t)}(1-a^2t)} \end{eqnarray*}
where in all four cases use has been made of the substitution \(t=\sin ^2 \theta \) and \(s=\sin \phi \). In the last three
equations one takes \(a_4=1, b_4 =0\) so that the fourth factor in \(s(t)\) is unity and in the fourth equation one
also takes \(a_1=1,b_1=0\).
20.20.7 Some Numerical Utility Functions
Five utility functions are provided:
nome2mod(q)
nome2mod!\('\)(q)
nome2!K(q)
nome2!K!\('\)(q)
nome(k)
These are only operative when the switch rounded is on and their argument is
numerical. The first pair relate the nome \(q\) of the theta functions with the moduli \(k\) and \(k'=\sqrt {1-k^2}\) of
the associated Jacobi elliptic functions.
The second pair return the quarter periods K and K\('\) respectively of the Jacobi elliptic
functions associated with the nome \(q\).
Finally, nome(k) returns the nome \(q\) associated with the modulus \(k\) of a Jacobi elliptic
function and is essentially the inverse of nome2mod.
20.20.8 Jacobi Theta Functions
These theta functions differ from those originally defined by Lisa Temme in a number
of respects. Firstly four separate functions of two arguments are defined:
elliptictheta1(u,tau)\(\qquad \vartheta _1(u, \tau )\)
elliptictheta2(u,tau)\(\qquad \vartheta _2(u, \tau )\)
elliptictheta3(u,tau)\(\qquad \vartheta _3(u, \tau )\)
ellipticthetas(u,tau)\(\qquad \vartheta _4(u, \tau )\)
rather than a single function with three arguments (with the first argument taking integer
values in the range 1 to 4). Secondly the periods are \(2\pi , 2\pi , \pi \) and \(\pi \) respectively (NOT 4K, 4K, 2K
and 2K). Thirdly the second argument is the modulus \(\tau = a+i b\) where \(b=\Im \tau >0\) and hence the quasi-period
of the theta functions is \(\pi \tau \).
The second parameter was previously the nome \(q=e^{i\pi \tau }\) where \(|q|<1\) since \(\Im \tau >0\). As a consequence
elliptictheta1 and elliptictheta2 were multi-valued owing to the appearance
of \(q^{1/4}\) in their defining expansions. elliptictheta3 and elliptictheta4 were,
however, single-valued functions of \(q\).
Regarded as functions of \(\tau \), elliptictheta1 and elliptictheta2 are
single-valued functions. The nome is given by \(q = \exp (i\pi \tau )\) so that the condition \(\Im (\tau )>0\) ensures that
\(|q| < 1\). Note also in this case \(q^{1/4}\) is interpreted as \(\exp (i\pi \tau /4)\) rather than the principal value of \(q^{1/4}\).
Thus, \(\tau \), \(2+\tau \), \(4+\tau \) and \(6+\tau \) produce four different values of both elliptictheta1 and
elliptictheta2 although they all correspond to the same nome \(q\).
The four theta functions are defined by their Fourier series:
\begin{align*} \vartheta _1(z,\tau ) & = 2 e^{i\pi \tau /4}\sum _{n=0}^\infty (-1)^nq^{n^2+n} \sin (2n+1)z\\ \vartheta _2(z,\tau ) & = 2 e^{i\pi \tau /4}\sum _{n=0}^\infty q^{n^2+n} \cos (2n+1)z\\ \vartheta _3(z,\tau ) & = 1 +2\sum _{n=1}^\infty q^{n^2} \cos 2n z\\ \vartheta _4(z,\tau ) & = 1 +2\sum _{n=1}^\infty (-1)^n q^{n^2} \cos 2n z. \end{align*}
Utilising the periodicity and quasi-periodicity of the theta functions some generalised
shift rules are implemented to shift their first argument into the base period parallelogram
with vertices
\[(\pi /2, \pi \tau /2),\quad (-\pi /2, \pi \tau /2),\quad (-\pi /2, -\pi \tau /2), \quad (\pi /2, -\pi \tau /2).\]
Together with the relation \(\vartheta _1(0,\tau )=0\), these shift rules serve to simplify all four theta
functions to zero when appropriate.
Numerical Evaluation
When the switch rounded is ON and the arguments are purely numerical, the theta
functions are evaluated numerically. Note that as the modulus \(\tau \) is necessarily
complex with a positive imaginary part, the switch complex should normally be
ON.
However, there is one exception if the second parameter is purely real and lies between 0
and 1, it is interpreted as the nome \(q\) rather than the modulus. This exception facilitates the
evaluation and graphing of theta functions in the case when both the argument and the
nome are real which occurs frequently in practical applications. Theta functions are not
defined if the second parameter has a negative imaginary part or if it is real and \(<= 0\) or \(>=1\) and
the numerical evaluation fails.
In what follows \(a\) and \(b\) will denote the real and imaginary parts of \(\tau \) respectively and
so \(|q| = \exp (-\pi b)\) and \(\arg q =\pi a\). The series for the theta functions are fairly rapidly convergent due
to the quadratic growth of the exponents of the nome \(q\) – except for values of \(q\)
for which \(|q|\) is near to 1 (i.e. \(b=\Im \tau \) close to zero). In such cases the direct algorithm
would suffer from slow convergence and rounding errors. For such values of \(|q|\),
Jacobi’s transformation \(\tau '=-1/\tau \) can be used to produce a smaller value of the nome and
so increase the rate of convergence. This works very well for real values of \(q\),
or equivalently for \(\tau \) purely imaginary since \(q'= q^{1/b^2}\), but for complex values the gains
are somewhat smaller. The Jacobi transformation produces a nome \(q'\) for which
\(|q'| = |q|^{1/(a^2+b^2)}\).
When \(\Re q < 0\), the Jacobi transformation is preceded by either the modular transformation \(\tau ' = \tau +1\) when
\(\Im q < 0\), or \(\tau ' = \tau -1\) when \(\Im q > 0\), which both have the effect of multiplying \(q\) by \(-1\), so that the new nome has a
non-negative real part and \(|a| \leq 1/2\). Thus the worst case occurs for values of the nome \(q\) near to \(\pm i\)
where \(|q'| \approx |q|^4\).
By using a series of Jacobi transformations preceded, if necessary by \(\tau \)-shifts to ensure \(|a| <= 1/2\), \(|q|\)
may be reduced to an acceptable level. Somewhat arbitrarily these Jacobi’s
transformations are used until \(b > 0.6\) (i.e. \(|q| < 0.15\)). This seems to produce reasonable behaviour.
In practice more than two applications of Jacobi transformations are rarely
necessary.
The previous version of the numerical code returned the principal values of \(\vartheta _1\) and \(\vartheta _2\), that is
the ones obtained by taking the principal value of \(q^{1/4}\) in their series expansions. The current
version replaces \(q^{1/4}\) by \(\exp (i\pi \tau /4)\). If the principal value is required, it is easily obtained by
multiplying by the ‘correcting’ factor \(q^{1/4}\exp (-i\pi \tau /4)\).
Derivatives of Theta Functions
Four functions are provided:
theta1d(u,ord,tau)
theta2d(u,ord,tau)
theta3d(u,ord,tau)
theta4d(u,ord,tau)
These return the \(d\)th derivatives of the respective theta functions with respect to
their first argument \(u\); the third parameter \(\tau \) is as for theta functions; either the
modulus if \(\Im \tau >0\) or the nome if \(\tau \) is real and lies between 0 and 1. These functions
are only triggered when the switch rounded is ON and their arguments are
numeric with \(d\) being a positive integer. They are provided mainly to support the
implementation the Weierstrassian and Sigma functions discussed in the following
subsection.
The numeric code simply sums the Fourier series for the required derivatives. Unlike the
theta functions themselves the code does not use the quasi-periodicity nor modular
transformations to speed up the convergence of the series by reducing the sizes of \(\Im u\)
and \(|q|\). In the numerical evaluation of the Weierstrassian and Sigma functions
these functions are only called after the necessary shifts of the argument \(u\) and
modular transformations of \(\tau \) have been performed. These are much simpler in this
context.
Nevertheless they may be used from top level and numerical experiments reveal that the
rounding errors are not significant provided \(|q|\) is not near one (say \(|q|<0.9\)) and \(u\) is real or at least
has a relatively small imaginary part.
20.20.9 Weierstrass Elliptic & Sigma Functions
Three main functions of three arguments are defined:
- \(\wp (u, \omega _1, \omega _3)\) —
weierstrass(u,omega1,omega3)
- \(\zeta _w(u, \omega _1, \omega _3)\) —
weierstrassZeta(u,omega1,omega3)
- \(\sigma (u, \omega _1, \omega _3)\) —
weierstrass_sigma(u,omega1,omega3)
The notation used is broadly similar used by Lawden [Law89] which is also used in the
NIST Digital Library of Mathematical Functions DLMF:NIST. However, \(\zeta _w\) is used
for the Weierstrassian Zeta function to distinguish it from the Riemann Zeta
function and the usual symbol \(\wp \) is used for the Weierstrassian elliptic function
itself.
The two primitive periods of the Weierstrass function are \(2\omega _1\) and \(2\omega _3\) and these must satisfy \(\Im (\omega _3/\omega _1) \neq 0\).
The two periods are normally numbered so that \(\tau = \omega _3/\omega _1\) has a positive imaginary part and hence
the nome \(q = exp(i\pi \tau )\) satisfies \(|q| <1\).
Any linear combination \(\Omega _{m,n} = 2m\omega _1 +2n\omega _3\) where \(m\) and \(n\) are integers (not both zero) is also a period. The set
of all such periods plus the origin form a lattice. In the literature \(-(\omega _1+\omega _3)\) is often denoted by \(\omega _2\) and
\(2\omega _2\) is clearly also a half-period; this accounts for the gap in the numbering of primitive
periods. The period \(\omega _2\) is little used in the REDUCE rule sets for the Weierstrassian elliptic
and related functions.
The primitive periods are not unique; indeed any periods \(2\Omega _1\) and \(2\Omega _3\) defined by the unimodular
integer bilinear transformation:
\[\Omega _1 = a\omega _1 + b\omega _3,\qquad \Omega _3 = c\omega _1 + d\omega _3,\qquad \text { where }ad-bc = 1\]
are also primitive. This fact is very useful in the
numerical evaluation of the Weierstrassian and Sigma functions as a sequence of such
transformations may be used to increase the size \(\Im \tau \) and so reduce the size of \(|q|\). Thus the
Fourier series for the theta functions and their derivatives will converge rapidly. In theory
these transformations may be used to reduce the size of \(|q|\) until \(\Im \tau \geq \sqrt 3/2\) when \(|q|<0.06\). However, in
numerical evaluations in REDUCE it is sufficient to use these transformations only until \(\Im \tau > 0.7\),
i.e. until \(|q| < 0.11\). In practice only two or three iterations are required and usually very much
smaller values of \(|q|\) are achieved particularly when \(\tau \) is purely imaginary i.e. \(q\) is
real.
The Weierstrassian function is even and has a pole of order 2 at all lattice points. The
Zeta and Sigma functions are only quasi-periodic on the lattice. Zeta is odd and has
simple poles of residue 1 at all lattice points. The basic Sigma function \(\sigma (u,\omega _1,\omega _3)\) is odd and
regular everywhere as is the function \(\vartheta _1(u,\tau )\) to which it is closely related. It has zeros at all
lattice points. All three functions \(\wp \), \(\zeta _w\) and \(\sigma \) are homogenous of degrees \(-2\), \(-1\) and \(+1\) respectively.
The functions are related by
\[ \wp (u) = -\zeta _w'(u),\qquad \qquad \zeta _w(u) = \sigma '(u)/\sigma (u),\]
where the lattice parameters have been omitted for
conciseness.
Rule sets are provided which implement all the properties such as double periodicity
discussed above. The rule for the derivative of the Weierstrass function has recently been
corrected; it involves the square root of \(4\wp (u)^3-g_2\wp (u)-g_3\) and to allow for the ambiguity of the sign of the
square root it includes an operator epsilon_w whose value is either \(+1\) or \(-1\). Its value
changes at the poles of the Weierstrass function, at points where the value of the
Weierstrass function is one of the three lattice roots \(e_1\), \(e_2\) or \(e_3\) and also where \(4\wp (u)^3-g_2\wp (u)-g_3\) becomes
negative as its square root has a branch cut there.
For numerical evaluation the switches rounded and complex must both be ON and
all three parameters of the Weierstrassian functions must be numeric. It is not, however,
necessary to ensure \(\Im (\omega _3/\omega _1) >0\) as the second and third parameters will be swapped if required.
Numerical evaluation of \(\epsilon _{w}\) has now been implemented for all complex values of its three
parameters, but currently it should be regarded as somewhat experimental; it uses fairly
crude numerical approximations to the derivative in order to choose the sign of thecorrect
branch of the square root. By and large the literature does not discuss how the sign may
be determined for general complex values of the arguments. One exception is the
Dover Handbook of Mathematical Functions[AS72] which gives a prescription
(without proof) only for the special cases where the lattice associated with the
Weierstrass function is either rectangular or rhombic, i.e. when the invariants \(g_2\) and \(g_3\) are
real.
Three commonly used alternative forms of the Weierstrassian functions in which
they are regarded as functions of the lattice invariants \(g_2\) and \(g_3\) rather than the primitive
half-periods \(\omega _1\) and \(\omega _3\) are provided:
- \(\wp (u \mid g_2, g_3)\) —
weierstrass1(u,g2,g3)
- \(\zeta _w(u \mid g_2, g_3)\) —
weierstrassZeta1(u,g2,g3)
- \(\sigma (u \mid g_2, g_3)\) —
weierstrass_sigma0(u,g2,g3).
Note the trailing zero rather than 1 in the third of these functions to avoid a clash with the
notation for the other sigma functions below. Note also that for output they are
distinguished by separating the first and second arguments by a vertical bar rather than a
comma.
The rule for differentiation of this variant of the Weierstrass function involves the
operator epsilon_w1 whose value is either \(+1\) or \(-1\) and is analogous to the operator
epsilon_w.
If the arguments of the Weierstrassian and sigma are all numeric, the functions will only
be evaluated if both rounded and complex are ON. The only exceptions are the
functions weierstrass1, weierstrassZeta1 and weierstrass_sigma0; these are real-valued
for real arguments and if all three arguments are real it suffices that the switch
rounded is ON. Similar remarks apply to those lattice functions (see below) which are
functions of the lattice invariants \(g_2\) and \(g_3\). Numerical evaluation of \(\epsilon _{w1}\) has now been
implemented for all complex values of its three parameters as for \(\epsilon _{w}\) although
for real values of its three parameters it suffices that the switch ROUNDED is
on.
Numerical Evaluation
If the discriminant \(g_2^3-27g_3^2\) vanishes, the functions Weierstrass1, WeierstrassZeta and
Weierstrass_sigma0 become elementary:
\begin{eqnarray*} \wp (u | 4b^4/3, 8b^6/27) &=& b^2(\mathrm {cosec}^2(b u)-1/3)\\ \wp (u | 4b^4/3, -8b^6/27) &=& b^2(\mathrm {cosech}^2(b u)+1/3)\\ \zeta _w(u |4b^4/3, 8b^6/27) &=& b(\cot (b u)+bu/3)\\ \zeta _w(u |4b^4/3, -8b^6/27) &=& b(\coth (b u)-bu/3)\\ \sigma (u |4b^4/3, 8b^6/27) &=& exp(b^2u^2/6)\sin (bu) \\ \sigma (u |4b^4/3, -8b^6/27) &=& exp(-b^2u^2/6)\sinh (bu) \end{eqnarray*}
where \(b\) is an arbitrary complex constant chosen for typographical simplicity. In these
cases the Weierstrass function becomes singly periodic, the other period becoming
infinite.
Laurent Series Expansions
Currently the Weierstrass and Weierstrass Zeta functions may be expanded as power
series using the extendible power series package tps. The expansions about
zero and the points \(m\omega _1+n\omega _3\) for integer \(m\) and \(n\) are informative whilst those about other
points whilst correct are not. Power series expansion of the Weierstrass function
about points where it vanishes (there are two such points in each fundamental
parallelogram) is also of interest. Currently the Weierstrass Sigma functions
(described below) cannot be expanded using tps and result in rather inelegant
errors.
Currently using the taylor package the expansions of the Weierstrass, Weierstrass Zeta
and Weierstrass Sigma about zero and \(2m\omega _1+2n\omega _3\) are incorrect whilst those about the mid-lattice
points \(m\omega _1+n\omega _3\) are formally correct, but not fully simplified.
Other Sigma Functions
Three further Sigma functions are also provided:
- \(\sigma _1(u, \omega _1, \omega _3)\) —
weierstrass_sigma1(u,omega1,omega3)
- \(\sigma _2(u, \omega _1, \omega _3)\) —
weierstrass_sigma2(u,omega1,omega3)
- \(\sigma _3(u, \omega _1, \omega _3)\) —
weierstrass_sigma3(u,omega1,omega3)
These are all even functions, regular everywhere, homogenous of degree zero and doubly
quasi-periodic. They are closely related to the theta functions \(\vartheta _2\), \(\vartheta _3\) and \(\vartheta _4\) respectively; but
note the difference in numbering. For more information on the properties these sigma
functions, see Lawden [Law89]; they do not appear in the NIST Digital Library of
Mathematical Functions, but are included here for completeness. These three functions
have no corresponding versions where the second and third arguments are the lattice
invariants \(g_2\) and \(g_3\).
Quasi-Period Factors & Lattice Functions
Ten functions are provided:
- \(e_1(\omega _1, \omega _3)\) —
lattice_e1(omega1, omega3);
- \(e_2(\omega _1, \omega _3)\) —
lattice_e2(omega1, omega3);
- \(e_3(\omega _1, \omega _3)\) —
lattice_e3(omega1, omega3);
- \(g_2(\omega _1, \omega _3)\) —
lattice_g2(omega1, omega3);
- \(g_3(\omega _1, \omega _3)\) —
lattice_g3(omega1, omega3);
- \(\Delta (\omega _1, \omega _3)\) —
lattice_delta(omega1, omega3);
- \(\mathrm {G}(\omega _1, \omega _3)\) —
lattice_g(omega1, omega3);
- \(\eta _1(\omega _1, \omega _3)\) —
eta1(omega1, omega3);
- \(\eta _2(\omega _1, \omega _3)\) —
eta2(omega1, omega3);
- \(\eta _3(\omega _1, \omega _3)\) —
eta3(omega1, omega3).
These are operative when the switches rounded and complex are ON and their
arguments are numerical. The first three are referred to as lattice roots and are related to
the invariants \(g_2, g_3\), the discriminant \(\Delta = g_2^3-27g_3^2\) and a closely related invariant \(\mathrm {G} = g_2^3/(27 g_3^2)\) of the Weierstrassian
elliptic function \(\wp \). The lattice roots also appear in the numerical evaluation of the
Weierstrass function. These lattice roots satisfy:
\[e_1+e_2+e_3=0,\qquad g_2=2(e_1^2+e_2^2+e_3^2),\qquad g_3= 4e_1e_2e_3.\]
If the discriminant \(\Delta \) vanishes or
equivalently if \(\mathrm {G} = 1\), there are at most two distinct lattice roots and the elliptic function
degenerates to an elementary one. The advantage of the invariant \(\mathrm {G}\) is that it is a function
of \(\tau = \omega _3/\omega _1\) only.
The remaining three functions eta1, eta2 & eta3 appear in the rules for the
quasi-periodicity of the four sigma functions and of the Weierstrassian Zeta function.
They are also used in the numerical evaluation of these functions when the
switches rounded and complex are ON. The quasi-period relations are:
\begin{align*} \zeta _w(u+2\omega _j) & = \zeta _w(u)+2\eta _j\\ \sigma (u+2\omega _j) & = -exp(2\eta _j(u+\omega _j))\sigma (u)\\ \sigma _k(u+2\omega _j) & = exp(2\eta _j(u+\omega _j))\sigma _k(u) \quad \text { if }j\neq k\\ \sigma _j(u+2\omega _j) & = -exp(2\eta _j(u+\omega _j))\sigma _j(u)\\ \zeta _w(\omega _j) & = \eta _j\\ \sigma _j(\omega _j) & = 0, \end{align*}
where the lattice parameters have been omitted for conciseness and \(j,k = 1\ldots 3\). The quasi-period
factors satisfy
\[\eta _1+\eta _2+\eta _3=0,\qquad \eta _1\omega _3-\eta _3\omega _1=\eta _2\omega _1-\eta _1\omega _2=\eta _3\omega _2-\eta _2\omega _3=i\pi /2.\]
Nine functions are provided which depend on the lattice invariants \(g_2\) and \(g_3\) rather than the
lattice half-periods:
- \(e_1(|g_2, g_3)\) —
lattice1_e1(g2, g3);
- \(e_2(|g_2, g_3)\) —
lattice1_e2(g2, g3);
- \(e_3(|g_2, g_3)\) —
lattice_e3(g2, g3);
- \(\omega _1(|g_2,g_3)\) —
lattice_omega1(g2, g3);
- \(\omega _2(|g_2,g_3)\) —
lattice_omega2(g2, g3);
- \(\omega _3(|g_2,g_3)\) —
lattice_omega3(g2, g3);
- \(\eta _1(|g_2,g_3)\) —
eta_1(g2,g3);
- \(\eta _2(|g_2,g_3)\) —
eta_2(g2,g3);
- \(\eta _3(|g_2,g_3)\) —
eta_3(g2,g3);
The first three return the lattice roots, the second triple return the three lattice
half-periods and the final three return the quasi-period factors. As with the alternative
Weierstrass functions which depend on the invariants, these nine functions are
distinguished from the similar functions depending on the periods by preceding their
arguments with a vertical bar.
As well as the scalar-valued functions discussed above in this section, there are six
functions which return a list as their value:
lattice_roots(omega1,omega3) — returns \(\{e_1,\ e_2,\ e_3\}\)
lattice_invariants(omega1,omega3) — returns \(\{g_2,\ g_3,\ \Delta ,\ \mathrm {G}\}\);
quasi_period_factors(omega1,omega3) — returns \(\{\eta _1,\ \eta _2,\ \eta _3\}\);
lattice_generators(g2,g3) — returns \(\{\omega _1,\ \omega _3\}\).
lattice1_roots(g2,g3) — returns \(\{e_1,\ e_2,\ e_3\}\)
quasi_periods(g2,g3) — returns \(\{\eta _1,\ \eta _2,\ \eta _3\}\).
The first three depend on the lattice half-periods \(\omega _1\) and \(\omega _3\) whilst the fourth, fifth
and sixth are functions of the invariants \(g_2\) and \(g_3\). These list-valued functions are
actually more efficient than calling the corresponding scalar-valued functions
individually. These functions are only useful when the switches rounded and
complex are ON and their arguments are all numerical. Note that the call
sequence:
lattice_generators(g2,g3);
lattice_invariants(first ws, second ws);
{first ws, second ws};
should reproduce the list {g2, g3}, perhaps with small rounding errors. The
corresponding sequence with the function lattice_generators being called after
lattice_invariants (and g2 & g3 replaced by w1 & w3), in general, will not
produce the same pair of lattice generators since the generators are only defined up to a
unimodular bilinear integer transformation.
For details of the algorithm used to calculate the lattice generators from the invariants see
the DLMF:NIST chapter on Lattice Calculations.
20.20.10 Inverse Jacobi Elliptic Functions
The following inverses of the 12 Jacobi elliptic functions are available:-
- arcsn(u,k)
- arcdn(u,k)
- arccn(u,k)
- arccd(u,k)
- arcsd(u,k)
- arcnd(u,k)
- arcdc(u,k)
- arcnc(u,k)
- arcsc(u,k)
- arcns(u,k)
- arcds(u,k)
- arccs(u,k)
Thus, for example,
jacobisn(arcsn(x, k), k) --> x
jacobisc(arcsc(x, k), k) --> x
A rule list is provided to simplify these functions for special values of their arguments
such \(x=0\), \(k=0\) and \(k=1\), to implement the inverse function simplification formulae illustrated
immediately above and for differentiation of these functions with respect to their two
arguments.
Note that \(\mathrm {arccs}\) is not defined to be an odd function of its first argument unlike \(\mathrm {cs}\). Instead
it is taken to satisfy:
\[ \mathrm {arccs}(-x, k) = 2\mathrm {K}(k)-\mathrm {arccs}(x, k).\]
This is analogous to the situation in Reduce for \(\mathrm {acot}\) where
\[ \mathrm {arctan}(-x) = -\mathrm {arctan}(x),\qquad \mathrm {arccot}(-x) = \pi -\mathrm {arccot}(x). \]
This choice means that the range of (real) principal values of \(\mathrm {arccs}\) is connected – it is the
open set \((0, 2K(k))\).
When their arguments are numerical, these functions will be evaluated numerically if the
rounded switch is ON. Note that in some cases the result may have an imaginary part
even if both arguments are real, hence the switch complex is always turned ON
temporarily during numerical evaluation and afterwards restored to its original
value.
Note also that for \(\mathrm {arcdn}\) and \(\mathrm {arcnd}\) a zero value of the modulus \(k\) is excluded (since \(\mathrm {dn}(x,0) = \mathrm {nd}(x,0) = 1 \quad \forall x\)).
As the Jacobi elliptic functions are doubly periodic, their inverse functions are
multi-valued. The numerical value returned is the principal value \(v\) which lies in the
parallelogram in the complex plane whose vertices are given in the table below. Other
values of the inverse functions are indicated in the fifth column of the table below where \(m\)
and \(n\) are arbitrary integers.
| Function | Quarter | Periods | Principal | Other |
| | \(p\) | \(q\) | Parallelogram | Values |
| | | | |
| \(\mathop {\mathrm {arcsn}}\) | \(\mathrm {K}\) | \(i\mathrm {K}'\) | \(-(p+q)\), \(-p+q\), | \(2m p+2n q +(-1)^mv\) |
| | | | \(p+q\), \(p-q\) | |
| \(\mathop {\mathrm {arccn}}\) | \(\mathrm {K}\) | \(\mathrm {K}+i\mathrm {K}'\) | \(-q\), \(q\), \(2p+q\), \(2p-q\) | \(4m p +2n q \pm v\) |
| \(\mathop {\mathrm {arcdn}}\) | \(i\mathrm {K}'\) | \(\mathrm {K}\) | \(0\), \(2p\), \(2(p+q)\), \(2q\) | \(2m q+4n p \pm v\) |
| \(\mathop {\mathrm {arcns}}\) | \(\mathrm {K}\) | \(i\mathrm {K}'\) | \(-(p+q)\), \(-p+q\), | \(2m p+2n q +(-1)^mv\) |
| | | | \(p+q\), \(p-q\) | |
| \(\mathop {\mathrm {arcnc}}\) | \(\mathrm {K}\) | \(\mathrm {K}+i\mathrm {K}'\) | \(-q\), \(q\), \(2p+q\), \(2p-q\) | \(4m p +2n q \pm v\) |
| \(\mathop {\mathrm {arcnd}}\) | \(i\mathrm {K}'\) | \(\mathrm {K}\) | \(0\), \(2p\), \(2(p+q)\), \(2q\) | \(2m q+4n p \pm v\) |
| \(\mathop {\mathrm {arccd}}\) | \(\mathrm {K}\) | \(i\mathrm {K}'\) | \(-q\), \(q\), \(2p+q\), \(2p-q\) | \(4m p +2n q \pm v\) |
| \(\mathop {\mathrm {arcdc}}\) | \(\mathrm {K}\) | \(i\mathrm {K}'\) | \(-q\), \(q\), \(2p+q\), \(2p-q\) | \(4m p +2n q \pm v\) |
| \(\mathop {\mathrm {arcsd}}\) | \(\mathrm {K}\) | \(\mathrm {K}+i\mathrm {K}'\) | \(-(p+q)\), \(-p+q\), | \(2m p +2n q +(-1)^mv\) |
| | | | \(p+q\), \(p-q\) | |
| \(\mathop {\mathrm {arcds}}\) | \(\mathrm {K}\) | \(\mathrm {K}+i\mathrm {K}'\) | \(-(p+q)\), \(-p+q\), | \(2m p +2n q +(-1)^mv\) |
| | | | \(p+q\), \(p-q\) | |
| \(\mathop {\mathrm {arcsc}}\) | \(i\mathrm {K}'\) | \(\mathrm {K}\) | \(-(p+q)\), \(-p+q\), | \(2m q +2n q +(-1)^nv\) |
| | | | \(p+q\), \(p-q\) | |
| \(\mathop {\mathrm {arccs}}\) | \(i\mathrm {K}'\) | \(\mathrm {K}\) | \(-p\), \(p\), \(p+2q\), \(-p+2q\) | \(2m q +2n q +(-1)^nv\) |
When both arguments are real and \(|k|<=1\) and when there are certain restrictions on the
range of the first parameter \(x\) (see the table below), then the principal value of
the inverse function is real. It lies in the range given in the third column of
the table below. (c.f. the inverse trigonometric functions). In these cases the
integrand in the defining integral has no branch points and a faster algorithm which
uses only real arithmetic is adopted. Other real values of the inverse functions
are indicated in the fourth column of the table below where \(n\) is an arbitrary
integer.
| Fn | Domain | Principal Value \(v\) | Other real values |
| \(\mathop {\mathrm {arcsn}}\): | \( |x| \le 1 \) | \(-\mathrm {K}(k) \le v \le \mathrm {K}(k)\) | \(2 n\mathrm {K}(k)+(-1)^nv\) |
| \(\mathop {\mathrm {arccn}}\): | \( |x| \le 1 \) | \(0 \le v \le 2\mathrm {K}(k)\) | \(4 n\mathrm {K}(k) \pm v\) |
| \(\mathop {\mathrm {arccd}}\): | \( |x| \le 1 \) | \(0 \le v \le 2\mathrm {K}(k)\) | \(4 n\mathrm {K}(k) \pm v\) |
| \(\mathop {\mathrm {arcns}}\): | \( |x| \ge 1 \) | \(-\mathrm {K}(k) \le v \le \mathrm {K}(k)\) & \(v \neq 0\) | \(2 n\mathrm {K}(k)+(-1)^nv\) |
| \(\mathop {\mathrm {arcnc}}\): | \( |x| \ge 1 \) | \(0 \le v \le 2\mathrm {K}(k)\) & \(v \neq \mathrm {K}(k)\) | \(4 n\mathrm {K}(k) \pm v\) |
| \(\mathop {\mathrm {arcdc}}\): | \( |x| \ge 1 \) | \(0 \le v \le 2\mathrm {K}(k)\) & \(v \neq \mathrm {K}(k)\) | \(4 n\mathrm {K}(k) \pm v\) |
| \(\mathop {\mathrm {arcdn}}\): | \( k' \le x \le 1\) | \(0 \le v \le \mathrm {K}(k)\) | \(2 n\mathrm {K}(k) \pm v\) |
| \(\mathop {\mathrm {arcnd}}\): | \( 1 \le x \le 1/k'\) | \(0 \le v \le \mathrm {K}(k)\) | \(2 n\mathrm {K}(k) \pm v\) |
| \(\mathop {\mathrm {arcds}}\): | \( |x| \ge k'\) | \(-\mathrm {K}(k) \le v \le \mathrm {K}(k)\) & \(v \neq 0\) | \(2 n\mathrm {K}(k)+(-1)^nv \) |
| \(\mathop {\mathrm {arcsd}}\): | \( |x| \le 1/k'\) | \(-\mathrm {K}(k) \le v \le \mathrm {K}(k)\) | \(2 n\mathrm {K}(k)+(-1)^nv \) |
| \(\mathop {\mathrm {arcsc}}\): | \( x \in \mathbb {R}\) | \(-\mathrm {K}(k) < v < \mathrm {K}(k)\) | \(2 n\mathrm {K}(k) + v\) |
| \(\mathop {\mathrm {arccs}}\): | \( x \in \mathbb {R}\) | \(0 < v < 2\mathrm {K}(k)\) | \(2 n\mathrm {K}(k) + v\ (v \neq 0)\) |
The numerical values of the inverse functions are calculated by expressing them in terms
of the symmetric elliptic integral:
\[ R_F(x,y,z)=\int _0^\infty 1/\sqrt {(t-x)(t-y)(t-z)}\,\mathrm {d}t. \]
For more details see the DLMF website: Inverse
Jacobian Elliptic Functions.
20.20.11 Table of Elliptic Functions and Integrals
| Function | Operator |
| | |
| \(\mathrm {am}(u,k)\) | jacobiam(u,k) |
| \(\mathrm {sn}(u,k)\) | jacobisn(u,k) |
| \(\mathrm {dn}(u,k)\) | jacobidn(u,k) |
| \(\mathrm {cn}(u,k)\) | jacobicn(u,k) |
| \(\mathrm {cd}(u,k)\) | jacobicd(u,k) |
| \(\mathrm {sd}(u,k)\) | jacobisd(u,k) |
| \(\mathrm {nd}(u,k)\) | jacobind(u,k) |
| \(\mathrm {dc}(u,k)\) | jacobidc(u,k) |
| \(\mathrm {nc}(u,k)\) | jacobinc(u,k) |
| \(\mathrm {sc}(u,k)\) | jacobisc(u,k) |
| \(\mathrm {ns}(u,k)\) | jacobins(u,k) |
| \(\mathrm {ds}(u,k)\) | jacobids(u,k) |
| \(\mathrm {cs}(u,k)\) | jacobics(u,k) |
| Inverse Functions of the above: |
| \(\mathrm {arcsn}(u,k)\) | arcsn(u,k) |
| \(\mathrm {arccn}(u,k)\) | arccn(u,k) |
 |
| \(\mathrm {arccs}(u,k)\) | arccs(u,k) |
| | |
| Complete Integral
|
| (1st kind) \(\mathrm {K}(k)\) | ellipticK(k) |
| \(\mathrm {K}'(k)\) | ellipticK!’(k) |
| Incomplete Integral
|
| (1st kind) \(\mathrm {F}(\phi ,k)\) | ellipticF(phi,k) |
| Complete Integral |
| (2nd kind) \(\mathrm {E}(k)\) | ellipticE(k) |
| \(\mathrm {E}'(k)\) | ellipticE!’(k) |
| Legendre Incomplete
|
| Integral (2nd kind) \(\mathrm {E}(u,k)\) | ellipticE(u,k) |
| Alternative Incomplete
|
| Integral (2nd kind) \(\mathrm {D}(u,k)\) | ellipticD(u,k) |
| Jacobi Incomplete
|
| Integral (2nd kind) \(\mathcal {E}(u,k)\) | jacobiE(u,k) |
| Jacobi’s Zeta \(\mathrm {Z}(u,k)\) | jacobiZeta(u,k) |
| \(\vartheta _1(u,\tau )\) | elliptictheta1(u,tau) |
| \(\vartheta _2(u,\tau )\) | elliptictheta2(u,tau) |
| \(\vartheta _3(u,\tau )\) | elliptictheta3(u,tau) |
| \(\vartheta _4(u,\tau )\) | elliptictheta4(u,tau) |
| \(\wp (u,\omega _1, \omega _3)\) | weierstrass(u,omega1,omega3) |
| \(\zeta _w(u,\omega _1, \omega _3)\) | weierstrassZeta(u,omega1,omega3) |
| \(\sigma (u,\omega _1, \omega _3)\) | weierstrass_sigma(u,omega1,omega3) |
| \(\sigma _1(u,\omega _1, \omega _3)\) | weierstrass_sigma1(u,omega1,omega3) |
| \(\sigma _2(u,\omega _1, \omega _3)\) | weierstrass_sigma2(u,omega1,omega3) |
| \(\sigma _3(u,\omega _1, \omega _3)\) | weierstrass_sigma3(u,omega1,omega3) |
| \(\wp (u \mid g_2, g_3)\) | weierstrass1(u,g2,g3) |
| \(\zeta _w(u \mid g_2, g_3)\) | weierstrassZeta1(u,g2,g3) |
| \(\sigma (u \mid g_2, g_3)\) | weierstrass_sigma0(u,g2,g3) |
| | |
Hosted by 
