§22.20 Methods of Computation

A powerful way of computing the twelve Jacobian elliptic functions for real or complex values of both the argument $z$ and the modulus $k$ is to use the definitions in terms of theta functions given in §22.2, obtaining the theta functions via methods described in §20.14.

Given real or complex numbers $a_0,b_0$, with $b_0/a_0$ not real and negative, define

for $n\ge 1$, where the square root is chosen so that $\mathrm{ph}{b}_n=\frac{1}{2}(\mathrm{ph}{a}_{n-1}+\mathrm{ph}{b}_{n-1})$, where $\mathrm{ph}{a}_{n-1}$ and $\mathrm{ph}{b}_{n-1}$ are chosen so that their difference is numerically less than $\pi$. Then as $n\to \mathrm{\infty }$ sequences $\{a_n\}$, $\{b_n\}$ converge to a common limit $M=M(a_0,b_0)$, the arithmetic-geometric mean of $a_0,b_0$. And since

convergence is very rapid.

For $x$ real and $k\in (0,1)$, use (22.20.1) with $a_0=1$, $b_0=k^{\prime }\in (0,1)$, $c_0=k$, and continue until $c_N$ is zero to the required accuracy. Next, compute ${\varphi }_N,{\varphi }_{N-1},\mathrm{\dots },{\varphi }_0$, where

and the inverse sine has its principal value (§4.23(ii)). Then

and the subsidiary functions can be found using (22.2.10). This formula for $\mathrm{dn}$ becomes unstable near $x=K$. If only the value of $\mathrm{dn}(x,k)$ at $x=K$ is required then the exact value is in the table 22.5.1. If both $k$ and $x$ are real then $\mathrm{dn}$ is strictly positive and $\mathrm{dn}(x,k)=\sqrt{1-k^2{\mathrm{sn}}^2(x,k)}$ which follows from (22.6.1). If either $k$ or $x$ is complex then (22.2.6) gives the definition of $\mathrm{dn}(x,k)$ as a quotient of theta functions.

See also Wachspress (2000).

By application of the transformations given in §§22.7(i) and 22.7(ii), $k$ or $k^{\prime }$ can always be made sufficently small to enable the approximations given in §22.10(ii) to be applied. The rate of convergence is similar to that for the arithmetic-geometric mean.

If either $\tau$ or $q={\mathrm{e}}^{\mathrm{i}\pi \tau }$ is given, then we use $k={\theta }_2^2(0,q)/{\theta }_3^2(0,q)$, $k^{\prime }={\theta }_4^2(0,q)/{\theta }_3^2(0,q)$, $K=\frac{1}{2}\pi {\theta }_3^2(0,q)$, and $K^{\prime }=-\mathrm{i}\tau K$, obtaining the values of the theta functions as in §20.14.

If $k,k^{\prime }$ are given with $k^2+k_{}^{\prime }{}_{}{}^2=1$ and , then $K,K^{\prime }$ can be found from

using the arithmetic-geometric mean.

See Wachspress (2000).

$\mathrm{am}(x,k)$ can be computed from its definition (22.16.1) or from its Fourier series (22.16.9). Alternatively, Sala (1989) shows how to apply the arithmetic-geometric mean to compute $\mathrm{am}(x,k)$.

Jacobi’s epsilon function can be computed from its representation (22.16.30) in terms of theta functions and complete elliptic integrals; compare §20.14. Jacobi’s zeta function can then be found by use of (22.16.32).

For additional information on methods of computation for the Jacobi and related functions, see the introductory sections in the following books: Lawden (1989), Curtis (1964b), Milne-Thomson (1950), and Spenceley and Spenceley (1947).