Elfun18 A collection of Matlab functions for the computation of Elliptical Integrals and Jacobian elliptic functions of real arguments

In the article we outline the set of Matlab functions that enable the computation of elliptic Integrals and Jacobian elliptic functions for real arguments. Correctness, robustness, efficiency and accuracy of the functions are discussed in some details. An example from the elasticity theory illustrates use of the collection.


Introduction
Elfun18 is a collection of the Matlab functions that enable the computation of elliptic Integrals and Jacobian elliptic functions for real arguments. Its purpose is to provide a sufficiently comprehensive set of elliptic integrals and elliptic functions which one can meet in practice. Altogether the collection contains 75 different functions (see Appendix). In particular, the collection covers all elliptic integrals and Jacobian elliptic functions given in [1][2][3][4][5]. Most of these functions are for example included in computer algebra programs Maple and Mathematica and computational environment MATLAB. Also, several software libraries include elliptic integrals and elliptic functions (AMath/DAMath0 F 1 , BOOST1 F 2 , CERN2 F 3 , elliptic3 F 4 , DATAPLOT4 F 5 , SLATEC5 F 6 , Mathematics Source Library6 F 7 , MPMATH7 F 8 , MATHCW8 F 9 ). However, each of these programs and libraries has its limitations; either in a set of available functions or in the choice of function arguments. The reason for that is probably a maze of possibilities by which one can define elliptic integrals and functions. Thus the elliptic integrals can be given in three different forms [6]. Also, an argument of every elliptic integral or function is either the modular angle, the modulus or the parameter [3,4]. To that, we add that for the Legendre form of elliptic integrals some these programs and libraries ignore the fact that these integrals are quasi-periodic functions and therefore give incorrect results.
In the next section, we will outline the elliptic integrals and functions. In section 3 we will describe the collection structure and its functions. Section 4 we will discuss the collection correctness, robustness, efficiency and accuracy. Section 5 provides two examples from elasticity theory.

On elliptic integrals and elliptic functions
General considerations. As already mentioned an argument of any elliptic integral or any elliptic function is either the modular angle α, the modulus sin k α = , or the parameter 2 m k = . Therefore three different notations are used [3] ( ) where Q is an elliptic integral or elliptic function. We shall call these a-functions, k-functions, and mfunctions. Because our goal is a computation, we will as an elliptic argument use the parameter m [3]. The reasons are the following. First, using the modular angle reduce values of k to 1 1 k − ≤ ≤ and values of m to 0 1 m ≤ ≤ . To obtain the values of an elliptic function for k and m outside these intervals, we must use reciprocal-modulus transformation formulas [1,3]. Second, any elliptic integral and any Jacobian elliptic function is symmetric with respect to a real modulus k, that is, However, if the modulus is a pure imaginary number the symmetry does not hold. In this case, ik k m = − = , i.e. the parameter m is real. Therefore we have Using the parameter m thus covers all real moduli and all pure imaginary moduli. Further, if we have a function ( ) | Q x m then we can define Besides these three elliptic functions, there are also nine Glaisher elliptic functions which are defined as quotients of these functions [4,5]. In the collection, all Jacobian elliptic functions are calculated by procedure sncndn [17,18]. The procedure is based on Gauss transformation and solves ( ) The inverse of Jacobian elliptic functions are but a special case of the elliptic integral of the first kind [3].

