Elliptic functions for Matlab and Octave

TheMatlab script implementations ofElliptic integrals of three types,Jacobi's elliptic functions andJacobi theta functions of four types.

The mainGOAL of the project is to provide the natural Matlab scriptsWITHOUT external library calls like Maple and others. All scripts are developed to accept tensors as arguments and almost all of them have their complex versions. Performance and complete control on the execution are the main features.

Citations and references

If you've used any of the routines in this package please cite and support the effort. Here is the example of the BibTeX entry

@misc{elliptic,  author = {Moiseev I.},  title = {Elliptic functions for Matlab and Octave},  year = {2008},  publisher = {GitHub},  journal = {GitHub repository},  howpublished = {\url{https://github.com/moiseevigor/elliptic}},  commit = {98181c4c0d8992746bcc6bea75740bb11b74b51b},  doi = {10.5281/zenodo.48264},  url = {http://dx.doi.org/10.5281/zenodo.48264}}

or simply

Moiseev I., Elliptic functions for Matlab and Octave, (2008), GitHub repository, DOI: http://dx.doi.org/10.5281/zenodo.48264

Contents of the package

• Elliptic Functions
  • ELLIPJ: Jacobi's elliptic functions
  • ELLIPJI: Jacobi's elliptic functions of the complex argument
  • JACOBITHETAETA: Jacobi's Theta and Eta Functions
  • THETA: Theta Functions of Four Types

• Elliptic Integrals

  • ELLIPTIC12: Incomplete Elliptic Integrals of the First, Second Kind and Jacobi's Zeta Function
  • ELLIPTIC12I: Incomplete Elliptic Integrals of the First, Second Kind and Jacobi's Zeta Function of the complex argument
  • ELLIPTIC3: Incomplete Elliptic Integral of the Third Kind
  • ELLIPTIC123: Complete and Incomplete Elliptic Integrals of the First, Second, and Third Kind
  • INVERSELLIPTIC2: INVERSE Incomplete Elliptic Integrals of the Second Kind

• Weierstrass's elliptic functions (in development)

• Elliptic Related Functions
  • AGM: Artihmetic Geometric Mean
  • NOMEQ: The Value of Nome q = q(m)
  • INVERSENOMEQ: The Value of Nome m = m(q)
• Contributors
• References

Elliptic Functions

The Jacobi's elliptic functions are a set of basic elliptic functions, and auxiliary theta functions, that have historical importance with also many features that show up important structure, and have direct relevance to some applications (e.g. the equation of the pendulum). They also have useful analogies to the functions of trigonometry, as indicated by the matching notationSN forSIN. They are not the simplest way to develop a general theory, as now seen: that can be said for the Weierstrass elliptic functions. They are not, however, outmoded. They were introduced by Carl Gustav Jakob Jacobi, around 1830.

Theta functions are special functions of several complex variables. They are important in several areas, including the theories of abelian varieties and moduli spaces, and of quadratic forms. They have also been applied to soliton theory. When generalized to a Grassmann algebra, they also appear in quantum field theory, specifically string theory and D-branes.

ELLIPJ: Jacobi's elliptic functions

ELLIPJ evaluates theJacobi's elliptic functions andJacobi's amplitude.

[Sn,Cn,Dn,Am] = ELLIPJ(U,M) returns the values of the Jacobi elliptic functionsSN,CN,DN andAM evaluated for corresponding elements of argument U and parameter M. The arrays U and M must be of the same size (or either can be scalar). As currently implemented, M is limited to0 <= M <= 1.

General definition:

u = Integral(1/sqrt(1-m^2*sin(theta)^2), 0, phi);Sn(u) = sin(phi);Cn(u) = cos(phi);Dn(u) = sqrt(1-m^2*sin(phi)^2);

Depends onAGM,ELLIPKE.
Used byTHETA.
See alsoELLIPKE.

ELLIPJI: Jacobi's elliptic functions of the complex argument

ELLIPJI evaluates the Jacobi elliptic functions of complex phaseU.

[Sni,Cni,Dni] = ELLIPJ(U,M) returns the values of the Jacobi elliptic functionsSNI,CNI andDNI evaluated for corresponding elements of argumentU and parameterM. The arraysU andM must be of the same size (or either can be scalar). As currently implemented,M is real and limited to0 <= M <= 1.

Example:

[phi1,phi2] = meshgrid(-pi:3/20:pi, -pi:3/20:pi);phi = phi1 + phi2*i;[Sni,Cni,Dni]= ellipji(phi, 0.99);

Depends onAGM,ELLIPJ,ELLIPKE
See alsoELLIPTIC12,ELLIPTIC12I

JACOBITHETAETA: Jacobi's Theta and Eta Functions

JACOBITHETAETA evaluates Jacobi's theta and eta functions.

[Th, H] = JACOBITHETAETA(U,M) returns the values of the Jacobi's theta and eta elliptic functionsTH andH evaluated for corresponding elements of argumentU and parameterM. The arraysU andM must be the same size (or either can be scalar). As currently implemented,M is real and limited to0 <= M <= 1.

Example:

[phi,alpha]= meshgrid(0:5:90, 0:2:90);                  [Th, H] = jacobiThetaEta(pi/180*phi, sin(pi/180*alpha).^2);

Depends onAGM,ELLIPJ,ELLIPKE
See alsoELLIPTIC12,ELLIPTIC12I,THETA

THETA: Theta Functions of Four Types

THETA evaluates theta functions of four types.

Th = THETA(TYPE,V,M) returns values of theta functionsevaluated for corresponding values of argumentV and parameterM.TYPE is a type of the theta function, there are four numbered types. The arraysV andM must be the same size (or either can be scalar). As currently implemented,M is limited to0 <= M <= 1.

Th = THETA(TYPE,V,M,TOL) computes the theta and eta elliptic functions to the accuracyTOL instead of the defaultTOL = EPS.

The parameterM is related to the nomeQ asQ = exp(-pi*K(1-M)/K(M)). Some definitions of the Jacobi's elliptic functions use the modulusk instead of the parameterm. They are related bym = k^2.

Example:

[phi,alpha] = meshgrid(0:5:90, 0:2:90);                  Th1 = theta(1, pi/180*phi, sin(pi/180*alpha).^2);  Th2 = theta(2, pi/180*phi, sin(pi/180*alpha).^2);  Th3 = theta(3, pi/180*phi, sin(pi/180*alpha).^2);  Th4 = theta(4, pi/180*phi, sin(pi/180*alpha).^2);

Depends onAGM,ELLIPJ,ELLIPKE,JACOBITHETAETA
See alsoELLIPTIC12,ELLIPTIC12I

Elliptic Integrals

Elliptic integrals originally arose in connection with the problem of giving the arc length of an ellipse. They were first studied by Giulio Fagnano and Leonhard Euler. Modern mathematics defines an elliptic integral as any function f which can be expressed in the form

f(x) = Integral(R(t,P(t), c, x)dt,

whereR is a rational function of its two arguments,P is the square root of a polynomial of degree3 or4 with no repeated roots, andc is a constant.In general, elliptic integrals cannot be expressed in terms of elementary functions. Exceptions to this general rule are whenP has repeated roots, or whenR(x,y) contains no odd powers ofy. However, with the appropriate reduction formula, every elliptic integral can be brought into a form that involves integrals over rational functions and the three canonical forms (i.e. the elliptic integrals of the first, second and third kind).

ELLIPTIC12: Incomplete Elliptic Integrals of the First, Second Kind and Jacobi's Zeta Function

ELLIPTIC12 evaluates the value of the Incomplete Elliptic Integrals of the First, Second Kind and Jacobi's Zeta Function.

[F,E,Z] = ELLIPTIC12(U,M,TOL) uses the method of the Arithmetic-Geometric Mean and Descending Landen Transformation described in1 Ch. 17.6, to determine the value of the Incomplete Elliptic Integrals of the First, Second Kind and Jacobi's Zeta Function (see1,2).

General definition:

F(phi,m) = int(1/sqrt(1-m*sin(t)^2), t=0..phi);E(phi,m) = int(sqrt(1-m*sin(t)^2), t=0..phi);Z(phi,m) = E(u,m) - E(m)/K(m)*F(phi,m).

Tables generating code (see1, pp. 613-621):

[phi,alpha] = meshgrid(0:5:90, 0:2:90);                  % modulus and phase in degrees[F,E,Z] = elliptic12(pi/180*phi, sin(pi/180*alpha).^2);  % values of integrals

Depends onAGM
See alsoELLIPKE,ELLIPJ,ELLIPTIC12I,ELLIPTIC3,THETA.

ELLIPTIC12I: Incomplete Elliptic Integrals of the First, Second Kind and Jacobi's Zeta Function of the complex argument

ELLIPTIC12i evaluates the Incomplete Elliptic Integrals of the First, Second Kind and Jacobi's Zeta Function for the complex value of phaseU. ParameterM must be in the range0 <= M <= 1.

[Fi,Ei,Zi] = ELLIPTIC12i(U,M,TOL) whereU is a complex phase in radians,M is the real parameter andTOL is the tolerance (optional). Default value for the tolerance iseps = 2.220e-16.

ELLIPTIC12i uses the functionELLIPTIC12 to evaluate the values of corresponding integrals.

Example:

[phi1,phi2] = meshgrid(-2*pi:3/20:2*pi, -2*pi:3/20:2*pi);phi = phi1 + phi2*i;[Fi,Ei,Zi] = elliptic12i(phi, 0.5);

Depends onELLIPTIC12,AGM
See alsoELLIPKE,ELLIPJ,ELLIPTIC3,THETA.

ELLIPTIC3: Incomplete Elliptic Integral of the Third Kind

ELLIPTIC3 evaluates incomplete elliptic integral of the third kindPi = ELLIPTIC3(U,M,C) whereU is a phase in radians,0 < M < 1 is the module and0 < C < 1 is a parameter.

ELLIPTIC3 uses Gauss-Legendre 10 points quadrature template described in [3] to determine the value of the Incomplete Elliptic Integral of the Third Kind (see [1, 2]).

General definition:

Pi(u,m,c) = int(1/((1 - c*sin(t)^2)*sqrt(1 - m*sin(t)^2)), t=0..u)

Tables generating code (1, pp. 625-626):

[phi,alpha,c] = meshgrid(0:15:90, 0:15:90, 0:0.1:1);Pi = elliptic3(pi/180*phi, sin(pi/180*alpha).^2, c);  % values of integrals

ELLIPTIC123: Complete and Incomplete Elliptic Integrals of the First, Second, and Third Kind

ELLIPTIC123 is a wrapper around the different elliptic integral functions, providing a unified interface and greater range of input parameters. (Unlike ELLIPKE, ELLIPTIC12 and ELLIPTIC3, which all require a phase between zero and pi/2 and a parameter between zero and one.)

[F,E] = ELLIPTIC123(m) — complete Elliptic Integrals of the first and second kind.

[F,E] = ELLIPTIC123(b,m) — incomplete Elliptic Integrals of the first and second kind.

[F,E,PI] = ELLIPTIC123(m,n) — complete Elliptic Integrals of the first to third kind.

[F,E,PI] = ELLIPTIC123(b,m,n) — incomplete Elliptic Integrals of the first to third kind.

The order of the input arguments has been chosen to be consistent with the pre-existingelliptic12 andelliptic3 functions.

This function is still under development and its results are not always well-defined or even able to be calculated (especially for the third elliptic integral withn>1). Please see the documentation for further details.

INVERSELLIPTIC2: INVERSE Incomplete Elliptic Integrals of the Second Kind

INVERSELLIPTIC2 evaluates the value of the INVERSE Incomplete Elliptic Integrals of the Second Kind.

INVERSELLIPTIC2 uses the method described by Boyd J. P. to determine the value of the inverse Incomplete Elliptic Integrals of the Second Kind using the “Empirical” initialization to the Newton’s iteration method7.

Elliptic integral of the second kind:

E(phi,m) = int(sqrt(1-m*sin(t)^2), t=0..phi);

“Empirical” initialization7:

T0(z,m) = pi/2 + sqrt(r)/(theta − pi/2)

where

z \in [E(pi/2,m)](−E(pi/2,m),)x[1](0,) - value of the entire parameter spacer = sqrt((1-m)^2 + zeta^2)zeta = 1 - z/E(pi/2,m)theta = atan((1 - m)/zeta)

Example:

% modulus and phase in degrees[phi,alpha] = meshgrid(0:5:90, 0:2:90);% values of integrals[F,E] = elliptic12(pi/180*phi, sin(pi/180*alpha).^2);% values of inverse invE = inverselliptic2(E, sin(pi/180*alpha).^2);% the difference between phase phi and invE should close to zerophi - invE * 180/pi

Weierstrass's elliptic functions (in development)

!IN DEVELOPMENT, help needed!.

Weierstrass's elliptic functions are elliptic functions that take a particularly simple form (cf Jacobi's elliptic functions); they are named for Karl Weierstrass. This class of functions are also referred to as p-functions and generally written using the symbol℘ (a stylised letter p called Weierstrass p).

The Weierstrass elliptic function can be defined in three closely related ways, each of which possesses certain advantages. One is as a function of a complex variable z and a lattice Λ in the complex plane. Another is in terms of z and two complex numbers ω1 and ω2 defining a pair of generators, or periods, for the lattice. The third is in terms z and of a modulus τ in the upper half-plane. This is related to the previous definition byτ = ω2 / ω1, which by the conventional choice on the pair of periods is in the upper half-plane. Using this approach, for fixed z the Weierstrass functions become modular functions of τ.

Elliptic Related Functions

AGM: Artihmetic Geometric Mean

AGM calculates theArtihmetic Geometric Mean ofA andB (see1).

[A,B,C,N]= AGM(A0,B0,C0,TOL) carry out the process of the arithmetic geometric mean, starting with a given positive numbers triple(A0, B0, C0) and returns in(A, B, C) the generated sequence.N is a number of steps (returns in the valueuint32).

The general scheme of the procedure:

A(i) = 1/2*( A(i-1)+B(i-1) );     A(0) = A0;B(i) = sqrt( A(i-1)*B(i-1) );     B(0) = B0;C(i) = 1/2*( A(i-1)+B(i-1) );     C(0) = C0;

Stop at theN-th step whenA(N) = B(N), i.e., whenC(N) = 0.

Used byELLIPJ andELLIPTIC12.
See alsoELLIPKE,ELLIPTIC3,THETA.

NOMEQ: The Value of Nomeq = q(m)

NOMEQ gives the value of Nomeq = q(m).

NomeQ = nomeq(M,TOL), where0<=M<=1 is the module andTOL is the tolerance (optional). Default value for the tolerance iseps = 2.220e-16.

Used byELLIPJ.
Depends onELLIPKE
See alsoELLIPTIC12I,ELLIPTIC3,THETA.

INVERSENOMEQ: The Value of Nomem = m(q)

INVERSENOMEQ gives the value of Nomem = m(q).

M = inversenomeq(q), whereQ is the Nome of q-series.

WARNING. The functionINVERSENOMEQ does not return correct values ofM forQ > 0.6, because of computer precision limitation. The functionNomeQ(m) has an essential singularity atM = 1, so it cannot be inverted at this point and actually it is very hard to find and inverse in the neigborhood also.

More preciesly:

nomeq(1) = 1nomeq(1-eps) = 0.77548641878026

Example:

nomeq(inversenomeq([0.3 0.4 0.5 0.6 0.7 0.8](0.001)))

Used byELLIPJ.
Depends onELLIPKE
See alsoELLIPTIC12I,ELLIPTIC3,THETA.

Contributors

Contributors

• @moiseevigor - maintainer (Igor Moiseev)
• @drbitboy (Brian Carcich)
• @wspr (Will Robertson)

References

1. NIST Digital Library of Mathematical Functions
2. M. Abramowitz and I.A. Stegun, "Handbook of Mathematical Functions" Dover Publications", 1965, Ch. 17.1 - 17.6.
3. D. F. Lawden, "Elliptic Functions and Applications" Springer-Verlag, vol. 80, 1989
4. S. Zhang, J. Jin, "Computation of Special Functions" (Wiley, 1996).
5. B. Carlson, "Three Improvements in Reduction and Computation of Elliptic Integrals", J. Res. Natl. Inst. Stand. Technol. 107 (2002) 413-418.
6. N. H. Abel, "Studies on Elliptic Functions", english translation from french by Marcus Emmanuel Barnes. Original "Recherches sur les fonctions elliptiques", Journal fr die reine und angewandte Mathematik, Vol. 2, 1827. pp. 101-181.
7. B. C. Berndt, H. H. Chan, S.-S. Huang, "Incomplete Elliptic Integrals in Ramanujan's Lost Notebook" in q-series from a Contemporary Perspective, M. E. H. Ismail and D. Stanton, eds., Amer. Math. Soc., 2000, pp. 79-126.
8. J. P. Boyd, "Numerical, Perturbative and Chebyshev Inversion of the Incomplete Elliptic Integral of the Second Kind", Applied Mathematics and Computation (January 2012)