About the Project
22 Jacobian Elliptic FunctionsComputation

§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 a0,b0, with b0/a0 not real and negative, define

22.20.1 an =12(an-1+bn-1),
bn =(an-1bn-1)1/2,
cn =12(an-1-bn-1),

for n1, where the square root is chosen so that phbn=12(phan-1+phbn-1), where phan-1 and phbn-1 are chosen so that their difference is numerically less than π. Then as n sequences {an}, {bn} converge to a common limit M=M(a0,b0), the arithmetic-geometric mean of a0,b0. And since

22.20.2 max(|an-M|,|bn-M|,|cn|)(const.)×2-2n,

convergence is very rapid.

For x real and k(0,1), use (22.20.1) with a0=1, b0=k(0,1), c0=k, and continue until cN is zero to the required accuracy. Next, compute ϕN,ϕN-1,,ϕ0, where

22.20.3 ϕN=2NaNx,
22.20.4 ϕn-1=12(ϕn+arcsin(cnansinϕn)),

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

22.20.5 sn(x,k) =sinϕ0,
cn(x,k) =cosϕ0,
dn(x,k) =cosϕ0cos(ϕ1-ϕ0),

and the subsidiary functions can be found using (22.2.10). This formula for dn becomes unstable near x=K. If only the value of 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 dn is strictly positive and dn(x,k)=1-k2sn2(x,k) which follows from (22.6.1). If either k or x is complex then (22.2.6) gives the definition of dn(x,k) as a quotient of theta functions.

See also Wachspress (2000).


To compute sn, cn, dn to 10D when x=0.8, k=0.65.

Four iterations of (22.20.1) lead to c4=6.5×10⁻¹². From (22.20.3) and (22.20.4) we obtain ϕ1=1.40213 91827 and ϕ0=0.76850 92170. Then from (22.20.5), sn(0.8,0.65)=0.69506 42165, cn(0.8,0.65)=0.71894 76580, dn(0.8,0.65)=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 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.


To compute dn(x,k) to 6D for x=0.2, k2=0.19, k=0.9.

From (22.7.1), k1=119 and x/(1+k1)=0.19. From the first two terms in (22.10.6) we find dn(0.19,119)=0.999951. Then by using (22.7.4) we have dn(0.2,0.19)=0.996253.

If needed, the corresponding values of sn and 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 τ or q=eiπτ is given, then we use k=θ22(0,q)/θ32(0,q), k=θ42(0,q)/θ32(0,q), K=12πθ32(0,q), and K=-iτK, obtaining the values of the theta functions as in §20.14.

If k,k are given with k2+k2=1 and k/k<0, then K,K can be found from

22.20.6 K =π2M(1,k),
K =π2M(1,k),

using the arithmetic-geometric mean.

Example 1

If k=k=1/2, then three iterations of (22.20.1) give M=0.84721 30848, and from (22.20.6) K=π/(2M)=1.85407 46773 — in agreement with the value of (Γ(14))2/(4π); compare (23.17.3) and (23.22.2).

Example 2

If k=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

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

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