Collection description
The functions in the collection are available in two levels: • low-level functions with scalar arguments • higher level functions with matrix arguments It is assumed that the input arguments of low-level functions are real scalars without check. Actual computation is conducted by the low-level functions which as input has the parameter m (mfunctions). Most of these computational functions are based on Matlab translations of Algol procedures from [17,19,20], Fortran functions from [18,21], and Pascal procedures from [22]. More precisely, for computation of elliptic integrals and inverse of Jacobian elliptic functions either Bulirsch's integrals el1, el2, el3 or Carlson's integrals rc, rd, rg, rj, are used [1]. The core function for computation of Jacobian elliptic function is procedure sncndn from [1]. Most of the low-level functions which use the module k as input argument are wrappers, i.e. the functions which calls appropriate lowlevel m-function by setting 2 m k = .
All high-level functions are wrappers which mimics an elemental behavior of a function by calling either of the functions ufun1, ufun2, ufun3, ufun4. Here the term elemental is browed from Fortran, i.e. it means that high level functions may be called with matrix arguments of the same size (or any of them can be scalar) in which case an coresponded low-level function is applied element-wise, with a conforming matrix return value. All higher level functions check its input data; matrix class and a number of arguments.
Naming convention: All high-level functions use Maple style naming, i.e. the functions begin with an uppercase letter. In this way, they cannot be a mismatch with Matlab symbolic toolbox functions which begin with a lowercase letter. Each function is available either with the modulus k or parameter m as an argument. In the later case, the function name start with letter 'm'. Incomplete elliptic integrals are given either in Jacobi form, Legendre form or Jacobi's second form. Legendre form begins with letter 'p' (from the argument which is 'phi'). For the Jacobi's second form of elliptic integrals we use special names: JacobiEpsilon and JacobiLambda.
We enhance the above description with examples. If the modulus k is used as an argument then the low-level name for Jacobian elliptic function sn is 'jsn' and its higher level name is 'JacobiSN'. If the parameter m is used as argument then the names are 'mjsn' and 'mJacobiSN'. Low-level name for a function which calculates the Jacobi form of elliptic integral of first kind with k as argument is 'elF' and the name of the higher level function is EllipticF. If the parameter m is used as an argument then the names are 'melF' and mEllipticF. If Legendre form of integral is used then the names are either 'pelF' and pEllipticF or 'mpelF' and mpEllipticF.
The names of some function are not standard. Thus for convenience, we named ( ) as JacobiLambda and JacobiOmega, trough Jacobi does not use these functions [13].

Robustness
In the present collection, the number and type of input data are checked in higher-level routines. Lowlevel routines intercept inputs which are NaN and out of function domain.
All computational routines were tested for very small and for huge arguments. Some additional checks are added to the function which use successive Landen's transformation to prevent them from hanging if any variable become NaN.

Efficiency
As was already stated the functions in the collection mimics the elemental-behavior on the matrix arguments i.e. the code is not vectorized. This has an impact to functions efficiency if a vector or matrix arguments are used. From Table 1-2 we can see that vectorised MATLAB function are on matrix argument two to five times faster than present functions. However, this is not critical because as we can see from Table 3 the computation time for a million of function evolutions is a matter of half of the second. Also from the Table1-3, we can also see that the present functions are faster from Matlab functions in the case when a function is called with a scalar argument. As expected the present functions are much faster than elliptic functions from Matlab symbolic toolbox.

Accuracy
The

Examples
As simple application of present formulas consider Euler's flexural elastica which has parametric form [23,24] ( ) ( where ω is load parameter and C is constant. The following program produces the shapes of the elstica shown on  As a second example, we consider finite-strain elastic cantilever under follower force from [24]. The coordinates X, Y of deformed cantilever base curve are and ( ) The non-dimensional parameters in the above equations are the load parameter ω, the generalized slenderness ratio λ and the stiffness ratio. For cantilever under follower force we have sn , 2sin 1 cn , where 1 ψ is given the clockwise angle between the direction of force and inward normal to the cantilever deformed cross section at the free end. When the cantilever is shearless i.e. when 1 ν = , its deformed length L is given by The program which implements these formulas is given in Appendix B. The shape of deformed cantilever shown in   Table 1 [24].

Conclusion
The described collection contains a complete set of Matlab functions for the calculation of real elliptic integrals and real Jacobian elliptic functions in their all possible forms. The functions were tested for correctness, robustness, efficiency, and accuracy. The collection is freely available from https://www.mathworks.com/matlabcentral/fileexchange/65915-elfun18 . The reference manual with the examples is available from https://www.researchgate.net/publication/323399074_elfun18_-_Elliptic_Integrals_Jacobian_Elliptic_Functions_and_Theta_Functions_-_Reference_Manual_v11.