§22.20 Methods of Computation

§22.20(i) Via Theta Functions

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.

§22.20(ii) Arithmetic-Geometric Mean

Given real or complex numbers $a_{0},b_{0}$, with $b_{0}/a_{0}$ not real and negative, define

 22.20.1 $\displaystyle a_{n}$ $\displaystyle=\tfrac{1}{2}\left(a_{n-1}+b_{n-1}\right),$ $\displaystyle b_{n}$ $\displaystyle=\left(a_{n-1}b_{n-1}\right)^{1/2},$ $\displaystyle c_{n}$ $\displaystyle=\tfrac{1}{2}\left(a_{n-1}-b_{n-1}\right),$ ⓘ Symbols: $a_{n}$: numbers, $b_{n}$: numbers, $c_{n}$: numbers and $n$: positive Referenced by: §22.20(ii), §22.20(iv), §22.20(iv), §22.20(ii) Permalink: http://dlmf.nist.gov/22.20.E1 Encodings: TeX, TeX, TeX, pMML, pMML, pMML, png, png, png See also: Annotations for §22.20(ii), §22.20 and Ch.22

for $n\geq 1$, where the square root is chosen so that $\operatorname{ph}b_{n}=\tfrac{1}{2}(\operatorname{ph}a_{n-1}+\operatorname{ph}% b_{n-1})$, where $\operatorname{ph}a_{n-1}$ and $\operatorname{ph}b_{n-1}$ are chosen so that their difference is numerically less than $\pi$. Then as $n\to\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

 22.20.2 $\max\left(\left|a_{n}-M\right|,\left|b_{n}-M\right|,\left|c_{n}\right|\right)% \leq\text{(const.)}\times 2^{-2^{n}},$ ⓘ Symbols: $a_{n}$: numbers, $b_{n}$: numbers, $c_{n}$: numbers, $n$: positive and $M(a_{0},b_{0})$: limit Permalink: http://dlmf.nist.gov/22.20.E2 Encodings: TeX, pMML, png See also: Annotations for §22.20(ii), §22.20 and Ch.22

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 $\phi_{N},\phi_{N-1},\dots,\phi_{0}$, where

 22.20.3 $\phi_{N}=2^{N}a_{N}x,$ ⓘ Symbols: $x$: real, $a_{n}$: numbers and $\phi_{N}$: approximations Referenced by: §22.20(ii) Permalink: http://dlmf.nist.gov/22.20.E3 Encodings: TeX, pMML, png See also: Annotations for §22.20(ii), §22.20 and Ch.22
 22.20.4 $\phi_{n-1}=\frac{1}{2}\left(\phi_{n}+\operatorname{arcsin}\left(\frac{c_{n}}{a% _{n}}\sin\phi_{n}\right)\right),$

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

 22.20.5 $\displaystyle\operatorname{sn}\left(x,k\right)$ $\displaystyle=\sin\phi_{0},$ $\displaystyle\operatorname{cn}\left(x,k\right)$ $\displaystyle=\cos\phi_{0},$ $\displaystyle\operatorname{dn}\left(x,k\right)$ $\displaystyle=\frac{\cos\phi_{0}}{\cos\left(\phi_{1}-\phi_{0}\right)},$ ⓘ Symbols: $\operatorname{cn}\left(\NVar{z},\NVar{k}\right)$: Jacobian elliptic function, $\operatorname{dn}\left(\NVar{z},\NVar{k}\right)$: Jacobian elliptic function, $\operatorname{sn}\left(\NVar{z},\NVar{k}\right)$: Jacobian elliptic function, $K\left(\NVar{k}\right)$: Legendre’s complete elliptic integral of the first kind, $\cos\NVar{z}$: cosine function, $\sin\NVar{z}$: sine function, $x$: real, $k$: modulus and $\phi_{N}$: approximations Referenced by: (22.20.5), §22.20(ii), 3rd item Permalink: http://dlmf.nist.gov/22.20.E5 Encodings: TeX, TeX, TeX, pMML, pMML, pMML, png, png, png Note (effective with 1.0.10): A note was added after (22.20.5) to deal with cases when computation of $\operatorname{dn}\left(x,k\right)$ with (22.20.5) becomes numerically unstable near $x=K$. Suggested 2014-10-20 by Hartmut Henkel See also: Annotations for §22.20(ii), §22.20 and Ch.22

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

Example

To compute $\operatorname{sn}$, $\operatorname{cn}$, $\operatorname{dn}$ to 10D when $x=0.8$, $k=0.65$.

Four iterations of (22.20.1) lead to $c_{4}=\Sci{6.5}{-12}$. From (22.20.3) and (22.20.4) we obtain $\phi_{1}=1.40213\;91827$ and $\phi_{0}=0.76850\;92170$. Then from (22.20.5), $\operatorname{sn}\left(0.8,0.65\right)=0.69506\;42165$, $\operatorname{cn}\left(0.8,0.65\right)=0.71894\;76580$, $\operatorname{dn}\left(0.8,0.65\right)=0.89212\;34349$.

§22.20(iii) Landen Transformations

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.

Example

To compute $\operatorname{dn}\left(x,k\right)$ to 6D for $x=0.2$, $k^{2}=0.19$, $k^{\prime}=0.9$.

From (22.7.1), $k_{1}=\tfrac{1}{19}$ and $x/(1+k_{1})=0.19$. From the first two terms in (22.10.6) we find $\operatorname{dn}\left(0.19,\tfrac{1}{19}\right)=0.999951$. Then by using (22.7.4) we have $\operatorname{dn}\left(0.2,\sqrt{0.19}\right)=0.996253$.

If needed, the corresponding values of $\operatorname{sn}$ and $\operatorname{cn}$ can be found subsequently by applying (22.10.4) and (22.7.2), followed by (22.10.5) and (22.7.3).

§22.20(iv) Lattice Calculations

If either $\tau$ or $q=e^{i\pi\tau}$ is given, then we use $k={\theta_{2}^{2}}\left(0,q\right)/{\theta_{3}^{2}}\left(0,q\right)$, $k^{\prime}={\theta_{4}^{2}}\left(0,q\right)/{\theta_{3}^{2}}\left(0,q\right)$, $K=\frac{1}{2}\pi{\theta_{3}^{2}}\left(0,q\right)$, and $K^{\prime}=-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 $\Im k^{\prime}/\Im k<0$, then $K,K^{\prime}$ can be found from

 22.20.6 $\displaystyle K$ $\displaystyle=\frac{\pi}{2M(1,k^{\prime})},$ $\displaystyle K^{\prime}$ $\displaystyle=\frac{\pi}{2M(1,k)},$

using the arithmetic-geometric mean.

Example 1

If $k=k^{\prime}=1/\sqrt{2}$, then three iterations of (22.20.1) give $M=0.84721\;30848$, and from (22.20.6) $K=\pi/(2M)=1.85407\;46773$ — in agreement with the value of $\left(\Gamma\left(\tfrac{1}{4}\right)\right)^{2}/\left(4\sqrt{\pi}\right)$; compare (23.17.3) and (23.22.2).

Example 2

If $k^{\prime}=1-i$, then four iterations of (22.20.1) give $K=1.23969\;74481+i0.56499\;30988$.

§22.20(v) Inverse Functions

See Wachspress (2000).

§22.20(vi) Related Functions

$\operatorname{am}\left(x,k\right)$ 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 $\operatorname{am}\left(x,k\right)$.

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).

§22.20(vii) Further References

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).